BE3024 - Sistemas de Control 1 (Biomédica)

2do ciclo, 2024

Lección 11: Control de sistemas LTI multivariable

¿Cómo se ve el control en el espacio de estados?

Un ejemplo clave

modelo mínimo de Bergman

glucemia \(y(t)\)

tasa de ingreso de insulina

ingesta de glucosa

control PID

glucemia basal saludable \(r(t)\)

+
-

¿Qué ocurre con la glucemia conforme avanza el tiempo?

y(t)\to r(t)

la salida rastrea a la referencia

y(t)\to r(t)

la salida rastrea a la referencia

Como ahora sabemos, esto no corresponde a la perspectiva completa.

y(t)\to r(t)

la salida rastrea a la referencia

Como ahora sabemos, esto no corresponde a la perspectiva completa.

 

¿Qué ocurre con el estado del sistema?

modelo mínimo de Bergman

glucemia \(y(t)\)

tasa de ingreso de insulina

ingesta de glucosa

control PID

glucemia basal saludable \(r(t)\)

+
-

¿Qué ocurre internamente en el sistema?

¿Cómo evolucionan las variables de estado conforme la salida rastrea a la referencia?

\begin{aligned} \dot{x}_1(t) = & -p_1(x_1(t)-G_b)-x_1(t)x_2(t)+D(t)+w(t)\\ \dot{x}_2(t) = & -p_2x_2(t)+p_3(x_3(t)-I_b)\\ \dot{x}_3(t) = & -n(x_3(t)-I_b)+\gamma\max(0, x_1(t)-\sigma)t+u(t) \end{aligned}

¿Qué podemos observar?

x_1(t)
x_2(t)
x_3(t)

¿Qué podemos observar?

x_1(t)
x_2(t)
x_3(t)

condiciones iniciales distintas a cero

\mathbf{x}_0=\begin{bmatrix} x_1(0) \\ x_2(0) \\ x_3(0) \end{bmatrix}=\begin{bmatrix} 180 \\ 0 \\ 60 \end{bmatrix}

¿Qué podemos observar?

x_1(t)
x_2(t)
x_3(t)

Pareciera que \(y(t)=x_1(t)\), por lo que \(x_1(t)\) está rastreando a \(r(t)\).

¿Qué podemos observar?

x_1(t)
x_2(t)
x_3(t)

\(x_2(t)\) y \(x_3(t)\) también están llegando a cierto valor,

pero no se sabe cuál es o de dónde salió.

¿De qué otra forma podemos describir el comportamiento de las variables de estado sin hablar del valor al cual están llegando?

x_1(t)
x_2(t)
x_3(t)
x_1(t)
x_2(t)
x_3(t)

¿De qué otra forma podemos describir el comportamiento de las variables de estado sin hablar del valor al cual están llegando?

x_1(t)
x_2(t)
x_3(t)
\Rightarrow \quad \begin{matrix} \dot{x}_1 = 0 \\ \dot{x}_2 = 0 \\ \dot{x}_3 = 0 \end{matrix}

¿De qué otra forma podemos describir el comportamiento de las variables de estado sin hablar del valor al cual están llegando?

¿De qué otra forma podemos describir el comportamiento de las variables de estado sin hablar del valor al cual están llegando?

x_1(t)
x_2(t)
x_3(t)
\Rightarrow \quad \begin{matrix} \dot{x}_1 = 0 \\ \dot{x}_2 = 0 \\ \dot{x}_3 = 0 \end{matrix}

Las variables de estado dejan de evolucionar

\(\equiv\) las variables de estado se estabilizan

\(\equiv\) se estabiliza el estado.

En el espacio de estados, las consecuencias del control son la estabilización del estado y el rastreo de referencias en la salida.

Parte I: Estabilización de estado

En esencia, la estabilización de estado implica

\displaystyle \lim_{t\to\infty} \dot{\mathbf{x}}(t) = \mathbf{0}

lo cual se garantiza al hacer que el sistema

\begin{cases} \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u}\end{cases}

sea asintóticamente estable.

