Estabilidad de sistemas LTI¶

Entendamos como un sistema se comporta

Dinámica del sistema¶

  • Empecemos ignorando la entrada del sistema, el estado será $x$
$$\array{\dot{x}=A\, x &&& x(t_0)=x_0}$$

¿Cuál es la respuesta temporal?

Resolviendo la ecuación diferencial ordinaria¶

Si todo fuera escalar:

$$\array{\dot{x}=a\, x &,& x(t_0)=x_0 & \to & x(t)=e^{a(t-t_0))}x_0 }$$

Como sabemos esto? verificamos la posible solución:

$$\array{ x(t_0)=e^{a(t_0-t_0))}x_0=e^{0}x_0=x_0 & \text{condición inicial verificada}\\ \frac{d}{dt}x(t_0)=ae^{a(t-t_0))}x_0=ax & \text{dinámica verificada} } $$

Para sistemas de orden superior, tenemos una versión matricial de esto:

$$\array{\dot{x}=A\, x &,& x(t_0)=x_0 &\to& x(t)=e^{A(t-t_0))}x_0}$$

Exponencial de una matriz¶

La definición es similar a la definición de escalares:

$$e^{at}=\sum_{k=0}^\infty\frac{a^kt^k}{k!} \qquad\qquad e^{At}=\sum_{k=0}^\infty\frac{A^kt^k}{k!}$$

cuya derivada es:

$$\frac{d}{dt} \sum_{k=0}^\infty\frac{A^kt^k}{k!} = 0 + \sum_{k=0}^\infty\frac{kA^kt^{k-1}}{k!} = A\sum_{k=0}^\infty\frac{kA^{k-1}t^{k-1}}{(k-1)!} = A\sum_{k=0}^\infty\frac{A^kt^k}{k!} $$

es decir que:

$$\frac{d}{dt}e^{At}=Ae^{At}$$

Resolviendo la ecuación del controlador¶

La exponencial de una matriz tiene un papel fundamental y por eso tiene su propio nombre.

Matriz de transición de estado

$$e^{A(t-t_0))}=\Phi(t,t_0)$$

para una dinámica de estado $\dot{x}= A x$ la solución será $x(t)=\Phi(t,\tau)x(\tau)$, tiene las siguientes propiedades:

$$\cases{\frac{d}{dt}\Phi(t,t_0)=A\Phi(t,t_0) \\ \Phi(t,t)=I}$$

Cuando tenemos un sistema controlado, tendremos:

$$x(t)= \Phi(t,t_0)x(t_0)+\int_{t_0}^{t}\Phi(t,\tau)B\,u(\tau)d\tau$$

Resolviendo la ecuación del controlador¶

Verifiquemos que la ecuación cumple con la condición inicial:

$$x(t_0)= \Phi(t_0,t_0)x(t_0)+\int_{t_0}^{t_0}\Phi(t_0,\tau)B\,u(\tau)d\tau = I x(t_0) + 0 = x(t_0)$$
$$x(t)= \Phi(t,t_0)x(t_0)+\int_{t_0}^{t}\Phi(t,\tau)B\,u(\tau)d\tau$$

y cumple para la derivada:

$$\frac{d}{dt}x(t)= A\Phi(t,t_0)x(t_0)+\frac{d}{dt}\int_{t_0}^{t}\Phi(t,\tau)B\,u(\tau)d\tau$$$$\frac{d}{dt}\int_{t_0}^{t}f(t,\tau)d\tau = f(t,t) +\int_{t_0}^{t} \frac{d}{dt}f(t,\tau)d\tau = \Phi(t,t)B\,u(t) + \int_{t_0}^{t} A\Phi(t,\tau)B\,u(\tau)d\tau$$$$\frac{d}{dt}x(t)= A\Phi(t,t_0)x(t_0)+B\,u(t) + \int_{t_0}^{t} A\Phi(t,\tau)B\,u(\tau)d\tau$$$$\frac{d}{dt}x = Ax + B u$$

En resumen¶

tenemos que

