Entendamos como un sistema se comporta
¿Cuál es la respuesta temporal?
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}$$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}$$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$$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)$$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$$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)}$$Lo primero que debemos hacer es entender por que un sistema explota o no.
Recordemos que los objetivos en control son:
con $a>0$
t = numpy.linspace(0,5,num=100)
y = numpy.exp(1*t)
plt.plot(t,y);
con $a<0$
t = numpy.linspace(0,5,num=100)
y = numpy.exp(-1*t)
plt.plot(t,y);
con $a=0$
t = numpy.linspace(0,5,num=100)
y = numpy.exp(0*t)
plt.plot(t,y);
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} }$$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)
Para la siguiente matriz
$$A = \left[\array{1&0\\0&-1}\right]$$encontrar los valores y vectores propios ⚠️
Encontremos las matrices para un péndulo regular y un péndulo invertido. Sin fricción. ⚠️
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?
Analicemos la estabilidad de un enjambre de robots para resolver el problema de Rendezvous
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$$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} \}$$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ó.
display(ani)
Si hay más de dos robots, deberiamos pensar en llevarlos a todos al centroide de sus vecinos (o algo parecido)
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.
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}$$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?
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)$$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)
Analizando los valores propios del sistema tenemos que:
display(mA.eigenvals())
{-I: 1, I: 1}
El sistema es criticamente estable.
Con el controlador propuesto, la particula se comporta como se muestra.
display(ani)
¿Por qué no se queda en el origen si este es el objetivo?
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.
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)
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}
Tomemos los valores para el controlador $k_1=1$ y $k_2=2$:
display(ani)
Es claro que algunos valores propios son mejores que otros.
Próximamente veremos como selecionar los valores propios.