En esencia, la estabilización de estado implica

\displaystyle \lim_{t\to\infty} \dot{\mathbf{x}}(t) = \mathbf{0}

lo cual se garantiza al hacer que el sistema

\begin{cases} \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u}\end{cases}

sea asintóticamente estable.

Todos los eigenvalores de \(\mathbf{A}\) deben estar en el lado izquierdo del plano complejo.

La estrategia principal en el espacio de estados implica utilizar retroalimentación de estado.

\mathbf{u}=-\mathbf{K}\mathbf{x}, \quad \mathbf{K}\in\mathbb{R}^{m \times n}

Esto transforma la ecuación de estado.

\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\left(-\mathbf{K}\mathbf{x}\right) =\mathbf{A}\mathbf{x}-\mathbf{B}\mathbf{K}\mathbf{x} \\ \dot{\mathbf{x}}=\left(\mathbf{A}-\mathbf{B}\mathbf{K}\right)\mathbf{x}=\mathbf{A}_{c\ell}\mathbf{x}

La estrategia principal en el espacio de estados implica utilizar retroalimentación de estado.

\mathbf{u}=-\mathbf{K}\mathbf{x}, \quad \mathbf{K}\in\mathbb{R}^{m \times n}

Esto transforma la ecuación de estado.

\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\left(-\mathbf{K}\mathbf{x}\right) =\mathbf{A}\mathbf{x}-\mathbf{B}\mathbf{K}\mathbf{x} \\ \dot{\mathbf{x}}=\left(\mathbf{A}-\mathbf{B}\mathbf{K}\right)\mathbf{x}=\mathbf{A}_{c\ell}\mathbf{x}

matriz de constantes

Para sistemas LTI, la estabilización de estado se reduce a seleccionar los valores de \(\mathbf{K}\) tal que todos los eigenvalores de \(\mathbf{A}_{c\ell}\) se encuentren en el lado izquierdo del plano complejo.

De forma manual mediante pole placement.

Dos técnicas para estabilización

De forma automática mediante el Regulador Lineal Cuadrático (LQR).

sys = ss(A, B, C, D);
K = place(sys.A, sys.B, p)
K = lqr(sys.A, sys.B, Q, R)

Dos técnicas para estabilización

De forma manual mediante pole placement.

De forma automática mediante el Regulador Lineal Cuadrático (LQR).

sys = ss(A, B, C, D);
K = place(sys.A, sys.B, p)
K = lqr(sys.A, sys.B, Q, R)

CUIDADO, esto funcionará si y sólo si el sistema es completamente controlable (C.C.).

Un sistema LTI es completamente controlable si y sólo si su matriz de controlabilidad

Controlabilidad

tiene rango

Gamma = ctrb(sys)
rank(Gamma)

igual al orden (número de variables de estado = \(n\)) del sistema.

Ambas técnicas emplean la siguiente arquitectura:

Regresando a estabilización

sistema

¿Cómo pueden seleccionarse estos polos deseados? Mediante regiones de diseño.

Pole placement

Nota: el algoritmo no permite que dos polos sean exactamente iguales. Ej: {-1, -1} vs {-1, -1.01}.

K = place(sys.A, sys.B, p)

vector de \(n\) polos deseados

\mathrm{Re}(s)
\mathrm{Im}(s)
\sigma
\omega_d
-\omega_d
\omega_n
\theta
\theta=\sin^{-1}(\zeta)
\mathrm{Re}(s)
\mathrm{Im}(s)
\omega_n \ge \dfrac{1.8}{t_r}
\omega_n

Región donde se cumple con los requerimientos.

t_r \le \text{spec.}
\mathrm{Re}(s)
\mathrm{Im}(s)
\zeta \ge f^{-1}\left(M_p\right)
\theta
\theta=\sin^{-1}(\zeta)
M_p \le \text{spec.}

Región donde se cumple con los requerimientos.

\mathrm{Re}(s)
\mathrm{Im}(s)
\sigma \ge \dfrac{4.6}{t_s}
\sigma
(\epsilon=1\%)
t_s \le \text{spec.}