$$\dot{x}=A\,x+B\,u \qquad y \qquad y=C\,x$$

luego,

$$y(t) = C \Phi(t,t_0)x(t_0) + C \int_{t_0}^{t}\Phi(t,\tau)B\,u(\tau)d\tau$$

sabiendo que

$$\Phi(t,\tau)= e^{A(t-\tau)}$$

Estabilidad¶

Lo primero que debemos hacer es entender por que un sistema explota o no.

Recordemos que los objetivos en control son:

  • Ser estable
  • Seguir el objetivo
  • Ser robusto
  • Entre otros

Estabilidad en sistemas escalares¶

  • Es útil empezar con sistemas escalares para ganar intuición.
$$\dot{x}=a\,x \quad \to\quad x(t)= e^{at}x(0)$$

con $a>0$

t = numpy.linspace(0,5,num=100)
y = numpy.exp(1*t)
plt.plot(t,y);

Estabilidad en sistemas escalares¶

  • Es útil empezar con sistemas escalares para ganar intuición.
$$\dot{x}=a\,x \quad \to\quad x(t)= e^{at}x(0)$$

con $a<0$

t = numpy.linspace(0,5,num=100)
y = numpy.exp(-1*t)
plt.plot(t,y);

Estabilidad en sistemas escalares¶

  • Es útil empezar con sistemas escalares para ganar intuición.
$$\dot{x}=a\,x \quad \to\quad x(t)= e^{at}x(0)$$

con $a=0$

t = numpy.linspace(0,5,num=100)
y = numpy.exp(0*t)
plt.plot(t,y);

Tres casos¶

  • Asintóticamente estable: $x(t)\to0, \forall x(0)$
  • Inestable: $\exists x(0) : ||{x(t)}||\to\infty$
  • Críticamente estable: en medio, no explota pero tampoco llega a cero

De

$$\dot{x}=a\,x \quad \to\quad x(t)= e^{at}x(0)$$

tenemos entonces:

$$\cases{ a>0 \quad: & \text{inestable} \\ a<0 \quad: & \text{asintóticamente estable}\\ a=0 \quad: & \text{críticamente estable} }$$

Estabilidad en Matrices¶

$$\dot{x}=A\,x \quad \to\quad x(t)= e^{At}x(0)$$

No podemos decir que $A>0$, pero podemos usar los valores propios.

$$A\,v = \lambda\,v$$

donde $\lambda \in \mathscr{C}$ son los valores propios y $v \in \mathscr{R}^n$ son los vectores propios.

Los valores propios no diran como la matriz $A$ actua en cada dirección.

En MATLAB

>> eig(A)

Ejemplo¶

Para la siguiente matriz

$$A = \left[\array{1&0\\0&-1}\right]$$

encontrar los valores y vectores propios ⚠️

$$\array{\lambda_1=1 &,& \lambda_2=-1 &,& v_1=\left[\array{1\\0}\right] &,& v_1=\left[\array{0\\1}\right]}$$

Estabilidad de Matrices¶

$$\dot{x}=A\,x \quad \to\quad x(t)= e^{At}x(0)$$
  • Asintóticamente estable (si y solo si): $$Re(\lambda)<0, \forall\lambda \in \texttt{eig}(A)$$
  • Inestable (si): $$\exists\lambda\in\texttt{eig}(A) : Re(\lambda)>0$$
  • Críticamente estable (solo si): $$Re(\lambda)\le0, \forall\lambda \in \texttt{eig}(A)$$
  • Críticamente estable (si): Un valor propio es cero y el resto tienen parte real negativa o Dos valores propios son puramente imaginarios y el resto tienen parte real negativa.

El cuento de dos pendulos¶

Encontremos las matrices para un péndulo regular y un péndulo invertido. Sin fricción. ⚠️

$$A_{\text{péndulo regular}} = \left[\array{0&1\\-1&0}\right] \qquad A_{\text{péndulo invertido}} = \left[\array{0&1\\1&0}\right]$$

encontremos los valores propios ⚠️

  • Para el péndulo regular tenemos: $$\lambda_1 = j \qquad \lambda_2=-j$$

  • Para el péndulo invertido tenemos: $$\lambda_1 = -1 \qquad \lambda_2=1$$

¿Sistemas estables?

Enjambre de robots¶

Analicemos la estabilidad de un enjambre de robots para resolver el problema de Rendezvous

  • Tenemos una colección de robots que pueden medir su posición relativa a sus vecinos.
  • Problema: hacer que todos los robots se encuentren en el mismo lugar (no especificada).

Caso de dos robots simples¶

Anteriormente revisamos el caso de dos robots en una linea.

Si los robots se dirigen el uno hacia el otro tenemos:

$$\cases{u_1=x_2-x_1\\u_2=x_1-x_2}$$

La matriz dinámica sera entonces:

$$\dot{x}=\left[\array{-1&1\\1&-1}\right]x$$

Estabilidad en el caso de dos robots simples¶

$$A = \left[\array{-1&1\\1&-1}\right] \qquad \lambda_1 = 0 \qquad \lambda_2=-2$$

En este sistema tenemos un valor propio 0 y todos los demas tienen parte real negativa. Aqui el estado del sistema terminara en algo llamado el *Espacio nulo (null-space)* de $A$:

$$null(A)=\{x: Ax = 0\}$$

para el caso particular de esta $A$, el espacio nulo es:

$$null(A)=\{x: x = \left[\array{\alpha\\\alpha}\right] , \alpha \in \mathscr{R} \}$$

Estabilidad en el caso de dos robots simples¶

Si $x_1\to\alpha$ y $x_2\to\alpha$ entonces tenemos que $(x_1-x_2)\to0$ por lo que el Rendezvous se logró.

# https://www.geeksforgeeks.org/how-to-install-ffmpeg-on-windows/ 
display(ani)
Your browser does not support the video tag.

Si hay más de dos robots, deberiamos pensar en llevarlos a todos al centroide de sus vecinos (o algo parecido)

Realimentación del sistema en espacio de estados¶

Sabemos que lo primero que debemos hacer para controlar un sistema es llevarlo a ser asintóticamente estable. Es decir que la parte real de todos sus valores propios sea negativa.

Una particula en una linea recta sin fricción¶

Diseñemos un control proporcional para el sistema propuesto, utilizando el espacio de estados.

$$m \ddot{x} = F$$

La ecuación de estado quedaría asi (con $m=1$).

$$\dot{\mathbf{x}}=\left[\begin{array}{}0&1\\0&0\end{array}\right]\mathbf{x}+\left[\begin{array}{}0\\1\end{array}\right]\mathbf{u}$$

La ecuación de salida sería.

$$\mathbf{y}=\left[\begin{array}{}1&0\end{array}\right]\mathbf{x}$$

Control de la particula¶

Si queremos controlar el sistema con un lazo cerrado debemos conectar de alguna forma la salida $\mathbf{y}$ con la entrada $\mathbf{u}$

plt.plot([-2,2],[0,0])
plt.plot([0,0],[-1,1],color='r')
plt.plot([-1],[0],color='k',marker='.',markersize=30)
plt.xlim(-2,2);

El objetivo de control es que se mueva al origen. ¿Cómo lo logramos?

$$\left.\begin{array}{}u>0 \text{ si } y<0 \\ u<0 \text{ si } y>0\end{array}\right\}\qquad \to \qquad u=-y$$

Modificación de la dinamica del sistema¶

En general tenemos:

$$\mathbf{u}=-K\mathbf{y} = -KC\mathbf{x} $$

entonces:

$$\dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{u}=A\mathbf{x}-BKC\mathbf{x} = \left(A-BKC\right)\mathbf{x}$$

tenemos aquí un nuevo sistema, el sistema en lazo cerrado.

$$\dot{\mathbf{x}} = \left(A-BKC\right)\mathbf{x}=\hat{A}\mathbf{x}$$

nuestro trabajo ahora es seleccionar $K$ de tal forma que los valores propios de la matriz $\hat{A}$ den por lo menos un sistema estable.