Región donde se cumple con los requerimientos.

Ejemplo: pole placement

Estabilice el sistema, garantizando un tiempo de subida menor a 1 segundo y un overshoot menor al 10%.

Ejemplo: pole placement

Listado de tareas:

  1. Encontrar el modelo del sistema.
  2. Verificar estabilidad.
  3. Verificar controlabilidad.
  4. Encontrar polos deseados.
  5. Encontrar \(\mathbf{K}\) mediante pole placement.
  6. Crear un nuevo sistema con el control aplicado.
  7. Evaluar el rendimiento (polos nuevos y respuesta al escalón).

>> clase11_ejemplo_poleplacement.m

Ejemplo: pole placement

Listado de tareas:

  1. Encontrar el modelo del sistema.
  2. Verificar estabilidad.
  3. Verificar controlabilidad.
  4. Encontrar polos deseados.
  5. Encontrar \(\mathbf{K}\) mediante pole placement.
  6. Crear un nuevo sistema con el control aplicado.
  7. Evaluar el rendimiento (polos nuevos y respuesta al escalón).

¿Qué ocurrió? Los polos se colocaron correctamente pero hay algo raro con el comportamiento en estado estacionario.

 

¿Qué pasa si se escogen otros polos?

Regulador Lineal Cuadrático (LQR)

Se recomienda iniciar con la matrices identidad e ir modificando los valores de la diagonal para penalizar el comportamiento no deseado en variables de estado y entradas/controles específicas.

K = lqr(sys.A, sys.B, Q, R)

matriz de penalización de estado \(\mathbf{Q}\in\mathbb{R}^{n \times n}\)

matriz de penalización de control \(\mathbf{R}\in\mathbb{R}^{m\times m}\)

Regulador Lineal Cuadrático (LQR)

A valores más altos, más se penaliza.

\mathbf{Q}=\begin{bmatrix} \square & 0 & \cdots & 0 \\ 0 & \square & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \square \end{bmatrix}
x_1
x_1
x_2
x_2
x_n
x_n
\cdots
\vdots

penalización para \(x_2\)

Regulador Lineal Cuadrático (LQR)

A valores más altos, más se penaliza.

\mathbf{R}=\begin{bmatrix} \square & 0 & \cdots & 0 \\ 0 & \square & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \square \end{bmatrix}
u_1
u_1
u_2
u_2
u_m
u_m
\cdots
\vdots

penalización para \(u_1\)

Ejemplo: LQR

Estabilice mediante el regulador lineal cuadrático.

Ejemplo: LQR

Listado de tareas:

  1. Encontrar el modelo del sistema.
  2. Verificar estabilidad.
  3. Verificar controlabilidad.
  4. Establecer las matrices de penalización.
  5. Encontrar \(\mathbf{K}\) mediante LQR.
  6. Crear un nuevo sistema con el control aplicado.
  7. Evaluar el rendimiento.

>> clase11_ejemplo_lqr.m

Ejemplo: LQR

Listado de tareas:

  1. Encontrar el modelo del sistema.
  2. Verificar estabilidad.
  3. Verificar controlabilidad.
  4. Establecer las matrices de penalización.
  5. Encontrar \(\mathbf{K}\) mediante LQR.
  6. Crear un nuevo sistema con el control aplicado.
  7. Evaluar el rendimiento.

De nuevo, la respuesta en estado estable se comportó extraño.

Conclusión: debe considerarse el rastreo además de la estabilización.

Nota: en el laboratorio aprenderemos un "truco" que puede servir para ajustar el error en estado estacionario (en ciertos casos), sin usar el rastreo explícitamente como veremos más adelante.

Implementación

Regresemos a la primera solución para el ejemplo del circuito.

K = [4, 3]

¿Qué implica esto?

\mathbf{u}=-\mathbf{K}\mathbf{x} \quad \Rightarrow u=-\begin{bmatrix} 4 & 3 \end{bmatrix}\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}

Una simple combinación lineal, versus una ecuación de diferencias "complicada" para el PID.