$$Re(\lambda)<0 \forall\lambda \in \texttt{eig}(A-BKC)$$

La matriz de estado del sistema en lazo cerrado¶

Remplazamos los valores de las matrices y de $K=1$:

$$\hat{A}=\left(A-BKC\right)=\left(\left[\begin{array}{}0&1\\0&0\end{array}\right]-\left[\begin{array}{}0\\1\end{array}\right]1\left[\begin{array}{}1&0\end{array}\right]\right)$$ 
mA = sympy.Matrix([[0,1],[0,0]])-sympy.Matrix([0,1])*1*sympy.Matrix([[1,0]])
display(mA)
$\displaystyle \left[\begin{matrix}0 & 1\\-1 & 0\end{matrix}\right]$

Analizando los valores propios del sistema tenemos que:

display(mA.eigenvals())

{-I: 1, I: 1}

El sistema es criticamente estable.

Respuesta de la particula al lazo cerrado¶

Con el controlador propuesto, la particula se comporta como se muestra.

display(ani)
Your browser does not support the video tag.

¿Por qué no se queda en el origen si este es el objetivo?

  • No tenemos en cuenta la velocidad.
  • Necesitamos la información del estado (posición,velocidad) del sistema para estabilizarlo.

Estabilizando la particula¶

Para estabilizar la particula necesitamos conocer la información de todos los estados del sistema. Salvo que en nuestro sistema solo tenemos un sensor de posición:

$$\mathbf{y}=\left[\begin{array}{}1&0\end{array}\right]\mathbf{x}$$

El estado desconocido es la velocidad, el cual puede ser estimado de la posición. Por ahora supongamos que podemos medir ambos estados.

$$\mathbf{y}_{supuesto}=\left[\begin{array}{}1&0\\0&1\end{array}\right]\mathbf{x}$$

con esta nueva matriz $C$, proponemos un controlador $K$:

$$K = \left[\begin{array}{}k_1&k_2\end{array}\right]$$

Encontremos la nueva matriz $\hat{A}$ para el sistema en lazo cerrado.

Nueva matriz de estado del sistema en lazo cerrado¶

Recordemos que aquí estamos suponiendo que podemos medir ambos estados del sistema (posición y velocidad). Remplazamos los valores de las matrices y de $K$:

$$\hat{A}=\left(A-BKC\right)=\left(\left[\begin{array}{}0&1\\0&0\end{array}\right]-\left[\begin{array}{}0\\1\end{array}\right]\left[\begin{array}{}k_1&k_2\end{array}\right]\left[\begin{array}{}1&0\\0&1\end{array}\right]\right)$$

Luego $\hat{A}=$

k1,k2 = sympy.symbols('k_1 k_2')
mA3 = sympy.Matrix([[0,1],[0,0]])-sympy.Matrix([0,1])*sympy.Matrix([[k1,k2]])*sympy.Matrix([[1,0],[0,1]])
display(mA3)
$\displaystyle \left[\begin{matrix}0 & 1\\- k_{1} & - k_{2}\end{matrix}\right]$

Los valores propios o polos del sistema son:

eigen = mA3.eigenvals()
display(eigen)

{-k_2/2 - sqrt(-4*k_1 + k_2**2)/2: 1, -k_2/2 + sqrt(-4*k_1 + k_2**2)/2: 1}

Verifiquemos el comportamiento de la particula¶

  • Tomemos los valores para el controlador $k_1=1$ y $k_2=1$
  • Tomemos los valores para el controlador $k_1=0.1$ y $k_2=1$

Verifiquemos el comportamiento de la particula¶

Tomemos los valores para el controlador $k_1=1$ y $k_2=2$:

display(ani)
Your browser does not support the video tag.

Los valores propios importante¶

  • Es claro que algunos valores propios son mejores que otros.

    • Algunos causan oscilaciones
    • Algunos hacen que el sistema responda lentamente
    • etc.

  • Próximamente veremos como selecionar los valores propios.