u=-4x_1-3x_2=-4v_c-3\imath_L

Implementación

u=-4x_1-3x_2=-4v_c-3\imath_L

Esto, sin embargo, asume que tenemos acceso a las variables de estado del sistema (!!!)

Implementación

Una simple combinación lineal, versus una ecuación de diferencias "complicada" para el PID.

Parte II: Rastreo de referencias

\begin{cases} \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u}\end{cases}

Fuimos capaces de estabilizar el estado

\displaystyle \lim_{t\to\infty} \dot{\mathbf{x}}(t) = \mathbf{0}
\begin{cases} \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u}\end{cases}

Fuimos capaces de estabilizar el estado

\displaystyle \lim_{t\to\infty} \dot{\mathbf{x}}(t) = \mathbf{0}
\begin{cases} \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u}\end{cases}

empleando retroalimentación de estado.

-\mathbf{K}\mathbf{x}

Fuimos capaces de estabilizar el estado

\displaystyle \lim_{t\to\infty} \dot{\mathbf{x}}(t) = \mathbf{0}

¿Qué hace falta?

\begin{cases} \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u}\end{cases}

empleando retroalimentación de estado.

-\mathbf{K}\mathbf{x}

Fuimos capaces de estabilizar el estado

\displaystyle \lim_{t\to\infty} \dot{\mathbf{x}}(t) = \mathbf{0}
\begin{cases} \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u}\end{cases}

Que la salida rastree a la referencia.

empleando retroalimentación de estado.

-\mathbf{K}\mathbf{x}

Para lograr esto se requiere estabilizar el error

\dot{\mathbf{e}} = \mathbf{r}-\mathbf{y}=\mathbf{r}-\mathbf{C}\mathbf{x}-\mathbf{D}\mathbf{u}

Es decir, el rastreo de referencias se hace mediante la estabilización del error.

\displaystyle \lim_{t\to\infty} \dot{\mathbf{e}}(t) = \mathbf{0} \triangleq \mathbf{r}-\mathbf{y}, \quad \mathbf{r}\in\mathbb{R}^p

Para lograr esto se requiere estabilizar el error

\displaystyle \lim_{t\to\infty} \dot{\mathbf{e}}(t) = \mathbf{0} \triangleq \mathbf{r}-\mathbf{y}, \quad \mathbf{r}\in\mathbb{R}^p

Es decir, el rastreo de referencias se hace mediante la estabilización del error.

\(\Rightarrow\) dos problemas de estabilización

\dot{\mathbf{e}} = \mathbf{r}-\mathbf{y}=\mathbf{r}-\mathbf{C}\mathbf{x}-\mathbf{D}\mathbf{u}

Estabilizar estado + estabilizar error

\Rightarrow \begin{bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{e}} \end{bmatrix} =\begin{bmatrix} \mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ -\mathbf{C}\mathbf{x}-\mathbf{D}\mathbf{u}+\mathbf{r} \end{bmatrix}

Estabilizar estado + estabilizar error

\Rightarrow \begin{bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{e}} \end{bmatrix} =\begin{bmatrix} \mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ -\mathbf{C}\mathbf{x}-\mathbf{D}\mathbf{u}+\mathbf{r} \end{bmatrix}
\begin{bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{e}} \end{bmatrix} =\begin{bmatrix} \mathbf{A} & \mathbf{0}_{n \times p} \\ -\mathbf{C} & \mathbf{0}_{p \times p} \end{bmatrix}\begin{bmatrix} \mathbf{x} \\ \mathbf{e} \end{bmatrix} +\begin{bmatrix} \mathbf{B} \\ -\mathbf{D} \end{bmatrix}\mathbf{u}+\begin{bmatrix} \mathbf{0}_{n \times p} \\ \mathbf{I}_p \end{bmatrix}\mathbf{r}

Estabilizar estado + estabilizar error

\Rightarrow \begin{bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{e}} \end{bmatrix} =\begin{bmatrix} \mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ -\mathbf{C}\mathbf{x}-\mathbf{D}\mathbf{u}+\mathbf{r} \end{bmatrix}
\begin{bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{e}} \end{bmatrix} =\begin{bmatrix} \mathbf{A} & \mathbf{0}_{n \times p} \\ -\mathbf{C} & \mathbf{0}_{p \times p} \end{bmatrix}\begin{bmatrix} \mathbf{x} \\ \mathbf{e} \end{bmatrix} +\begin{bmatrix} \mathbf{B} \\ -\mathbf{D} \end{bmatrix}\mathbf{u}+\begin{bmatrix} \mathbf{0}_{n \times p} \\ \mathbf{I}_p \end{bmatrix}\mathbf{r}
\dot{\mathbf{z}}
\mathbf{z}
\boldsymbol{\mathcal{A}}
\boldsymbol{\mathcal{B}}
\boldsymbol{\mathcal{M}}

Estabilizar estado + estabilizar error

\Rightarrow \begin{bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{e}} \end{bmatrix} =\begin{bmatrix} \mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \\ -\mathbf{C}\mathbf{x}-\mathbf{D}\mathbf{u}+\mathbf{r} \end{bmatrix}
\mathbf{y}=\begin{bmatrix} \mathbf{C} & \mathbf{0}_{p \times p} \end{bmatrix}\begin{bmatrix} \mathbf{x} \\ \mathbf{e} \end{bmatrix}
\mathbf{z}
\boldsymbol{\mathcal{C}}

Regulador Cuadrático Integral (LQI)

sys = ss(A, B, C, D);
sys_ext = ss(AA, BB, CC, D);

Regulador Lineal Cuadrático Integral (LQI)

\boldsymbol{\mathcal{A}}
\boldsymbol{\mathcal{B}}
sys = ss(A, B, C, D);
sys_ext = ss(AA, BB, CC, D);
\boldsymbol{\mathcal{C}}
KK = lqi(sys, QQ, R)
KK = lqr(sys_ext.A, sys_ext.B, QQ, R)
\boldsymbol{\mathcal{A}}
\boldsymbol{\mathcal{B}}
sys = ss(A, B, C, D);
sys_ext = ss(AA, BB, CC, D);
\boldsymbol{\mathcal{C}}

Regulador Lineal Cuadrático Integral (LQI)

KK = lqi(sys, QQ, R)
KK = lqr(sys_ext.A, sys_ext.B, QQ, R)
\boldsymbol{\mathcal{Q}}=\begin{bmatrix} \mathbf{Q} & \mathbf{0}_{n \times p} \\ \mathbf{0}_{p \times n} & \mathbf{Q}_y \end{bmatrix}
\boldsymbol{\mathcal{K}}=\begin{bmatrix} \mathbf{K} & \mathbf{K}_I \end{bmatrix}, \ \mathbf{K}_I \in \mathbb{R}^{m \times p}
\boldsymbol{\mathcal{A}}
\boldsymbol{\mathcal{B}}
sys = ss(A, B, C, D);
sys_ext = ss(AA, BB, CC, D);
\boldsymbol{\mathcal{C}}

Regulador Lineal Cuadrático Integral (LQI)

\mathbf{u}=-\mathbf{K}\mathbf{x}-\mathbf{K}_I\displaystyle\int_{t_0}^{t} \left[\mathbf{r}(\tau)-\mathbf{y}(\tau)\right] d\tau

i

\mathbf{u}=-\mathbf{K}\mathbf{x}-\mathbf{K}_I\displaystyle\int_{t_0}^{t} \left[\mathbf{r}(\tau)-\mathbf{y}(\tau)\right] d\tau

i

Acl = AA - BB*KK;
sys_cl = ss(Acl, MM, CC, [])
\mathbf{u}=-\mathbf{K}\mathbf{x}-\mathbf{K}_I\displaystyle\int_{t_0}^{t} \left[\mathbf{r}(\tau)-\mathbf{y}(\tau)\right] d\tau

i

Acl = AA - BB*KK;
sys_cl = ss(Acl, MM, CC, [])
\boldsymbol{\mathcal{M}}

>> clase11_ejemplo_lqi.m