Sergio C

Mar 192017
 

En esta entrada analizaremos el modelado fenomenológico de un reactor de Van de Vusse NO isotérmico, el cual será un proceso multivariable MIMO y será de gran ayuda para ejemplificar varios controladores Multivariables que trataremos en el Blog Control Automático Educación.

 

El modelado del Reactor de Van de Vusse No Isotérmico es básicamente el mismo modelado explicado en TOTAL DETALLE en la entrada pasada (Por favor dar Click Aquí para ver la entrada pasada si aún no la has visto) con la diferencia que vamos adicionar una tercera ecuación que representará el comportamiento de la temperatura del reactor, lo que nos permitirá ejercer control tanto en la concentración del producto B cuanto en la temperatura del reactor (Sistema MIMO). El esquema del proceso se muestra en la siguiente figura:

Reactor de Van de Vusse No Isotermico

A continuación son descritas las ecuaciones analíticas del balance de masa para los componentes A e B y el balance de energía en el reactor. Para esto es adoptada la siguiente nomenclatura:

  •  C_{a_0} concentración de A en la entrada del reactor;
  •  T_0 la temperatura de la carga;
  • C_a y C_b son las concentraciones de A y B en el reactor;
  • T es la temperatura del reactor;
  • F es el flujo de alimentación.

y es considerado:

  • que el reactor es perfectamente agitado;
  • que la densidad del líquido \rho y capacidad calorífica C_p constantes;
  • Flujos iguales a F en la entrada y en la salida del reactor;
  • reacción A\rightarrow B de 1a orden en relación a A;
  • reacción B\rightarrow C de 1a orden en relación a B;
  • reacción 2A\rightarrow D de 2a orden en relación a A;

 

Los parámetros de las tasas de reacción k_i,i = 1,2,3, dependen de la temperatura por medio de la ley de Arrhenius

 

K_i(T)=k_{i0} exp\left( \dfrac{-E_i/R}{T(C)+273.15} \right)

 

Donde E_i,i = 1,2,3, son las energías de activación de las tres diferentes reacciones que ocurren en el sistema y R es la constante universal de los gases.

 

  • En el flujo de entrada solo se considera el reagente A;
  • Se considera la dinámica de la camisa de refrigeración despreciable;
  • T_k temperatura de la camisa del reactor constante;
  • Calor transferido del reactor para la camisa como Q = K_w A_R (T - T_k), donde K_w es el coeficiente de transferencia de calor e A_R la área de superficie del reactor;
  • Las temperaturas de la reacción son llamados como (-\Delta H_{RAB}),(-\Delta H_{RBC}),(-\Delta H_{RAD}).

 

El proceso es descrito por el modelo no-adiabático representado por las ecuaciones de balance de masa por componente y energía en el reactor.

 

\dfrac{d(C_A)}{dt}=\dfrac{F}{V}(C_{Af}-C_A)-k_1(T)C_A-K_3(T)C_A^2

\dfrac{d(C_B)}{dt}=-\dfrac{F}{V}C_B+k_1(T)C_A-K_2(T)C_B

\dfrac{d(T)}{dt}=\dfrac{1}{\rho C_p} \left[k_1(T)C_a(-\Delta H_{RAB})+k_2(T)C_b(-\Delta H_{RBC})+k_3(T)C_a^2(-\Delta H_{RAD})\right]+\dfrac{F}{V} (T_0-T)+\dfrac{K_w A_R}{\rho C_p V} (T_k-T)

Los valores de los parámetros que describen este sistema fueron extraídos de (TRIERWEILER, 1997) e están descritos en la siguiente Tabla.

parametros reactor cstr

Aplicando el concepto de la Jacobiana, explicado en detalle en la entrada anterior (Click aquí para ver la entrada pasada) linealizamos el sistema sobre un punto de operación. El punto de operación escogido se da con la entrada de \dfrac{F}{V}=72h^{-1} y la entrada de temperatura T_k=128.95K.

 

Para este punto de Operación se llega a la siguiente matriz de funciones de transferencia MIMO 2×2:

P_n=\begin{bmatrix}  \dfrac{-1.1033 (s+26.52) (s+94.31)}{(s+162.2) (s^2 + 203s + 1.284e04)} & \dfrac{142.23 (s+31.3)}{(s+162.2) (s^2 + 203s + 1.284e04)}\\ \\  \dfrac{-6.7606 (s+159.2) (s-7.671)}{(s+162.2) (s^2 + 203s + 1.284e04)} & \dfrac{30.829 (s+131) (s+168.8)}{(s+162.2) (s^2 + 203s + 1.284e04)}  \end{bmatrix}

 

A continuación vemos una comparación entre el sistema Lineal y el Sistema NO Lineal haciendo una variación del tipo escalón igual a 5 sobre la variable de flujo, \dfrac{F}{V}=77h^{-1}.

Reactor CSTR Lineal vs no lineal

Reactor CSTR Lineal vs no lineal

A continuación se muestra el código en MATLAB que simula el modelo fenomenológico del reactor de Van de Vusse, realiza la linealización sobre el punto de operación y calcula la matriz de funciones de transferencia del sistema.

Archivo principal del programa, llamar como “mainVV.m”

 

Archivo función del modelo del Reactor, nombrar el archivo como “VanVusseModel.m”

Bibliografia

  • TRIERWEILER, J. (1997). A Systematic Approach to Control Structure Design, Tesis de Doctorado. Alemania: Universität Dortmund.
Mar 052017
 

Consideremos el siguiente modelo fenomenológico de un reactor de Van de Vusse, al cual se le desea obtener un modelo lineal en diferentes puntos de operación sobre el cual va a trabajar el proceso.

 

Descripción del proceso:

 

El CSTR (reactor perfectamente agitado, del ingles continuous stirred-tank reactor) de van de Vusse, mostrado en la figura a continuación, representa el sistema reacional utilizado para la síntesis de ciclopentenol (representado por B) a partir del ciclopentadieno (representado por A). Este reactor envuelve entonces reacciones en serie y en paralelo entre la entrada del sistema (Reacción A) y el producto obtenido (Reacción B). Sin embargo dado a la fuerte reatividade de los reagentes y produtos, se obtienen dos productos adicionales NO deseados, el diciclopentadieno (D) es produzido como subproduto por la reacción de Diels-Alder y ciclopentenodiol (C) es produzido como um produto consecutivo por la adición de otra molécula de água.

Reactor de Van de Vusse

A\overset{k_1}{\rightarrow}B\overset{k_2}{\rightarrow}C

2A\overset{k_3}{\rightarrow}D

donde los parámetros de la tasa de reacción, vienen dados por:

k_1=5/6 l/min

k_2=5/3 l/min

k_3=1/6 l/min

 

 

Balance General de Masa

Si se asume la densidad y el volumen constante en el interior del reactor tenemos:

\dfrac{dV}{dt}=0 y F\approx F_i

Balance del componente A:

Interior del reactor = Flujo de entrada – Flujo de Salida – sale reaccion 1 – sale reacción 2

\dfrac{d(VC_A)}{dt}=F(C_{Af}-C_A)-Vk_1C_A-VK_3C_A^2

Y debido a que consideramos el volumen constante, el balance de masa del componente A se resume a:

\dfrac{d(C_A)}{dt}=\dfrac{F}{V}(C_{Af}-C_A)-k_1C_A-K_3C_A^2

Balance del componente B:

De la misma forma podemos obtener el balance de masa para el componente B

\dfrac{d(C_B)}{dt}=-\dfrac{F}{V}C_B+k_1C_A-K_2C_B

Balance del componente C:

\dfrac{d(C_C)}{dt}=-\dfrac{F}{V}C_C+k_2C_B

Balance del componente D:

\dfrac{d(C_D)}{dt}=-\dfrac{F}{V}C_D+\dfrac{1}{2}k_3C_A^2

Como deseamos obtener el modelo del Componente B con relación al producto de entrada Componente A, observemos que estos dos balances de energia NO dependen de los balances C y D, por lo tanto, solo nos van a interesar estas dos ecuaciones diferenciales:

\dfrac{d(C_A)}{dt}=\dfrac{F}{V}(C_{Af}-C_A)-k_1C_A-K_3C_A^2

\dfrac{d(C_B)}{dt}=-\dfrac{F}{V}C_B+k_1C_A-K_2C_B

Estas dos ecuaciones representan el comportamiento fenomenológico no lineal de nuestro reactor de Van de Vusse. Ahora como deseamos linealizar este sistema en torno a un punto de equilibrio o punto de operación el próximo paso sera encontrar dicho punto. Para lograr esto hacemos todas las derivadas (cambios) iguales a cero, lo que significa que el sistema no está cambiando y se encuentra en un estado estacionario.

A partir de la ecuación de C_A obtenemos lo siguiente:

0=\dfrac{F_s}{V}(C_{Afs}-C_{As})-k_1C_{As}-K_3C_{As}^2

-K_3C_{As}^2+\left(-k_1-\dfrac{F_s}{V}\right)C_{As} +\dfrac{F_s}{V}C_{Afs}=0

Resolviendo la ecuación cuadrática anterior utilizando la formula general y escogiendo la raiz positiva, tenemos que:

C_{As}=\dfrac{-\left(k_1+\dfrac{F_s}{V}\right)}{2k_3}+\dfrac{\sqrt{\left(k_1+\dfrac{F_s}{V}\right)^2+4k_3\dfrac{F_s}{V}}}{2k_3}

Y para el componente B tenemos:

0=-\dfrac{F_s}{V}C_{Bs}+k_1C_{As}-K_2C_{Bs}

C_{Bs}=\dfrac{k_1C_{As}}{\dfrac{F_s}{V}+K_2}

La entrada de nuestro sistema es la tasa de disolución F/V l/min. Si hacemos un barrido de nuestra entrada desde 1 1/min hasta 10 1/min. Obtendremos las curvas que relacionan la entrada y la salida del sistema en su estado estacionario.

Van de Vusse Estacionario
Si observamos el comportamiento del Componente B, vemos que tiene una inversión de fase, que con una entrada de 1.292 1/min se obtiene la máxima producción del componente B (Lo más deseado) y que con una entrada superior la producción de B comienza a caer.

Linealización del sistema de ecuaciones no lineales en Espacio de Estados utilizando el Jacobiano

Vamos a obtener una representación en espacio de estados de nuestro reactor.

\dot{x}=Ax+Bu

y=Cx

Definimos cuales serán nuestras variables de estado:

x_1=C_A

\dot{x}_1=\dfrac{dC_A}{dt}

x_2=C_B

\dot{x}_2=\dfrac{dC_B}{dt}

u= F/V

Si renombramos las dos ecuaciones de balance de masa como función 1 (f_1) y función 2 (f_2).

\dfrac{d(C_A)}{dt}=f_1(C_A,C_B,F/V)=\dfrac{F}{V}(C_{Af}-C_A)-k_1C_A-K_3C_A^2

\dfrac{d(C_B)}{dt}=f_2(C_A,C_B,F/V)=-\dfrac{F}{V}C_B+k_1C_A-K_2C_B

Jacobiano:

A=\begin{bmatrix}  \dfrac{\partial f_1}{\partial x_1} & \dfrac{\partial f_1}{\partial x_2}\\ \\  \dfrac{\partial f_2}{\partial x_1} & \dfrac{\partial f_2}{\partial x_2}  \end{bmatrix},B=\begin{bmatrix}  \dfrac{\partial f_1}{\partial u} \\ \\  \dfrac{\partial f_2}{\partial u}  \end{bmatrix}

Matriz A:

A_{11}=\dfrac{\partial f_1}{\partial x_1}=-\dfrac{F_s}{V}-k_1-2k_3C_{As}

A_{12}=\dfrac{\partial f_1}{\partial x_2}=0

A_{21}=\dfrac{\partial f_2}{\partial x_1}=k_1

A_{22}=\dfrac{\partial f_2}{\partial x_2}=-\dfrac{F_s}{V}-k_2

Matriz B:

B_{11}=\dfrac{\partial f_1}{\partial u}=C_{Afs}-C_{As}

B_{21}=\dfrac{\partial f_2}{\partial u}=-C_{Bs}

El modelo en espacio de estado viene dado por:

\dot{x}=\begin{bmatrix}  -\dfrac{F_s}{V}-k_1-2k_3C_{As} &0\\ \\  k_1 & -\dfrac{F_s}{V}-k_2  \end{bmatrix}x+\begin{bmatrix}  C_{Afs}-C_{As} \\ \\  -C_{Bs}  \end{bmatrix}u

y=\begin{bmatrix}  0&1  \end{bmatrix}x

Teniendo la representación en el espacio de estado de nuestro sistema, vamos a seleccionar diferentes puntos de operación para analizar el comportamiento de este proceso. Veremos tres casos diferentes. En el caso 1 veremos como aparece una respuesta inversa o respuesta de fase no mínima en la concentración de B mientras que en el caso 2 esta respuesta inversa no aparece. Y en el caso 3 se ilustra el punto de operación donde se alcanza la máxima producción del componente B.

CASO 1

Para este punto de operación vamos a considerar la tasa de flujo \dfrac{F_s}{V}=\dfrac{4}{7}min^{-1}. Notemos que esta tasa de flujo se ubica en el lado izquierdo del máximo pico del gráfico de estado estacionario que relaciona la entrada salida.

Inicialmente encontramos el estado estable del componente A (C_{As}) y el componente B(C_{Bs}).

C_{As}=\dfrac{-\left(k_1+\dfrac{F_s}{V}\right)}{2k_3}+\dfrac{\sqrt{\left(k_1+\dfrac{F_s}{V}\right)^2+4k_3\dfrac{F_s}{V}}}{2k_3}

C_{As}=\dfrac{-\left(5/6+4/7\right)}{2(1/6)}+\dfrac{\sqrt{\left(5/6+4/7\right)^2+4(1/6)(4/7)}}{2(1/6)}

C_{As}=3\dfrac{mol}{litro}

C_{Bs}=\dfrac{k_1C_{As}}{\dfrac{F_s}{V}+K_2}

C_{Bs}=\dfrac{3(5/6)}{4/7+5/3}

C_{Bs}=1.117\dfrac{mol}{litro}

Con esto obtenemos la siguiente representación en el espacio de estado

\dot{x}=\begin{bmatrix}  -2.4048 &0\\ \\  0.8333 & -2.2381  \end{bmatrix}x+\begin{bmatrix}  7 \\ \\  -1.117  \end{bmatrix}u

y=\begin{bmatrix}  0&1  \end{bmatrix}x

Y podemos obtener la función de transferencia de la misma forma que lo hemos visto en el curso de sistemas dinámicos lineales usando G(s)=C(sI-A)^{-1}B.

G(s)=\dfrac{-1.117s+3.1472}{s^2+4.6429s+5.3821}

Si notamos, vemos que esta función de transferencia tiene un cero positivo, lo que provoca la respuesta inversa o respuesta de fase no mínima.

Cero=\dfrac{3.1472}{1.117}=2.8175

Al final del post, puedes encontrar el código de implementación en Matlab para que lo copies y lo pegues en tu computador.

 

Veamos una comparación de los dos sistemas en el punto de operación linealizado

Van de Vusse caso 1

Ahora una comparación alejándonos más del punto de operación linealizado. Veamos como el sistema lineal comienza a dejar de valer.

Van de Vusse caso 1

CASO 2

Para este punto de operación vamos a considerar la tasa de flujo \dfrac{F_s}{V}=2.8744min^{-1}. Notemos que esta tasa de flujo se ubica en el lado derecho del máximo pico del grafico de estado estacionario que relaciona la entrada salida.

Procedemos a encontrar el punto estable de ambos componentes de la misma forma que lo hicimos en el caso 1:

C_{As}=6.0870\dfrac{mol}{litro}

C_{Bs}=1.117\dfrac{mol}{litro}

Con esto obtenemos la siguiente representación en el espacio de estado

\dot{x}=\begin{bmatrix}  -5.7367 &0\\ \\  0.8333 & -4.5411  \end{bmatrix}x+\begin{bmatrix}  3.9130 \\ \\  -1.117  \end{bmatrix}u

y=\begin{bmatrix}  0&1  \end{bmatrix}x

Y podemos obtener la función de transferencia de la misma forma que lo hemos visto en el curso de sistemas dinámicos lineales usando G(s)=C(sI-A)^{-1}B.

G(s)=\dfrac{-1.117s-3.1472}{s^2+10.2778s+26.0508}

En este caso ya se tiene el cero en el semiplano negativo, por lo tanto la respuesta inversa desaparece.

Cero=-\dfrac{3.1472}{1.117}=-2.8175

Van de Vusse caso 2

CASO 3

El punto de operación óptimo para la producción del componente B se logra utilizando la siguiente tasa de disolución:

\dfrac{F_s}{V}=1.2921min^{-1}

C_{As}=4.4949\dfrac{mol}{litro}

C_{Bs}=1.266\dfrac{mol}{litro}

Con esto obtenemos la siguiente representación en el espacio de estado

\dot{x}=\begin{bmatrix}  -3.6337 &0\\ \\  0.8333 & -2.9588  \end{bmatrix}x+\begin{bmatrix}  5.5051 \\ \\  -1.2660  \end{bmatrix}u

y=\begin{bmatrix}  0&1  \end{bmatrix}x

Y podemos obtener la función de transferencia de la misma forma que lo hemos visto en el curso de sistemas dinámicos lineales usando G(s)=C(sI-A)^{-1}B.

G(s)=\dfrac{-1.2660s}{s^2+6.5825s+10.7217}

El cero de esta función de transferencia se ubica en el origen. En este caso la condición de estado estable no cambia en la concentración de B para pequeños cambios en la entrada. Esto es debido a la ganancia cero de la función de transferencia.

 

Van de Vusse caso 3

Bibliografica

  • Bequette, B. (1998). Process Dynamics: Modeling, Analysis, and Simulation. New Jersey: Prentice Hall

Código de Implementación

MATLAB:

SIMULINK:

A continuación se presenta el código de implementación para que lo puedas copiar y pegar en tu computador. Recuerda que para poder ver el código, debes compartir el contenido de esta entrada con alguno de los tres botones mostrados abajo. Asi ayudas a que más personas conozcan la página y se beneficien con esta información.

>>>>> Descargar Archivos de Simulación <<<<<<

Principal:

 

 

Función del Modelo:

 

Ene 172017
 

La decomposición canónica de las ecuaciones de estado nos ayuda a establecer la relación que existe entre la representación por variables de estado y la representación por matriz de función de transferencia.

 

Para comenzar, debemos recordar lo visto en la clase pasada de representación equivalente (Click aca para ver la clase pasada):

 

 

Teníamos que la transformación de similitud me permite representar mi sistema original en otro sistema equivalente, donde este concepto consiste en tomar mi sistema dado por:

\dot{x}(t)=Ax(t)+Bu(t)

y(t)=Cx(t)+Du(t)

Defino una matriz de similitud P, donde:

\xi(t)=Px(t),det(P)\neq0

Multiplico por P a la izquierda del sistema

P\dot{x}(t)=PAx(t)+PBu(t), x(0)

y(t)=Cx(t)+Du(t)

Recordando que: x(t)=P^{-1}\xi(t)

\dot{\xi}(t)=PAP^{-1}\xi(t)+PBu(t)

y(t)=CP^{-1}\xi(t)+Du(t)

Mi sistema transformado vendria dado entonces por:

\dot{\xi}(t)=\bar{A}\xi(t)+\bar{B}u(t)

y(t)=\bar{C}\xi(t)+Du(t)

donde \bar{A}=PAP^{-1}, \bar{B}=PB, \bar{C}=CP^{-1}.

 

Si determinamos la matriz de Controlabilidad para mi sistema original (tal como lo vimos en la entrada de controlabilidad y observabilidad) y determinamos la misma matriz de controlabilidad para mi sistema equivalente, podremos ver que la matriz de transformación de similitud (P), aparece en la relación de las dos matrices:

\bar{C_o}=PC_o

Y lo mismo ocurre en la matriz de Observabilidad:

\bar{O_b}=O_bP^{-1}

Como la matriz de similitud P es inversible (Osea que tiene rango completo) implica en que NOOOO se altera la matriz de Controlabilidad y Observabilidad cuando se hace una transformación de similitud del sistema.

 

Ahora cuando el sistema no sea controlable significa que algunos estados no son alcanzables para la entrada. Osea que no puedo modificar todos los estados generando una determinada entrada, entonces lo que podria hacer en esos casos es utilizar una transformación de similitud para separar el problema en dos partes, una parte donde se consiga controlar el sistema y otra parte donde NO lo consigo controlar a partir de la entrada. Con esa representación mi sistema queda particionado como:

 

\begin{bmatrix}\dot{\xi}\\ \dot{\xi}\end{bmatrix}=\begin{bmatrix}\bar{A}_c & \bar{A}_{12}\\0& \bar{A}_{\bar{c}}\end{bmatrix}  \begin{bmatrix}\xi\\ \xi\end{bmatrix}+\begin{bmatrix}\bar{B}\\ 0\end{bmatrix}u

y=\begin{bmatrix}\bar{C}_c & \bar{C}_{\bar{c}}\end{bmatrix}  \begin{bmatrix}\xi\\ \xi\end{bmatrix}+Du

 

En los estados la parcela controlable es la primera fila y la parcela no controlable es la segunda fila. Siempre que tenga estados no controlables va a aparecer un cero en la matriz de entrada B (como puede verse en la segunda fila de los estados).

 

Para obtener esa representación por medio de la matriz de similitud, lo que hago es saber como determinar la matriz P. Para eso se toma la matriz de controlabilidad de dimension n, dada por:

rank(C_o)=rank(\begin{bmatrix}B &AB &... &A^{n-1}B\end{bmatrix})=n_1<n

Como sabemos que mi matriz no va a ser controlable, pues nuestro rank nos va a dar [/latex]n_1[/latex], quiere decir que esa matriz de controlabilidad solo tendra [/latex]n_1[/latex] columnas linealmente independientes. Entonces se usan [/latex]n_1[/latex] columnas linealmente independientes para formar la inversa de la matriz de similitud P:

P^{-1}=\begin{bmatrix}q_1&...&q_{n_1}&...&q_n\end{bmatrix}

Pero notemos que para formar la matriz de similitud P, es necesario que tenga n columnas (dimension del sistema), entonces lo que hago es adicionar columnas arbitrariamente y que sean linealmente independientes en relación a las [/latex]n_1[/latex] columnas primeras de forma que la matriz P tenga rango completo.

Ahora si apenas se considera el SubSistema Controlable:

\dot{\xi}(t)=\bar{A}_c\xi(t)+\bar{B}_cu(t)

y(t)=\bar{C}_c\xi(t)+Du(t)

Sabemos que la dimension de ese subsistema sera [/latex]n_1[/latex] que era dado por el rango de la matriz de controlabilidad. Entonces aqui aparece un resultado muy interesante, si yo calculo la función de transferencia de ese subsistema menor y la función de transferencia del sistema grande (donde esta la parcela no controlable) vamos a encontrar que van a dar exactamente la misma función de transferencia:

G(s)=C(sI-A)^{-1}B+D=\bar{C}_c(sI-\bar{A}_c)^{-1}\bar{B}_c+D

EJEMPLO:

\dot{x}(t)=\begin{bmatrix}  1 & 1 & 0\\  0 & 1 & 0\\  0 & 1 & 1  \end{bmatrix}x(t)+  \begin{bmatrix}  0&1\\  1&0\\  0&1  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  1 & 1 & 1  \end{bmatrix}x(k)

Determino la matriz de controlabilidad, sabemos que el rango de una matriz nunca sera mayor que el menor de sus dimensiones, en otras palabras sabemos que el rango va a ser menor igual a 2 (el rango minimo de la matriz B es 2) :

rank(C_o)=rank(\begin{bmatrix}B &AB\end{bmatrix})=2<3

rank(C_o)=rank\begin{bmatrix}  0 & 1 & 1 &1\\  1 & 0 & 1 &0\\  0 & 1 & 1 &1  \end{bmatrix}=2<3

Con ese resultado vemos que en ese sistema tengo 2 estados controlables y 1 estado que NO es controlable.

 

Ahora si usamos la idea de tomar dos columnas linealmente independientes de la matriz de controlabilidad mas otra columna cualquiera para formar una matriz de similitud, para hacer la transformación del sistema, voy a conseguir separarlo en una parte controlable y una parte no controlable.

 

De la matriz de controlabilidad voy a escoger las dos primeras las cuales sumadas me dan la tercera columna, asi mi matriz de similitud viene dada por [/latex]P^{-1}=\begin{bmatrix}B &q\end{bmatrix}[/latex]. Aqui la tercera columna es arbitraria, la unica conición que debe cumplir es que sea linealmente independiente con respecto a las otras dos columnas.

P^{-1}=\begin{bmatrix}  0 & 1 & 1 \\  1 & 0 & 0\\  0 & 1 & 0  \end{bmatrix}

Aplico transformación de similitud usando esa matriz y obtengo la siguiente decomposición canonica, \bar{A}=PAP^{-1}, \bar{B}=PB[latex], \bar{C}=CP^{-1}[/latex].

\bar{A}=\begin{bmatrix}  1 & 0 & 0 \\  1 & 1 & 0\\  0 & 0 & 1  \end{bmatrix}, \bar{B}=\begin{bmatrix}  1 & 0 \\  0 & 1 \\  0 & 0  \end{bmatrix}, \bar{C}=\begin{bmatrix}  1 & 2 & 1 \\  \end{bmatrix}

Asi el subsistema controlable viene dado por:

\dot{\xi}(t)=\begin{bmatrix}  1 & 0 \\  1 & 1  \end{bmatrix}\xi(t)+  \begin{bmatrix}  1&0\\  0&1  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  1 & 2  \end{bmatrix}\xi(t)

Ahora si calculo la función de transferencia de ese subsistema y la función de transferencia del sistema original, voy a obtener la misma función de transferencia para ambos.

G(s)=\bar{C}_c(sI-\bar{A}_c)^{-1}\bar{B}_c+D=\begin{bmatrix}\dfrac{s+1}{s^2-2s+1} &\dfrac{2}{s-1}\end{bmatrix}

En matlab puedo calcular esa descomposicion canonica:

>> ctrbf(A,B,C)

Con esto sacamos una primera conclusión. Y es que todo lo que no sea controlable no va a aparecer en la función de transferencia del sistema. Veamos que lo mismo se cumple para la observabilidad.

rank(O_b)=rank(\begin{bmatrix}C \\CA \\\vdots \\CA^{n-1}\end{bmatrix})=n_2<n

Como sabemos que mi matriz no va a ser controlable, pues nuestro rank nos va a dar [/latex]n_2[/latex]. Al igual que el caso anterior la matriz P se calcula con la matriz de observabilidad pero esta vez lo hago con filas linealmente independientes.

P^{-1}=\begin{bmatrix}p_1\\\vdots \\p_{n_2}\\\vdots \\p_n\end{bmatrix}

Nuestro subsistema observable viene dado entonces por:

\begin{bmatrix}\dot{\xi}\\ \dot{\xi}\end{bmatrix}=\begin{bmatrix}\bar{A}_o & 0 \\ \bar{A}_{21}& \bar{A}_{\bar{o}}\end{bmatrix}  \begin{bmatrix}\xi\\ \xi\end{bmatrix}+\begin{bmatrix}\bar{B}_o \\ 0\end{bmatrix}u

y=\begin{bmatrix}\bar{C}_o & \bar{C}_{\bar{o}}\end{bmatrix}  \begin{bmatrix}\xi\\ \xi\end{bmatrix}+Du

Donde también obtendremos que nuestra función de transferencia del subsistema observable será igual a la funcion de transferencia del sistema completo:

G(s)=C(sI-A)^{-1}B+D=\bar{C}_o(sI-\bar{A}_o)^{-1}\bar{B}_o+D

Sistemas No Controlables y No Observables

Para este caso particular, consigo realizar otra transformación de similitud conocida como descomposición de Kalman, donde el sistema es separado en cuatro bloques

\begin{bmatrix}\dot{\xi}_{CO}\\ \dot{\xi}_{C\bar{O}}\\ \dot{\xi}_{\bar{C}O}\\ \dot{\xi}_{\bar{C}\bar{O}}\end{bmatrix}=\begin{bmatrix}\bar{A}_{CO} & 0 &\bar{A}_{13}&0\\ \bar{A}_{21}& \bar{A}_{C\bar{O}}&\bar{A}_{23}&\bar{A}_{24}\\ 0 & 0 &\bar{A}_{\bar{C}O}&0\\ 0& 0&\bar{A}_{43}&\bar{A}_{\bar{C}\bar{O}}\end{bmatrix}  \begin{bmatrix}\xi_{CO}\\ \xi_{C\bar{O}}\\ \xi_{\bar{C}O}\\ \xi_{\bar{C}\bar{O}}\end{bmatrix}+\begin{bmatrix}\bar{B}_{CO} \\ \bar{B}_{\bar{C}O} \\ 0 \\ 0 \end{bmatrix}u

y=\begin{bmatrix}\bar{C}_{CO} &0& \bar{C}_{\bar{C}O}&0\end{bmatrix}  \begin{bmatrix}\xi_{CO}\\ \xi_{C\bar{O}}\\ \xi_{\bar{C}O}\\ \xi_{\bar{C}\bar{O}}\end{bmatrix}+Du

La unica parte controlable y observable son los primeros terminos de cada matriz (\bar{A}_{CO}, \bar{B}_{CO}, \bar{C}_{CO}). Aqui se cumple la misma propiedad, es decir que si se crea un subsitema apenas con la parte controlable y observable y se calcula la función de transferencia, será exactamente igual a la función de transferencia del sistema completo.

 

En matlab existe una función que encuentra ese subsistema

>> sys=ss(A,B,C,D); %Creo el sistema completo

>>[subsys,P] = minreal(sys); %Determina el subsistema y la matriz de transformación

 

EJEMPLO

\dot{x}(t)=\begin{bmatrix}  -1 & 0 & 0 & 0\\  0 & -2 & 0 & 0\\  0 & 0 & -3 & 0 \\  0 & 0 & 0 & -4  \end{bmatrix}x(t)+  \begin{bmatrix}  1\\  0\\  0\\  1\\  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  1&1 & 0&0  \end{bmatrix}x(k)

Código en Matlab, debes compartir para poderlo ver:

 

Condiciones en la forma de Jordan

Cuando tengo un sistema MIMO (Varias entradas y varias salidas), la forma de identificar si mi sistema es Controlable y Observable mirando las matrices A,B y C puede hacerse llevando el sistema a la Forma diagonal, ó como ya lo habíamos visto, en el caso de tener autovalores repetidos en la forma de Jordan.

Para eso supongamos que tenemos el siguiente sistema en la forma diagonal, donde tengo tres bloques de Jordan y mi matriz B y C tienen varias columnas:

\dot{x}(t)=Jx(t)+Bu(t)

y(t)=Cx(t)+Du(t)

Donde J=diag[J_1,J_2] con J_1=diag[J_{11},J_{12},J_{13}] y J_2=diag[J_{21},J_{22}]

 

Lo que implica que la b_{ij} fila de la matriz B está asociada a la ultima fila de la matriz J_{ij} y que la columna c_{ij} de la matriz C está asociada a la primeira coluna da matriz J_{ij}.

 

Viendo en un ejemplo numerico seria como sigue donde tenemos 2 entradas y 2 salidas:

\dot{x}(t)=\begin{bmatrix}  -1 & 1 & 0 & 0\\  0 & -1 & 0 & 0\\  0 & 0 & -2 & 1 \\  0 & 0 & 0 & -2  \end{bmatrix}x(t)+  \begin{bmatrix}  b_{11}&b_{12}\\  b_{21}&b_{22}\\  b_{31}&b_{32}\\  b_{41}&b_{42}\\  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  c_{11}&c_{12} & c_{13}&c_{14}\\  c_{21}&c_{22} & c_{23}&c_{24}  \end{bmatrix}x(k)

1. La ultima fila del bloque de Jordan es esta, por eso la fila que me interesa de la matriz B es esta.

2. La columna de la matriz C asociada a la primera columna del bloque de Jordan. La primera columna de Jordan es esta, entonces la columna en C que me interesa es esta para la primera partición.

3. La segunda partición del bloque de Jordan es la ultima fila de la matriz B y la primera Columna de la matriz C.

 

Son las únicas filas y columnas que me interesan cuando deseo hacer el análisis sea de sistema controlable o observable cuando dicho sistema se encuentra en la representación de la forma diagonal ó de Jordan

 

Si las dos filas escogidas en la matriz B son linealmente independientes significa que el sistema es Controlable y si las dos columnas escogidas en la matriz C son linealmente independientes, quiere decir que mi sistema es Observable.
Miremos el siguiente sistema y tratemos de identificar si el sistema es Controlable y Observable:

\dot{x}(t)=\begin{bmatrix}  \lambda_1 & 1 & 0 & 0 & 0 & 0 & 0\\  0 &\lambda_1 & 0 & 0 & 0 & 0 & 0\\  0 & 0 & \lambda_1 & 0 & 0 & 0 & 0\\  0 & 0 & 0 & \lambda_1 & 0 & 0 & 0\\  0 & 0 & 0 & 0 & \lambda_2 & 1 & 0\\  0 & 0 & 0 & 0 & 0 & \lambda_2 & 1\\  0 & 0 & 0 & 0 & 0 & 0 & \lambda_2  \end{bmatrix}x(t)+  \begin{bmatrix}  0 & 0 & 0\\  1 & 0 & 0\\  0 & 1 & 0\\  1 & 1 & 1\\  1 & 2 & 3\\  0 & 1 & 0\\  1 & 1 & 1  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  1 & 1 & 2 & 0 & 0 & 2 & 1\\  1 & 0 & 1 & 2 & 0 & 1 & 1\\  1 & 0 & 2 & 3 & 0 & 2 & 0  \end{bmatrix}x(k)

Tomando las primeras filas de la Matriz B de cada bloque de Jordan

C=\begin{bmatrix}  1 & 0 & 0\\  0 & 1 & 0\\  1 & 1 & 1\\  1 & 1 & 1\\  \end{bmatrix}

Como solo tengo 3 Filas linealmente independientes, por lo tanto no es rank completo y NO es CONTROLABLE.

Tomando las primeras columnas de la Matriz C de cada bloque de Jordan

O=\begin{bmatrix}  1 & 2 & 0 & 0\\  1 & 1 & 2 & 0\\  1 & 2 & 3 & 0  \end{bmatrix}

Como tengo una columna lleno de ceros sabemos que solo tendremos 3 columnas linealmente independientes, por lo tanto no es rank completo y NO es OBSERVABLE.

 

Sistemas Discretos

Si tenemos el siguiente sistema discreto

\dot{x}[k+1]=Ax[k]+Bu[k] y(t)=Cx[k]+Du[k]

Se dice que el par (A,B) es controlable si para cualquier condicion inicial x[0]=x_0 y cualquier estado final x_1 existe una secuencia de entrada de tamaño finito que lleva el estado desde x_0 para x_1, caso contrario el sistema no es Controlable.

 

Recordemos que el caso discreto es por instanes, osea que puedo tener una entrada en el instante cero, uno, dos, …. u[0],u[1],...u[m].

Entonces si el sistema es de orden 5, voy a necesitar de 5 instantes de tiempo para poder determinar si consigo llevar el estado inicial al estado final y saber si es controlable.

 

Aplicamos el mismo concepto visto en el caso continuo, determinando las matrices de Controlabilidad y Observabilidad.

Ene 122017
 

En esta entrada veremos un segundo ejemplo de un controlador PID implementado en un microcontrolador PIC. Para esta entrada veremos el funcionamiento de este algoritmo aplicado sobre un proceso real. Recuerda que ya habíamos visto un primer ejemplo totalmente detallado de la implementación de un PID, si todavía no lo has visto da CLIC ACA para ir a la entrada del control PID PIC para un horno, pues usaremos casi que el mismo código que fue explicado en esa entrada.

 

El proceso escogido a ser controlado es conocido como un Motor-Generador y es mostrado en la siguiente fotografía:

Motor-Generador

El funcionamiento del proceso es el siguiente: Este proceso consta de dos motores los cuales están acoplados en sus ejes por medio de poleas. Uno de estos motores funciona propiamente como motor, entonces cuando dicho motor es excitado con un voltaje de alimentación comienza a girar a una velocidad proporcional al voltaje de entrada, el giro del motor hace que el eje del segundo motor (Conocido como generador) comience a girar debido al acople entre ejes. Cuando el eje del generador está en movimiento, este comienza a generar voltaje directo el cual puede ser tomado de los dos terminales del motor.

 

Aquí tenemos un sistema a ser Controlado, como variable manipulada tendremos el Voltaje de Entrada del Motor y como variable controlada tendremos el Voltaje generado por el Generador.

 

Adicionalmente puede ser observado en la fotografía, que el sistema cuenta con dos bombillos que simularan la carga del sistema, esto seria la perturbación que constantemente afectan los sistemas de control. Nuestro control debe ser lo suficientemente robusto para rechazar esas perturbaciones y llevar la variable controlada de nuevo para el Set-Point.

 

El circuito básico de este proceso MOTOR-GENERADOR puede ser descargado en el siguiente archivo de Word dando CLIC ACA.

 

Para este ejemplo inicialmente debemos obtener el modelo matemático del proceso (Identificar la planta) y posteriormente con la función de transferencia identificada haremos el controlador PID. Estos fueron los mismos pasos que utilizamos en la entrada anterior del PID del horno.

 

Identificación del Proceso

Mostraré un procedimiento para realizar la identificación del proceso, pero aquí les dejaré su primera tarea. Lo ideal es que lo siguiente que vamos a ver ustedes lo implementen con una comunicación serial (Click aca para ver entrada de comunicación serial), enviando los datos al computador para poder obtener la curva de reacción del proceso y posteriormente obtener la función de transferencia.

 

Como en este momento no poseo un FT232 para poder comunicar el PIC con mi Computador, decidí por hacer una identificación un poco más manual, utilizando un LCD 4×20. El circuito para realizar la identificación es el siguiente:

PID con PIC

NOTA:

En este circuito se aplica lo aprendido en la entrada del Teclado Matricial 4×4 y como se está utilizando un LCD4x20, ya habíamos visto que el compilador CCS C Compiler ya trae lista la librería para trabajar con este teclado, pero solo está lista para trabajarla conectandola al PUERTO B, por eso en el video de la identificación de Youtube, se muestra como hacer la modificación de la libreria para que trabaje con el PUERTO D.

 

 

La idea de una identificación del proceso consiste en lo siguiente:

Tenemos nuestro proceso MOTOR-GENERADOR, el cual no tenemos idea de cual puede ser su representación matemática para poder sintonizar nuestro controlador PID, pero la idea entonces es obtener ese modelo.

Aquí estamos recayendo sobre la teoría de Control Lineal, por eso estamos obteniendo una Función de Transferencia en el dominio de LaPlace.

 

Todos los procesos en la industria y en la vida real en si, son procesos NO LINEALES, recayendo entonces en la teoría no lineal que es mucho mas compleja. Ahora, dado que en la industria en la mayoria de los casos, los procesos siempre están trabajando sombre un punto de operación (Un mismo setpoint) y raramente son movidos para otro punto de operación, podemos aproximar ese proceso no lineal en un modelo Lineal por medio de una función de transferencia.

 

Entonces lo que haremos es exitar el sistema con un set-point o un escalón del mismo valor a nuestro punto de operación, para este ejemplo yo exité el sistema con un escalón del 40%, colocando entonces un PWM en la salida del CCP1 un ancho de pulso proporcional al 40% (Ver entrada del PWM con PIC). Adicionalmente establezco el tiempo de muestreo que quiero tomar datos del MOTOR-GENERADOR (Dado que este sistema es sumamente Rápido, coloque un periodo de muestreo de 10mS). De esa manera el PIC exita el motor por causa de los 40% del PWM y comienza a guardar los valores de voltaje generados por el generador cada 10mS. Estos datos los guarda en un vector de 45 posiciones. Una vez las 45 posiciones son llenadas, el PIC coloca el PWM en Cero y muestra todos los datos en la pantalla del LCD. Este código de identificación se encuentra al Final de este Post o también en el Video de Youtube, donde se explica en detalle.

 

Como la Planta MOTOR-GENERADOR solo acepta en su entrada voltaje de 0 a 5 voltios, y nosotros estamos generando un PWM, es necesario transformar ese PWM en un valor de voltaje continuo que dependa del ancho de pulso. Para eso se implementa el siguiente filtro pasa bajos que hace esa función.

Filtro Pasa Bajos

 

Los datos de voltaje mostrados en el LCD, los copie en un archivo de excel junto con los datos de tiempo (periodo de muestreo) y los datos de entrada o escalón (que son los 40%). Estos datos los muestro tanto en Voltaje como en Porcentaje. (Click aca para bajar el archivo en Excel)

 

La respuesta obtenida en porcentaje se muestra a continuación:

Identificacion del Motor-Generador

Tipicamente, todo proceso industrial es representado por funciones de transferencia de primer o segundo orden. A continuación es expuesto un modelo de primero orden sobre el cual vamos a trabajar para poder modelar el comportamiento de nuestra planta.

P(s)=\dfrac{Ke^{-\theta s}}{\tau s +1}

El modelo está representado en su forma continua, es decir en el dominio de Laplace. donde K es la ganancia del sistema, \tau es la constante de tiempo del proceso y \theta es el retardo del proceso.

Vemos que la gráfica el voltaje generado (linea azul) comienza en 0% y llega hasta 32% y que para conseguir esta respuesta tuvimos que alimentar el motor con una exitación del 40% (linea naranja). Así podemos obtener la ganancia del proceso con la siguiente formula:

K=\dfrac{Y_{final}-Y_{inicial}}{U_{final}-U_{inicial}}

K=\dfrac{32-0}{40-0}=0.8

Notemos que el voltaje generado (linea zul), al comienzo ella se demora en comenzar a subir, ese tiempo que ella se queda en cero, es conocido como el retardo del sistema y en este caso equivale a 0,03 segundos o 30mS

\theta=30

Observemos que el voltaje generado en la gráfica, se comienza a estabilizar mas o menos a los 0.26 segundos, pero debemos restarle los 0,03 segundos del retardo, entonces el tiempo de establecimiento seria 0,23 segundos, así podemos obtener la constante de tiempo de la variable temperatura con la siguiente formula:

T_{establecimieto}=4\tau

Osea que el \tau=57,5.

 

Así que en términos generales, nuestro proceso del MOTOR-GENERADOR está representado por la siguiente función de transferencia:

P(s)=\dfrac{0,8e^{-30 s}}{57,5 s +1}

Con nuestro proceso identificado podemos calcular los parámetros del controlador (k_p,t_i,t_d), para eso utilizamos cualquiera de las tres técnicas vistas en la entrada el PID del horno.

Si te está sirviendo de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Algoritmo del PID

Una ves tenemos nuestro modelo matemático que representa nuestro proceso MOTOR-GENERADOR en el escalón del 40%, podremos calcular nuestro control PID para que controle adecuadamente el proceso en ese punto de operación.

 

 

Para este caso PARTICULAR de mi MOTOR-GENERADOR, es un proceso de muy baja NO LINEALIDAD, es decir que tiene un comportamiento casi lineal, entonces voy a poder controlar con facilidad el proceso en cualquier punto de operación diferente a el 40%. Esto en la mayoría de los procesos NO sucede.

 

Utiliza el control discreto PID que vimos en la entrada pasada  dado por:

C(z^{-1})=\dfrac{u(k)}{e(k)}=\dfrac{q_0+q_1z^{-1}+q_2z^{-2}}{1-z^{-1}}

donde:

q_0=k_p\left [ 1+\dfrac{T}{2t_i}+\dfrac{t_d}{T} \right ] q_1=-k_p\left [ 1-\dfrac{T}{2t_i}+\dfrac{2t_d}{T} \right ] q_2=\dfrac{k_pt_d}{T}

Con esto, la ley de control que vamos a ingresar a nuestro PIC sale del control PID discreto (Despejando u)

u(k)(1-z^{-1})=q_0e(k)+q_1z^{-1}e(k)+q_2z^{-2}e(k)

u(k)-u(k)z^{-1}=q_0e(k)+q_1z^{-1}e(k)+q_2z^{-2}e(k)

u(k)=u(k)z^{-1}+q_0e(k)+q_1z^{-1}e(k)+q_2z^{-2}e(k)

aplicando transformada inversa Z obtenemos la ecuacion en diferencias:

u(k)=u(k-1)+q_0e(k)+q_1e(k-1)+q_2e(k-2)

Como tiempo de muestreo para el control PID utilizamos el criterio:

T=\dfrac{\tau}{20}=\dfrac{57,5}{20}=2,87\approx 3mS

NOTA

Como vemos el tiempo de muestreo debe ser de 3mS, osea que realmente debe ser rápido nuestro microcontrolador para poder alcanzar tal velocidad, para ello por la alta velocidad lo más recomendable es usar un cristal externo de 20Mhz. En mi caso personal, no cuento con este cristal, entonces voy a tener que implementarlo con un cristal de 4Mhz, por lo que no voy a poder cumplir el criterio del tiempo de muestreo. Esto me puede traer serios problemas de estabilidad, pero al hacer las pruebas, mi periodo de muestreo alcanzado fue de 11mS y el sistema consiguio responder adecuadamente, principalmente porque todavía es 5 veces menor que  la constante de tiempo. Si se hubieran presentado problemas de estabilidad por causa que no puedo respetar el periodo de muestreo, lo que hubiera hecho es ponderar un poco el error del controlador PID, es decir multiplicar el error por valores entre (0 – 1).

 

Código de Implementación:

A continuación se presenta el código para que lo copies y lo pegues en tu compilador y puedas reproducirlo. Recuerda que para ver el código debes compartir el contenido de este blog para que más personas se beneficien de esta información.

> Descargar los Archivos de PICC y PROTEUS 8 (Click aca) <<<

Descargar librería del LCD4x20 para puertoD: LCD420D.C

Código de la Identificación

 

Código del Control PID para MOTOR-GENERADOR

 

Oct 062016
 

Para sistemas lineales e invariantes en el tiempo el análisis de estabilidad es una tarea relativamente simple, sin embargo los métodos de estabilidad vistos anteriormente, no pueden ser aplicados en el análisis de sistemas lineales variantes en el tiempo ni tampoco en sistemas no lineales.

 

Una forma genérica para analizar la estabilidad de sistemas es a través de la teoría de estabilidad de Lyapunov, donde dicho análisis es determinado de forma indirecta. Porque no se usa directamente la dinámica del sistema si no que por el contrario se analiza como el sistema se está comportando.

 

Esta idea surgió en el siglo 19 donde todo se basa en energía, en aquella época solo existían sistemas mecánicos, entonces si la energía de ese sistema mecánico va para cero al momento del tiempo tender a infinito, se podria caracterizar un sistema como estable.

 

Se sabe que la energia siempre es una función escalar construida a partir de la dinámica del sistema

 

Las expresiones para definir la energía son infinitas, dependiendo del tipo de aplicación con la cual se desee trabajar (energia potencial cinetica,…)

 

Básicamente se asocia al sistema una función escalar V(||x||) de tal forma que capture la energía asociada a la trayectoria del sistema y voy a analizar si dicha función va para cero en el transcurrir del tiempo.

 

Esta función escalar asocia la norma de los estados del sistema.

 

Si para nuestro analisis comenzamos recordando cual es la energia de una señal en el tiempo:

 

E_x=\sqrt{\int_{0}^{\infty}x(t)^2dt}

E_x^2=\int_{0}^{\infty}x(t)^2dt

Ahí se analiza si esa función va para cero a lo largo del tiempo. Una característica de esa definición de energía es que ella siempre es positiva

Energia Lyapunov

Para la integral poder ir para cero, se analiza que la característica de la gráfica es siempre positiva y que la derivada es siempre negativa porque ella siempre decrece. Y si se cumplen esas dos características, ese sistema siempre va ir para cero.

 

Entonces se verifica que un sistema es internamente estable si la energía es

E_x>0

\dfrac{dE_x}{dt}<0

Ahora si se tiene más de un estado, osea que no es solo un escalar, se puede usar la definición de producto interno

E_x^2=\int_{0}^{\infty}x(t)'x(t)dt

El único problema es que si por ejemplo mi función de energía no tiene una derivada que sea siempre negativa

Energia de Lyapunov

Dicha función es estable porque va para cero, pero no cumple con las dos propiedades, entonces yo no consigo capturar la dinámica de esa energía. El problema de esa abordaje, es que la función que se vaya a utilizar debe ser monotonicamente decreciente. Para conseguir eso se coloca una matriz P en los estados y se intenta hallar la matriz tal que la función sea monotonicamente decreciente

E_x^2=\int_{0}^{\infty}x(t)'Px(t)dt

En conclusión, para determinar la estabilidad por Lyapunov, consiste en poder determinar una función que presente un comportamiento de energía monotono para decir que el sistema es estable, de lo contrario si no se consigue determinar esta función, en términos generales no quiere decir que el sistema sea inestable.

Propiedades de Matrices Simetricas

Para entender como analizar la estabilidad por Lyapunov, es necesario rever algunas propiedades de las matrices simétricas:

 

Si M es una matriz cuadrada (m\times m)

 

Si M=M^T se dice que es una matriz Simetrica.

Si M=-M^T se dice que es una matriz Antisimetrica.

1. Cualquier matriz cuadrada equivale a la suma de una matriz simetrica y una matriz antisimetrica

M=\dfrac{M+M^T}{2}+\dfrac{M-M^T}{2}

2. La función cuadrática x^TMx asociada a una matriz antisimetrica es siempre CERO.

x^TMx=-x^TM^Tx=(-x^TM^Tx)^T=-x^TMx=0

Luego, la función x^TMx con M simétrica o no, es igual a una función cuadrática con matriz simétrica.

Matrices definidas positivas:

Una matriz cuadrada (m\times m) es definida positiva si

 

\forall x \neq 0 \rightarrow x^TMx>0

 

Lo que quiere decir que la Matriz M es definida positiva si la funcion x^TMx es definida positiva.

 

Aplicación a sistemas Lineales

Si tenemos un sistema lineal \dot{x}=Ax y probamos con la función de Lyapunov V(x)=x^TMx>0

 

Primero derivamos la función de Lyapunov

 

\dot{V}(x)=\dot{x}^TMx+x^TM\dot{x}=x^TA^TMx+x^TMAx

\dot{V}(x)=x^T (A^TM+MA)x<0

 

La derivada tiene que ser menor que cero

Como (A^TM+MA) es una forma cuadrática si yo tomo ese resultado que va a ser una matriz que es Simetrica basta con analizar los autovalores de esa matriz, si todos los autovalores son negativos, la matriz es definida negativa lo que indica que la derivada de Lyapunov siempre dará negativa.

 

Por otro lado si en la ecuación cuadratica (x^TMx>0) si M tiene todos los autovalores positivos, esa ecuación siempre va a dar positivo garantizando la estabilidad

 

Ahora esto solo vale si yo conozco la matriz M,si yo no conozco dicha matriz tendria que encontrar una Matriz M simetrica con todos los autovalores positivos de forma que (A^TM+MA) tenga todos los autovalores negativos

 

\exists M=M^T>0 : A^TM+MA<0

 

Para poder resolver ese problema analíticamente, en lugar de trabajar con esa desigualdad, se puede transformar esa expresión en una ecuación lineal.

Si tengo esta ecuación:

 

\exists M=M^T>0, N=N^T>0:A^TM+MA=-N

 

Si esa matriz N tiene todos los autovalores positivos y colocamos un signo (-) alfrente de N, quiere decir que esa matriz es Negativa. De esa forma consigo demostrar que el termino (A^TM+MA) es menor que cero resolviendo una simple ecuación. Esto es conocido como la ECUACIÓN DE LYAPUNOV.

 

De esa forma consigo por medio de un sistema de ecuaciónes encontrar una matriz M. Solo que cuando hallamos esa M, ademas de ser una matriz Simetrica ella debe tener todos los autovalores Positivos.

 

Teorema 1

Si todos los autovalores de A tienen parte real negativa si y solo si para una matriz dada que sea simetrica positiva N a la ecuación de Lyapunov A^TM+MA=-N se consigue una unica solución simetrica y definida positiva M. Indicando que el sistema es Estable caso contrario es Inestable (esto es solo valido para sistemas LTI), para No lineales o variantes en el tiempo, simplemente no se puede concluir.

Teorema 2

Si todos los autovalores de A tienen parte real negativa entonces la ecuación de Lyapunov A^TM+MA=-N tiene una unica solución para una dada matriz N simetrica y definida positiva dada por:

M=\int_{0}^{\infty}e^{A^Tt}Ne^{At}dt

Se define a P=A^TQ+QA como la ecuación de Lyapunov. Y la estabilidad asintotica del sistema está asegurada si siendo P una matriz simetrica definida negativa, se puede encontrar Q como una matriz simetrica definida positiva (Solución de la ecuación de Lyapunov)

Si los autovalores de A cumplen lo siguiente:

– Q es definida positiva si Re(\lambda_i)>0 \forall \lambda_i de A

– Q es definida negativa si Re(\lambda_i)<0 \forall \lambda_i de A

– Q es indefinida si \exists Re(\lambda_i)>0 y \exists Re(\lambda_i)<0 de A

Si te está sirviendo de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




EJEMPLO

 

\dot{x}(t)=\begin{bmatrix}  0 & 1\\  -4 & -5  \end{bmatrix}x(t)

 

Sabemos que el sistema es asintoticamente estable si

 

V(||x||)>0 y [/latex]\dfrac{dV(||x||)}{dt}<0[/latex]

Tomamos como función de Lyapunov la expresión cuadratica:

V(x)=x^TMx>0 y [/latex]\dot{V}(x)=x^T (A^TM+MA)x<0[/latex]

Determinamos los Autovalores de A.

det(\lambda I-A)=0 \rightarrow (\lambda_1=-1, \lambda_2=-4)

Sea N=I y M=\begin{bmatrix}  a & b\\  b & c  \end{bmatrix}

A^TM+MA=-N

\begin{bmatrix}0 & 1\\-4 & -5\end{bmatrix}^T\begin{bmatrix}a & b\\b & c\end{bmatrix}+\begin{bmatrix}a & b\\b & c\end{bmatrix}\begin{bmatrix}0 & 1\\-4 & -5\end{bmatrix}=\begin{bmatrix}-1 & 0\\0 & -1\end{bmatrix}

\begin{bmatrix}0 & -4\\1 & -5\end{bmatrix}\begin{bmatrix}a & b\\b & c\end{bmatrix}+\begin{bmatrix}a & b\\b & c\end{bmatrix}\begin{bmatrix}0 & 1\\-4 & -5\end{bmatrix}=\begin{bmatrix}-1 & 0\\0 & -1\end{bmatrix}

\begin{bmatrix}-4b & -4c\\a-5b & b-5c\end{bmatrix}+\begin{bmatrix}-4b & a-5b\\-4c & b-5c\end{bmatrix}=\begin{bmatrix}-1 & 0\\0 & -1\end{bmatrix}

\begin{bmatrix}-8b & a-5b-4c\\a-5b-4c & 2b-10c\end{bmatrix}=\begin{bmatrix}-1 & 0\\0 & -1\end{bmatrix}

Ecuaciones:

-8b=-1 \rightarrow b=1/8=0.125

2b-10c=-1 \rightarrow c=\dfrac{2(1/8)+1}{10}=0.125

a-5b-4c=0 \rightarrow a=5(1/8)+4(1/8)=9/8=1.125

M=\begin{bmatrix}  1.125 & 0.125\\  0.125 & 0.125  \end{bmatrix}

V(x)=x^TMx>0\rightarrow \begin{bmatrix}x_1^T\\x_2^T\end{bmatrix}\begin{bmatrix}1.125 & 0.125\\0.125 & 0.125\end{bmatrix}\begin{bmatrix}x_1\\x_2\end{bmatrix}=1.125x_1^2+0.250x_1x_2+0.125x_2^2>0

Calculando los autovalores de la matriz M tenemos que:

det(\lambda I-M)=0 \rightarrow (\lambda_1=0.1096, \lambda_2=1.1404)

Como los autovalores son positivos, quiere decir que la Matriz M es simetrica y definida positiva, por lo tanto el sistema es estable.

En MATLAB:

> M=lyap(A,N)

 

Para Sistemas Discretos

Para el caso discreto la función de la energia de una señal viene dada por:

E_x^2=\sum_{k=0}^{\infty}x(k)^TPx(k)

Energia de Lyapunov

Para verificar si la función está decreciendo, siempre la diferencia entre un instante futuro con relación a un instante actual debe dar negativo.

V(x[k])>0 y [/latex]\Delta V(x[k])=V(x[k+1])-V(x[k])<0[/latex]

La ecuación de Lyapunov para el caso discreto viene dado por:

Si se tiene un sistema discreto

x[k+1]=Ax[k]

Y se utiliza la función de Lyapunov mas simple, conocida como función de Lyapunov Cuadrada:

V(x[k])=x[k]^TMx[k]>0

Y ya sabemos que como es una ecuación cuadrática, la matriz M tiene que tener todos los autovalores positivos y también ser simétrica.

Se calcula entonces la variación de la función:

\Delta V(x[k])=V(x[k+1])-V(x[k])<0

\Delta V(x[k])=x[k+1]^TMx[k+1]-x[k]^TMx[k]<0

\Delta V(x[k])=(Ax[k])^TM(Ax[k])-x[k]^TMx[k]<0

\Delta V(x[k])=x[k]^TA^TMAx[k]-x[k]^TMx[k]<0

\Delta V(x[k])=x[k]^T(A^TMA-M)x[k]<0

Y esta también es una forma cuadratica.

Para determinar si es estable basta con analizar si todos los autovalores de la matriz (A^TMA-M<0) son menosres que cero.

Para resolver analíticamente, se iguala la matriz a otra matriz negativa N, igual que en el caso continuo

A^TMA-M=-N

Si consigo determinar una matriz M definida positiva puedo decir que mi sistema dinámico Lineal e invariante en el tiempo es estable, recordando que si no es definida positiva seria inestable.

Para el caso de Lineal y Variante en el tiempo y No lineal en el caso de la matriz M ser definida negativa, no se puede determinar nada, y se tendria que seguir buscando otra función de Lyapunov.

 

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Sep 262016
 

La condición de estabilidad interna esta asociada a como los estados se comportan a partir de la imposición de una condición inicial. Lo que quiere decir, para este análisis de estabilidad interna, se desconsidera la entrada del sistema y únicamente se analizan los estados cuando se parte de una condición inicial.

 

Recordando que el sistema que hemos venido trabajando se representa por:

\dot{x}=Ax+Bu, y=Cx+Du

y que una matriz de función de transferencia es calculada por:

G(s)=C(sI-A)^{-1}B+D=\dfrac{1}{det(sI-A)}CAdj(sI-A)B+D

 

Lo que demuestra que idealmente POLOS y AUTOVALORES son para ser idénticos idealmente, pues se evidencia que en la función de transferencia los polos están relacionados al determinante de la matriz A, lo que me permite obtener la ecuación característica de la matriz.

 

Solo que los autovalores y los polos NO son identicos cuando existen cancelamientos por la izquierda por la mariz C y por la derecha por la matriz B:

CAdj(sI-A)B y det(sI-A)

Esos modos cancelados no apareceran en G(s) y por lo tanto no van a influenciar en la estabilidad Entrada – Salida.

Consideremos el siguiente ejemplo

\dot{x}(t)=\begin{bmatrix}  1 & 0\\  0 & -1  \end{bmatrix}x(t)+\begin{bmatrix}  1 \\  1  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  0 & 1  \end{bmatrix}x(t)

Donde la función de transferencia es:

G(s)=\dfrac{1}{s+1}-> Sistema BIBO estable.

Observe que la función de transferencia es de orden uno y que el sistema tiene dos estados, entonces ciertamente aquí hubo un cancelamiento. Donde uno de los autovalores es -1 que es el que aparece en la función de transferencia, solo que hay otro autvalor que NO aparece en la representación entrada-salida.

Si fueramos a analizar la respesta de los estados con un escalón unitario tendriamos que:

x(t)=\begin{bmatrix}  e^{t}-1 \\  1-e^{-t}  \end{bmatrix}u(t)

el e^{t} del primer estado va para infinito por causa de su exponencial ser positivo. Solo que ese estado x_1(t) no aparece en la representación entrada-salida del sistema.
Entonces si fuéramos a analizar únicamente la función de transferencia del sistema, dicho sistema es ESTABLE, pero, se puede ver que este mismo sistema es internamente inestable.
La estabilidad interna intenta capturar la evolución de todos los estados a lo largo del tiempo.
Teorema:
Para tener estabilidad interna se debe cumplir que toda condición inicial la trayectoria de los estados x(t) tienden a cero cuando el tiempo tiende a infinito (Estabilidad asintomática).

Para una dada condición inicial x(0)=x_0, todos los estados deben asintoticamente converger para cero a medida que el tiempo tiende a infinito.

Como solo me tengo que preocupar con la respuesta de la condición inicial del sistema, es mucho más facil de calcular:

x(t)=e^{\mathbf{A}t}x_0

Donde el exponencial debe tender a cero cuando el tiempo va para infinito, para que sea asintoticamente estable. Y esto esta relacionado directamente con la respuesta al impulso del sistema.

Para analizar la estabilidad interna del sistema analizamos apenas la matriz A.

Suponiendo que la matriz A esta en la forma diagonal la respuesta en el tiempo para una condición inicial x(0) es dada por:

x(0)=\begin{bmatrix}  x_1(0) \\\vdots \\  x_n(0)  \end{bmatrix}\rightarrow x(t)=\begin{bmatrix}  x_1(0)e^{\lambda_1t} \\\vdots \\  x_n(0)e^{\lambda_nt}  \end{bmatrix}

Todos los estados del sistema van a tender para cero cuando el tiempo tiende para infinito si y solamente si la parte real de los autovalores es negativa.

EJEMPLO.

Si retomamos el caso del ejemplo de la entrada de Estabilidad Entrada – Salida, teníamos que:

\dot{x}(t)=\begin{bmatrix}  -2 & 0\\  0 & 5  \end{bmatrix}x(t)+\begin{bmatrix}  1 \\  0  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  -0.1 & 4  \end{bmatrix}x(t)

Donde la función de transferencia es:

G(s)=C(sI-A)^{-1}B=\dfrac{-0.1}{s+2}\rightarrow Sistema BIBO Estable.

Como se evidencia en la función de transferencia, el sistema es BIBO Estable visto del punto de vista Entrada-Salida, pero también notamos que el sistema es de dos estados y que la función de transferencia es de orden uno, lo que nos indica que hubo un cancelamiento de uno de los autovalores.

Analicemos entonces la estabilidad interna del sistema, si consideremos la siguiente condición inicial x(0)=[1,1]'

De esa forma la respuesta en el tiempo para la condición inicial x(0) es dada por

x(t)=\begin{bmatrix}  1e^{-2t} \\  1e^{5t}  \end{bmatrix}

Aqui se puede apreciar que el sistema es INTERNAMENTE INESTABLE, porque el estado x_2 se vá para el infinito cuando el tiempo es infinito.

Si te está sirviendo de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




SISTEMAS DISCRETOS

Para sistemas discretos tenemos que un sistema representado por:

x[k+1]=Ax[k]+Bu[k]

Donde para analizar la estabilidad interna del sistema, únicamente se considera la respuesta al estado cero

x[k+1]=Ax[k]

Suponiendo que la matriz A esta en la forma diagonal la respuesta en el tiempo para una condición inicial x(0) es dada por

x(0)=\begin{bmatrix}  x_1(0) \\\vdots \\  x_n(0)  \end{bmatrix}\rightarrow x(k)=\begin{bmatrix}  x_1(0)\lambda_1^k \\\vdots \\  x_n(0)\lambda_n^k  \end{bmatrix}

Para el caso discreto todos los módulos de los autovalores tienen que ser menor que uno.

EJEMPLO

Determinar la estabilidad asintotica del siguiente sistema

x[k+1]=\begin{bmatrix}  0 & 1\\  -0.45 & 1.40  \end{bmatrix}x[k]

Como voy analizar la estabilidad interna, no me interesa la salida, únicamente me interesa la matriz dinámica A.

La ecuación característica:

det(\lambda I-A)=0

\lambda^2-1.4\lambda+0.45=0\rightarrow (\lambda-0.9)(\lambda-0.5)=0

Entonces es un sistema es internamente Estable, porque los autovalores son de modulo menor que uno. Lo que quiere decir que si yo aplico cualquier condición inicial, al momento en que k tienda a infinito todo va para cero.

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Sep 252016
 

Los sistemas dinámicos son proyectados para realizar determinado tipo de tareas o para procesar una determinada señal.

 

Pero debemos entender que si un sistema dinámico no es estable, este no podrá desempeñar su función cuando una determinada señal de entrada es aplicada al sistema.

Por lo tanto en la práctica un sistema dinámico debe operar siempre en una región de estabilidad.

 

Esta definición de estabilidad se realiza con base en la representación matemática usada:

  • Estabilidad Entrada-Salida
  • Estabilidad Interna.

Ya habíamos visto que la respuesta en el tiempo de los estados de un sistema LTI podia dividirse en dos partes:

x(t)=x_{zi}(t)+x_{zs}(t)

Donde x_{zi}(t) es conocida como respuesta a la entrada cero y x_{zs}(t) como respuesta al estado cero.

Y que para un sistema con la forma de \dot{x}=Ax+Bu, y=Cx+Du la respuesta de los estados:

x_{zi}(t)=e^{\mathbf{A}t}x(0) y x_{zs}(t)=\int_{0}^{t}e^{\mathbf{A}(t-\tau)}Bu(\tau)d\tau

La estabilidad del punto de vista entrada salida considera la respuesta al estado zero, ejemplo x(0)=0

y(t)=C\int_{0}^{t}e^{\mathbf{A}(t-\tau)}Bu(\tau)d\tau+Du(t)

Estabilidad Entrada Salida

Una entrada u(t) limitada implica en que la salida y(t) es limitada (bounded-input bounded-output BIBO estable)

El analisis de estabilidad BIBO puede ser analizado utilizando la respuesta al impulso

y(t)=\int_{0}^{t}g(t-\tau)u(\tau)d\tau

donde g(t) es la respuesta al impulso del sistema.

Teorema

Un sistema SISO es BIBO estable si y solamente si g(t) es absolutamente integrable en el intervalo de [0,\infty)

\int_{0}^{\infty}|g(t)|dt\leq M<\infty

Supongamos que la siguiente figura representa la respuesta del impulso del sistema:

Respuesta al impulso

Basta integrar desde cero a t del modulo del sistema

Respuesta al impulso

Al integrar, se puede apreciar que el sistema evidentemente es un sistema estable, pues el area de la curva es un area finita y limitada. Esto nos dice que de forma general, un sistema con una entrada impulso debe de tender para cero en el tiempo.

 

Si el sistema no va para cero al momento de hacer la integral va a dar ilimitada.

 

EJEMPLO 1

Considere el siguiente sistema Lineal

\dot{x}(t)=\begin{bmatrix}  -2 & 0\\  0 & -5  \end{bmatrix}x(t)+\begin{bmatrix}  1 \\  3  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  -0.1 & 4  \end{bmatrix}x(t)

y(t)=\begin{bmatrix}  -0.1 & 4  \end{bmatrix}  \begin{bmatrix}  \int_{0}^{t}e^{-2(t-\tau)}u(\tau)d\tau\\  \int_{0}^{t}e^{-5(t-\tau)}3u(\tau)d\tau  \end{bmatrix}

Para u(t)=1 se obtiene que y(t)=-0.05e^{-2t}(e^{2t}-1)+2.4e^{-5t}(e^{5t}-1)

Sistema estable, pues con una entrada limitada se tiene una salida limitada.

 

EJEMPLO 2

Considere el siguiente sistema Lineal

\dot{x}(t)=\begin{bmatrix}  -2 & 0\\  0 & 5  \end{bmatrix}x(t)+\begin{bmatrix}  1 \\  0  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  -0.1 & 4  \end{bmatrix}x(t)

y(t)=\begin{bmatrix}  -0.1 & 4  \end{bmatrix}  \begin{bmatrix}  \int_{0}^{t}e^{-2(t-\tau)}u(\tau)d\tau\\  \int_{0}^{t}e^{5(t-\tau)}0u(\tau)d\tau  \end{bmatrix}

Para u(t)=1 se obtiene que y(t)=-0.05e^{-2t}(e^{2t}-1). Este sistema es estable a pesar de que tenga un autovalor positivo. El estado de X2 va para infinito, solo que ese estado no aparece en la salida. En este caso se va a apreciar la DIFERENCIA que existe entre el concepto de estabilidad ENTRADA-SALIDA y el concepto de Estabilidad Interna.

Una de las ventajas de representar un sistema por variables de estado es que podemos tener todos los modos dinámicos del sistema inclusive aquellos que no tienen influencia en la salida.

 

Utilizando la transformada de Laplace podemos llegar a un teste de estabilidad BIBO facil

Si g(s)=L(g(t)) entonces

 

g(s)=\sum_{i=1}^{m}\frac{1}{(s-p_i)} +\sum_{i=n-m}^{n}\frac{b_i(s+\alpha_i)+c_i\beta_i}{(s+\alpha_i)^2+\beta_i^2}

Aplicando L^{-1}(g(s)u(s))

 

y(t)\approx\sum_{i=1}^{m}a_ie^{p_it}+\sum_{i=n-m}^{n}e^{\alpha_it}(c_isen(\beta_it)+b_icos(\beta_it))

Por lo tanto el sistema es BIBO estable si y solamente si todos los polos de g(s) tienen parte real negativa.

 

Para sistemas MIMO

  • Un sistema MIMO con respuesta al impulso G(t)=[g_{ij}(t)] es BIBO estable si y solamente si cada elemento [g_{ij}(t) es absolutamente integrable en el intervalo de [0,\infty)
  •  Un sistema MIMO con una matriz de función de transferencia propia G(s)=[g_{ij}(s)] es BIBO estable si y solamente si los polos de cada elemento g_{ij}(s) tienen parte real negativa.

Para sistemas Discretos

Para sistemas discretos la estabilidad entrada salida es caracterizada por la respuesta al impulso discreta.

y[k]=\sum_{m=0}^{k}g[k-m]u[m]

Una secuencia de entrada se dice que es limitada si

u[k]:|u[k]|\leq u_m<\infty,k=0,1,2,...

Un sistema discreto LTI es BIBO estable si toda secuencia de entrada limitada genera una secuancia de salida limitada (considerando condiciones iniciales nulas)

Un sistema discreto LTI con respuesta al impulso g[k] es BIBO estable si y solamente si g[k] es absolutamente sumable en el intervalo de [0,\infty).

\sum_{k=0}^{\infty}|g[k]|\leq M<\infty

utilizando el concepto de transformada Z podemos facilmente verificar la estabilidad BIBO de un sistema discreto.

Recordemos que:

y(z)=Z(g[k]u[k])\approx \sum_{j=1}^{n}\dfrac{1}{z-p_i}

aplicando Z^{-1}(y(z))

y[k]\approx \sum_{j=1}^{n}p_i^k

Por lo tanto esa secuencia es sumable si todos los polos de g(z) tienen un modulo menor que 1.

Para sistemas Discretos MIMO

  1. un sistema multivariable discreto con respuesta al impulso G[k]=[g_{ij}[k]] es BIBO estable si y solamente si cada elemento [g_{ij}[k] es absolutamente sumable en el intervalo de [0,\infty)
  2. Un sistema MIMO con una matriz de función de transferencia propia G(z)=[g_{ij}(z)] es BIBO estable si y solamente si los polos de cada elemento g_{ij}(z) tienen modulo menor que uno.

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Sep 192016
 

Un sistema es representado usando variables de estado a través de realizaciones, las cuales consisten de una ecuación de estado y una ecuación de salida (en términos vectoriales).

 

Las realizaciones básicamente son descripciones de un mismo sistema que dependen del método usado para obtenerlas y que en ciertos casos permiten visualizar o calcular facilmente algunas de las características o propiedades del sistema.

 

Cuando usamos transformaciones de similitud para representar nuestro sistema lineal de forma diferente, se puede comprobar que a pesar de cambiar la representación del sistema, la respuesta del sistema no cambia.

 

Empecemos comprobando esto con el siguiente ejemplo:

\dot{x}(t)=\begin{bmatrix}  0 & -1\\  1 & -1  \end{bmatrix}x(t)+\begin{bmatrix}  1 \\  0  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  0 & 1  \end{bmatrix}x(t)

Donde la función de transferencia viene dada por:

G(s)=C(sI-A)^{-1}B=\dfrac{1}{s^2+s+1}

Donde podemos redefinir los estados a traves de la siguiente matriz de Similitud.

\xi(t)=Px(t)

donde P es una matriz cualquiera de transformación de estados que debe ser inversible

P=\begin{bmatrix}1 & 0\\1 & -1\end{bmatrix}

La caracteristica de que esa matriz sea inversible, es para poder devolverme al estado original, osea volver al estado x(t):

x(t)=P^{-1}\xi(t)

La nueva representación del sistema en terminos del nuevo estado viene dado por:

\dot{\xi}(t)=\begin{bmatrix}1 & 0\\1 & -1\end{bmatrix}\xi(t)+\begin{bmatrix}  1 \\  1  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  1 & -1  \end{bmatrix}\xi(t)

La función de transferencia en esta nueva representación es:

G(s)=C_{\xi}(sI-A_{\xi})^{-1}B_{\xi}=\dfrac{1}{s^2+s+1}

Note que la funcion de transferencia (representación entrada-salida) no se altera cuando se usa una transformación de similitud. Lo unico que cambia es el valor de los estados. Obviamente para que esto se cumpla, la condición inicial de la representación inicial es diferente a la condición inicial del sistema transformado.
En resumen, la transformación de similitud consiste en tomar mi sistema dado por:

\dot{x}(t)=Ax(t)+Bu(t)

y(t)=Cx(t)+Du(t)

Defino una matriz:

\xi(t)=Px(t),det(P)\neq0

Multiplico por P a la izquierda del sistema

P\dot{x}(t)=PAx(t)+PBu(t), x(0)

y(t)=Cx(t)+Du(t)

Recordando que: x(t)=P^{-1}\xi(t)

\dot{\xi}(t)=PAP^{-1}\xi(t)+PBu(t)

y(t)=CP^{-1}\xi(t)+Du(t)

Esto en MATLAB puede ser calculado como:

>> [Ab,Bb,Cb,Db]= ss2ss(A,B,C,D,P)

Es posible también utilizar otros tipos de representaciones del sistema, como por ejemplo usar las formas canonicas del sistema que son formas especiales de nuevas representaciones.

 

Las realizaciones más usados son las formas canónicas: controlable, observable, de Jordan o diagonal o en paralelo, cascada o serie y las realizaciones físicas en las cuales los estados están asociados a la energía almacenada dinámicamente en el sistema.

 

Esas formas tienen una estructura conocida, que para determinado tipo de problemas es mucho más interesante usar esa representación, que la estructura original del sistema.

En MATLAB puedo calcularlo como:

>> [Ac,Bc,Cc,Dc]= canon(A,B,C,D,’tipo’)

La variable ‘tipo’ define la forma canonica que se quiere calcular.

Por ejemplo, considere un sistema con una matriz diagonal con autovalores reales \lambda_1,\lambda_2 y complejos [/latex]\alpha\pm j\beta[/latex]

J=\begin{bmatrix}\lambda_1 & 0& 0& 0\\0& \lambda_2 & 0& 0\\0& 0 & \alpha+ j\beta& 0\\0& 0 & 0& \alpha- j\beta\end{bmatrix}

Esa representación no es muy util en la practica, pues la matriz no es real.

Para esto existe una versión para cuando los autovalores son complejos conjugados conocida como representación modal, donde en la diagonal principal quedan la parte real de los autovalores y en la diagonal secundaria queda la parte imaginaria de los autovalores.

\bar{A}=\begin{bmatrix}\lambda_1 & 0& 0& 0\\0& \lambda_2 & 0& 0\\0& 0 & \alpha& \beta\\0& 0 & -\beta& \alpha\end{bmatrix}

En Matlab:
>> [At,Bt,Ct,Dt]= canon(A,B,C,D,’modal’)

La matriz de transformacion de similitud para la forma modal es:

P^{-1}=\begin{bmatrix}v_1 & v_2& Re(v_3)& Im(v_3)\end{bmatrix}

Esta matriz de transformación se basa en los autovectores de la matriz, donde el primer y segundo autovector son referentes al primer y segundo autovalor real, luego se toma la parte real del tercer autovalor y la parte imaginaria del tercer autovalor y si hago el producto [/latex]PAP^{-1}[/latex] llego a la matriz [/latex]\bar{A}[/latex].

Ejemplo

Obtener la forma modal del siguiente sistema:

\dot{x}(t)=\begin{bmatrix}  0 & 1\\  -5 & -2  \end{bmatrix}x(t)+\begin{bmatrix}  0 \\  1  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  1 & 0  \end{bmatrix}x(t)

ECUACIÓN CARACTERISTICA:

(\lambda I-A)=\lambda^2+2\lambda+5=0

\lambda_{1,2}=-1\pm j2

forma modal

\dot{\bar{x}}(t)=\begin{bmatrix}  -1 & 2\\  -2 & -1  \end{bmatrix}\bar{x}(t)+\bar{B}u(t)

y(t)=\bar{C}\bar{x}(t)+\bar{D}u

La matriz de transformacion de similitud para la forma modal es:

P^{-1}=\begin{bmatrix}Re(v_1)& Im(v_1)\end{bmatrix}

 

Si te está sirviendo de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




 

Realización de Funciones de Transferencia

Ya habiamos visto que podiamos pasar facilmente de espacio de estados para funcion de transferencia

\dot{x}(t)=Ax(t)+Bu(t)

y(t)=Cx(t)+Du(t)

usando

G(s)=C(sI-A)^{-1}B+D

El problema es la realización inversa (pasar de una función de transferencia para el espacio de estados).

Para sistemas SISO se puede obtener una realización de G(s) utilizando las llamadas formas compañeras.

G(s)=\dfrac{b_1s^{n-1}+...+b_{n-1}s^{n-1}s+b_n}{s^n+a_1s^{n-1}+...+a_{n-1}s+a_n}+d

Forma Compañera Observable es:

\dot{x}(t)=\begin{bmatrix}  0 & 0 & \cdots & 0 & -a_n\\  1 & 0 & \cdots & 0 & -a_{n-1}\\  0 & 1 & 0 & \cdots & -a_{n-2}\\  \vdots & \vdots &\ddots & \vdots & \vdots\\  0& \cdots & 0& 1 & -a_1  \end{bmatrix}x(t)+  \begin{bmatrix}  b_n\\  b_{n-1}\\  b_{n-2}\\  \vdots \\  b_1  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  0 & \cdots & 0 & 0 & 1  \end{bmatrix}x(t)+du(t)

Nota: No toda función de transferencia es realizable. Por ejemplo sistemas com Retardo.

Para sistemas MIMO tenemos que:

G(s)=\dfrac{1}{det(sI-A)}N(s)+D

donde N es una matriz cuyos elementos son polinomios.

N(s)=CAdj(sI-A)B

Descomponiendo N(s)

N(s)=N_1s^{n-1}+N_2s^{n-2}+...+N_{n-1}s+N_n

N_i\in R^{q\times p}

Una particularización del sistema MIMO

\dot{x}(t)=\begin{bmatrix}  -a_1I_p & -a_2I_p & \cdots & -a_{n-1}I_p & -a_nI_p\\  I_p & 0 & \cdots & 0 & 0\\  0 & I_p & 0 & \cdots & 0\\  \vdots & \vdots &\ddots & \vdots & \vdots\\  0& \cdots & 0& I_p & 0  \end{bmatrix}x(t)+  \begin{bmatrix}  I_p\\  0\\  0\\  \vdots \\  0  \end{bmatrix}u(t)

y(t)=\begin{bmatrix}  N_1 & N_2 & \cdots & \cdots & N_n  \end{bmatrix}x(t)+Du(t)

EJEMPLO

Considere el siguiente sistema MIMO.

G(s)=\begin{bmatrix}  \dfrac{1}{s+1} & -\dfrac{1}{s+5}\\  \dfrac{1}{s+1} & \dfrac{5}{2}\dfrac{s+2}{s+5}  \end{bmatrix}

Obtener una realización de G.

En la matriz de funciones de transferencia G, para poder calcular el determinate, es necesario separar la parte constante debido a que en la expresión del sistema MIMO, los polinomios N(s) debe ser de grado n-1.

Lo primero que se hace, es obtener la matriz D, y substraerla de la representacion matricial. Como puede observarse el unico elemento que tiene constante, o respuesta instantanea es el elemento G_{22}.

s+2\ \ \ \ \underline{|s+5}\\  >\ \ \ \ -3\ \ \ \ \ \ 1

\dfrac{5}{2}\left(1+\dfrac{-3}{s+5}\right)

G(s)=\begin{bmatrix}  \dfrac{1}{s+1} & -\dfrac{1}{s+5}\\  \dfrac{1}{s+1} & \dfrac{-7.5}{s+5}  \end{bmatrix}+\begin{bmatrix}  0 & 0\\  0 & \dfrac{5}{2}  \end{bmatrix}

Esta ultima matriz corresponde a la matriz D.

El siguiente paso es calcular el determinante de G

det(G(s))=\dfrac{-7.5}{(s+1)(s+5)}+\dfrac{1}{(s+1)(s+5)}

det(G(s))=\dfrac{-6.5}{(s+1)(s+5)}

Se puede ver que los polos de G son (s+1)(s+5) asi reescribo la matriz como

\tilde{G}(s)=\dfrac{1}{(s+1)(s+5)}\begin{bmatrix}  s+5 & -(s+1)\\  s+5 & -7.5(s+1)  \end{bmatrix}

\tilde{G}(s)=\dfrac{1}{s^2+6s+5}\begin{bmatrix}  s+5 & -(s+1)\\  s+5 & -7.5(s+1)  \end{bmatrix}

Donde N(s) es:

N(s)=\begin{bmatrix}  s+5 & -(s+1)\\  s+5 & -7.5(s+1)  \end{bmatrix}

La matriz N(s) se puede expresar como:

N(s)=N_1s+N_2

N(s)=\begin{bmatrix}  1 & -1\\  1 & -7.5  \end{bmatrix}s+\begin{bmatrix}  5 & -1\\  5 & -7.5  \end{bmatrix}
Expresando el sistema en la forma compañera

\dot{x}=\begin{bmatrix}  -6I_2 & -5I_2\\  I_2 & 0_2  \end{bmatrix}x+\begin{bmatrix}I_2\\0_2\end{bmatrix}u

y=\begin{bmatrix}N_1&N_2\end{bmatrix}x+\begin{bmatrix}0 & 0\\0 & 1\end{bmatrix}u

Expresado con todos los terminos expandidos seria representado como:

\dot{x}=\begin{bmatrix}  -6 &0& -5&0\\  0& -6& 0&-5\\  1&0 & 0&0\\  0&1&0&0  \end{bmatrix}x+\begin{bmatrix}1&0\\0&1\\0&0\\0&0\end{bmatrix}u

y=\begin{bmatrix}1 & -1&5 & -1\\1 & -7.5&5 & -7.5\end{bmatrix}x+\begin{bmatrix}0 & 0\\0 & 1\end{bmatrix}u

 

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Sep 182016
 

Una forma de obtener la representación discreta del sistema representado en variables de estados es usando la aproximación de Euler

\dot{x}(t)\approx \dfrac{x(t+T)-x(t)}{T}

Con esto podemos aproximar un sistema continuo tomando el próximo valor del estado menos el valor actual dividido por el periodo de muestreo. Note que esta ecuación es aproximada a una derivada. Esta aproximación se conoce como Euler para al frente.

Usando esa aproximación se tiene que:

x(t+T)=x(t)+TAx(t)+TBu(t)

y así estimando apenas los valores de x(t) y de y(t) apenas en los instantes t=KT se consigue llegar al siguiente modelo aproximado.

x[k+1]=(I+TA)x[k]+TBu[k]

y[k]=Cx[k]+Du[k]

Discretización por ZOH

ZOH

si la entrada u(t) es generada por un computador digital seguido por un conversor DA y manteniendo el valor durante el periodo de muestreo se tiene que:

u(t)=u(kT)=u[k], para kT< t \leq (k+1)T

para esa señal de entrada la solución del sistema continuo es dado por:

x[k]=e^{AkT}x(0)+\int_{0}^{kT}e^{A(kT-\tau)}Bu(\tau)d\tau

x[k+1]=e^{A(k+1)T}x(0)+\int_{0}^{(k+1)T}e^{A((k+1)T-\tau)}Bu(\tau)d\tau

Manipulando la expresión anterior

x[k+1]=e^{AT}\left [e^{AkT}x(0)+\int_{0}^{kT}e^{A(kT-\tau)}Bu(\tau)d\tau \right]+\int_{kT}^{(k+1)T}e^{A(kT+T-\tau)}Bu(\tau)d\tau

x[k+1]=e^{AT}x[k]+\left [\int_{0}^{T}e^{A\alpha}d\alpha\right]Bu[k]

Si las salidas también son determinadas en los instantes t=KT.

x[k+1]=A_dx[k]+B_du[k]

y[k]=C_dx[k]+D_du[k]

donde las matrices del modelo discretizado son:

A_d=e^{AT},B_d=\left (\int_{0}^{T}e^{A\alpha}d\alpha\right)B,C_d=C,D=D_d

La determinación de la matriz B puede simplificarse cuando A sea no singular.

B_d=A^{-1}(A_d-I)B

Este modelo no tiene ninguna aproximación y es una solución exacta del sistema en los instantes de tiempo t=KT

En MATLAB puede calcularse como:

>>[Ad,Bd]=c2d(A,B,T);

La expresión de la solución de sistemas discretos se puede obtener fácilmente como:

x[1]=Ax[0]+Bu[0]

x[2]=Ax[1]+Bu[1]=A^2x[0]+ABu[0]+Bu[1]

generalizando el termino k

x[k]=A^kx[0]+\sum_{m=0}^{k-1} A^{k-1-m}Bu[m]

y[k]=CA^kx[0]+\sum_{m=0}^{k-1} CA^{k-1-m}Bu[m]+Du[k]

 

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.




Sep 182016
 

En esta entrada aprenderemos como obtener la respuesta temporal de un sistema lineal utilizando la matriz de transicion de estados:

 

1. Utilizando transformada de laplace.

 

2. Utilizando el teorema de Cayley-Hamilton

 

3. Utilizando la representación en la forma diagonal y de Jordan.

 

 

Ya habiamos visto que un sistema lineal puede estar representado por la siguiente ecuación integral, (en la forma de representación Entrada-Salida)

 

y(t)=\int_{t0}^{t}g(t,\tau )u(\tau )d\tau

 

O tambien puede estar dado por la descripción en variables de estado:

 

\dot{x}(t)=Ax(t)+Bu(t)

y(t)=Cx(t)+Du(t)

 

Entre tanto, no existe una forma simple de obtener una solución simple analítica para determinar la ecuación de la integral. Numericamente seria facil de resolver si yo discretizo esa integral y la convierto en una sumatoria y podria resolver tranquilamente numericamente.

y(k\Delta)=\sum_{m=k_0}^{k}g(k\Delta,m\Delta)+u(m\Delta)\Delta

Para sistemas lineales e invariantes en el tiempo, yo puedo representar mi sistema lineal a través de la transformada de Laplace para de esa forma conseguir hallar una solución analítica.

y(s)=g(s)u(s)

Donde g(s) es una función de transferencia en el dominio complejo de Laplace.

Si en este caso, nuestra señal de entrada u(s), tambien es una función racional en s, la transformada de Laplace de la señal de salida (y(s)) puede ser determinada de la siguiente forma:

– Calculando los polos de la función de transferencia de g(s).

– Aplicando algún método de fracciones parciales.

– A través de la antitransformada de Laplace, para obtener y(t)

y(t)=L^{-1}(g(s) \times u(s)) = L^{-1}\left \{ \sum_{i}\frac{a_i}{(s+p_i)^{n_i}} +\sum_{j}\frac{b_is+c_j}{(s^2+d_js+e_j)}\right \}

No obstante, en el caso de polos repetidos y de polos complejos conjugados la determinación numerica de la transformada inversa de Laplace es generalmente sensible a errores de aproximación.

Así que desde el punto de vista numerico, la mejor solución para poder obtener la respuesta de un sistema dinámico es a través de la representación por variables de estado.

\dot{x}(t)=Ax(t)+Bu(t)

y(t)=Cx(t)+Du(t)

Ya hemos venido trabajando con esa representación, e incluso ya hemos obtenido la representación de estados de sistemas a través de las ecuaciones diferenciales, ahora la pregunta es:

¿Como puedo yo determinar la salida de un sistema, y(t), a partir de una condición inicial x_0 y una entrada u(t) conocidas?

Analicemos como podemos llegar a la representación de un sistema lineal teniendo condiciones NO nulas:

Esto se conoce como la resolución de la ecuación completa, para eso partimos de la representación del estado

\dot{x}(t)=Ax(t)+Bu(t)

Para poder calcular la solución de esa ecuación se utiliza el método de variación de las constantes, según el cual se va a suponer que existe una z(t) que permita que la solución completa pueda expresarse de la forma: (Ver: http://www.giematic.unican.es/edosso/ejercicios/Teorvarcoef.html)

x(t)=\Phi (t,t_0)z(t)

Esta hipótesis se puede verificar, si se encuentra una z(t) que cumpla esa suposición.

Si la hipótesis formulada es correcta, entonces el estado x(t) debe verificar la ecuación completa, entonces sustituyendo en la ecuación de estados:

\dot{\Phi}(t,t_0)z(t)+\Phi(t,t_0)\dot{z(t)}=A\Phi (t,t_0)z(t)+Bu(t)

\Phi (t,t_0)\dot{z}(t)+\underset{0}{\underbrace{[\dot{\Phi}(t,t_0)-A(t)\Phi (t,t_0)]}}z(t)=B(t)u(t)

La respuesta completa de un sistema lineal viene dado por:

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

donde x_0=x(t_0) es la condición inicial y \Phi () es una matriz de transición de estados.

Esta respuesta la podemos dividir en dos partes:

x(t)=x_{zi}(t)+x_{zs}(t)

Donde x_{zi}(t) es conocida como respuesta a la entrada cero y x_{zs}(t) como respuesta al estado cero.

En otras palabras, la respuesta a la entrada cero es la respuesta del sistema referente a una condición inicial no nula.

x_{zi}(t)=\Phi (t,t_0)x_0

La respuesta al estado cero es la respuesta impuesta por la señal de entrada u(t) cuando se tienen condiciones iniciales nulas:

x_{zs}(t)=\int_{t0}^{t}\Phi(t,\tau )B(\tau )d\tau

Sin embargo, determinar la matriz \Phi (t,t_0)x_0 no es fácil, pero esta puede ser obtenida utilizando el concepto de matriz fundamental del sistema.

Particularizando el problema para un sistema lineal e invariante en el tiempo (LTI) descrito por:

\dot{x}(t)=Ax(t)+Bu(t)

y(t)=Cx(t)+Du(t)

la expresión analítica de la matriz de transición de estados puede ser determinada utilizando el teorema de Cayley-Hamilton que establece que una matriz satisface su propia ecuación característica.
Recordemos cuanto es la derivada de euler cuando trabajamos con un numero escalra:

\dfrac{\mathrm{d} e^{\lambda t}}{\mathrm{d} t}=\lambda e^{\lambda t}

Esa misma derivada cuando trabajamos con matrices es:

\dfrac{\mathrm{d} e^{\mathbf{A}t}}{\mathrm{d} t}=\mathbf{A} e^{\mathbf{A} t}=e^{\mathbf{A} t}\mathbf{A}

Por definición del e elevado a alguna cosa, es una serie infinita, para el caso matricial esa serie es:

e^{\mathbf{A} t}=\mathbf{I}+\mathbf{A}t+\dfrac{\mathbf{A}^2t^2}{2!}+...

Si multiplicamos la ecuación de estados \dot{x}(t)=Ax(t)+Bu(t) por e^{-\mathbf{A}t} se tiene que:

e^{-\mathbf{A}t}\dot{x}(t)-e^{-\mathbf{A}t}Ax(t)=e^{-\mathbf{A}t}Bu(t)

\dfrac{\mathrm{d} e^{-\mathbf{A}t}x(t)}{\mathrm{d} t}=e^{-\mathbf{A}t}Bu(t)

si se integra desde 0 hasta t

\left [ e^{-\mathbf{A}\tau}x(\tau) \right ]_{\tau=0}^{\tau=t}=\int_{0}^{t}e^{-\mathbf{A}\tau}Bu(\tau)d\tau

e^{-\mathbf{A}t}x(t)-e^{0}x(0)=\int_{0}^{t}e^{-\mathbf{A}\tau}Bu(\tau)d\tau

e^{\mathbf{A}t}\times \left ( e^{-\mathbf{A}t}x(t)=x(0)+\int_{0}^{t}e^{-\mathbf{A}\tau}Bu(\tau)d\tau\right )

por lo tanto se tiene la expresión analitica de los estados es:

x(t)=e^{\mathbf{A}t}x(0)+\int_{0}^{t}e^{-\mathbf{A}(t-\tau)}Bu(\tau)d\tau

y que la señal de salida es:

y(t)=C\left ( e^{\mathbf{A}t}x(0)+\int_{0}^{t}e^{-\mathbf{A}(t-\tau)}Bu(\tau)d\tau \right ) +Du(t)

Asi la matriz de transición de estados para sistemas LTI es dada de la siguiente forma:

\Phi (t,t_0)=e^{\mathbf{A}(t-\tau)}\rightarrow \Phi (t)=e^{\mathbf{A}(t)}

El valor de e^{\mathbf{A}(t)} puede ser determinado de varias formas:

– Utilizando extenciones del teorema de Cayley Hamilton en la forma de determinar las funciones matriciales

– Metodo de respuesta frecuencial.

– Representación en forma frecuencial.

Ejemplo:

Podemos aplicar el teorema de Cayley Hamilton por definición

e^{\mathbf{A} t}=\mathbf{I}+\mathbf{A}t+\dfrac{\mathbf{A}^2t^2}{2!}+\dfrac{\mathbf{A}^3t^3}{3!}++...

O aplicando la transformada de Laplace:

L(e^{\mathbf{A} t})=(sI-A)^{-1}\rightarrow e^{\mathbf{A} t}=L^{-1}((sI-A)^{-1})
\textbf{Forma diagonal}
Puedo representar la matriz en forma diagonal.

Tomando una matriz diagonal de la siguiente forma:

e^{\begin{bmatrix}  \lambda_1 & 0\\  0 & \lambda_2  \end{bmatrix}}=\begin{bmatrix}  1 & 0\\  0 & 1  \end{bmatrix}+\begin{bmatrix}  \lambda_1t & 0\\  0 & \lambda_2t  \end{bmatrix}+\begin{bmatrix}  \frac{\lambda_1^2t^2}{2!} & 0\\  0 & \frac{\lambda_2^2t^2}{2!}  \end{bmatrix}+\cdots=\begin{bmatrix}  1+\lambda_1t+\frac{\lambda_1^2t^2}{2!}+...+ & 0\\  0 & 1+\lambda_2t+\frac{\lambda_2^2t^2}{2!}+...+  \end{bmatrix}=\begin{bmatrix}  e^{\lambda_1t} & 0\\  0 & e^{\lambda_2t}  \end{bmatrix}

Se puede ver que determinar la matriz e^{At} cuando la matriz A es diagonal, es super facil y sencillo, por lo tanto si se consigue llevar la Matriz A a la representación diagonal, va ser relativamente facil determinar la respuesta temporal del sistema.

 

Método de la Respuesta en Frecuencia

 

Si tengo un sistema dado por

\dot{x}=Ax(t)+Bu(t)

Aplicando transformada de Laplace obtengo la siguiente representación:

Sx(s)-x(0)=Ax(s)+Bu(s) x(s)=(s-A)^{-1}x(0)+(s-A)^{-1}Bu(s)

Por inspección, se logra determinar que

\phi(s)=(s-A)^{-1}

Por lo tanto la trayectoria de x(t) es dada por:

L^{-1}(\phi(s))x(0)+L^{-1}(\phi(s)Bu(s))

 

\textbf{Usando Teorema de Cayley Hamilton}

 

Otra forma de calcular el e^{At} que es utilizando el teorema de Cayley Hamilton, que basicamente nos dice:

 

“Que una matriz es capaz de satisfacer la ecuación característica de ella”.

 

Recordando que la ecuación característica de una matriz (A) es cuando por ejemplo se toma una variable, en el curso usamos “s”, y realizamos la siguiente operacion

det(sI-A)=0. Las raíces de esta ecuación son los autovalores de la matriz (Esta ecuación es llamada de ecuación característica).

¿Pero que quiere decir que la matriz satisface la ecuación característica?

Quiere decir que si por ejemplo yo tengo una ecuación característica s^2+2s+1=0 entonces A^2+2A+I va a dar Cero.

Veamos como se pueden representar polinomios de matrices cuadradas:

A^p=AAAA....A, (p terminos), A^0=I

Suponiendo que f(\lambda) es un polinomio monico de la forma: \lambda^2+5\lambda+6 implica que f(A) estaría definido por:

f(A)=A^2+5A+6=(A+2I)(A+3I)

Note que si aplicamos una transformación de similaridad dado por \hat{A}=Q^{-1}AQ, la ecuación matricial satisface

f(\hat{A})=f(Q^{-1}AQ)=Q^{-1}f(A)Q

El teorema de Cayley Hamilton está definido como:

\Delta(\lambda)=det(\lambda I-A)=\lambda^n+\alpha_1\lambda^{n-1}+...+\alpha_{n-1}\lambda+\alpha_n=0

La ecuación característica de A viene dado entonces por:

\Delta(A)=A^n+\alpha_1A^{n-1}+...+\alpha_{n-1}A+\alpha_nI=0

La idea principal es que yo puedo factorizar una expresión cualquiera en función de la ecuación característica.

Supongamos que queremos calcular A^{100} y que A es una matriz 2×2.

A^{100}=q(A)\Delta(A)+h(A)

En la ecuación de encima, se escribe la matriz como una función matricial multiplicado por la ecuación característica y sumado por un residuo que es de grado (n-1), en este caso 1 (n dimensión de la matriz).

Cual es el valor de la ecuacion caracteristica de A (\Delta(A))? es Cero, segun la ecuacion de Cayley Hamilon.

Con esto puedo llegar en la conclusión que A^{100}=h(A). Osea que va a ser una función de A que como máximo tiene grado n-1 donde n es la dimensión de la matriz.

De esta manera es mucho más fácil calcular una matriz elevado a alguna cosa, porque únicamente voy a necesitar calcular el resto (obviamente en calculadora es fácil hacer todo)

 

Pero como voy a encontrar h(A)?

Usando la expresión del teorema de Cayley Hamilton. Es decir que si yo consigo hacer eso para un autovalor, también debe funcionar.

\lambda^{100}=q(\lambda)\Delta(\lambda)+h(\lambda)

De esa forma encuentro el valor de la expresión h(\lambda) y por lo tanto voy a saber cual es el valor de h(A)

Con esto se llega a la conclusión que cualquier polinomio f(A) puede ser escrito como:

f(A)=\beta_0\mathbf{I}+\beta_1\mathbf{A}+...+\beta_{n-1}\mathbf{A}^{n-1}

Lo que tengo que hacer es intentar encontrar las constantes \beta

Ejemplo:

\mathbf{A}=\begin{bmatrix}  0 & -1\\  1 & 0  \end{bmatrix}

Ecuacion caracteristica:

\Delta(\lambda)=det(\lambda\mathbf{I}-\mathbf{A})^{-1}=\lambda^2-1=(\lambda+1)(\lambda-1)=0

Para determinar \mathbf{A}^{100}

f(\lambda)=\lambda^{100}\rightarrow h(\lambda)=\beta_0+\lambda\beta_1 ,para \lambda=\lambda_i \rightarrow f(\lambda_i)=h(\lambda_i)

(-1)^{100}=\beta_0-\beta_1

(1)^{100}=\beta_0+\beta_1

\mathbf{A}^{100}=\beta_0\mathbf{I}+\beta_1\mathbf{A}

Para Autovalores repetidos voy a tener problemas para calcular el valor de A de esa forma, debido a que necesito crear ecuaciones linealmente dependients, entonces simplemente cuando tenga autovalores repetidos hago lo siguiente>

f(\lambda_i)=h(\lambda_i)

Calculo la derivada

\left.\begin{matrix}  \dfrac{df(\lambda)}{d\lambda}  \end{matrix}\right|_{\lambda=\lambda_i}=\left.\begin{matrix}  \dfrac{dh(\lambda)}{d\lambda}  \end{matrix}\right|_{\lambda=\lambda_i}

\vdots

\left.\begin{matrix}  \dfrac{d^{n-1}f(\lambda)}{d^{n-1}\lambda}  \end{matrix}\right|_{\lambda=\lambda_i}=\left.\begin{matrix}  \dfrac{d^{n-1}h(\lambda)}{d\lambda^{n-1}}  \end{matrix}\right|_{\lambda=\lambda_i}

Ejemplo:

Considere el siguiente sistema

\dot{x}(t)=\begin{bmatrix}  0 & 0 & -2\\  0 & 1 & 0\\  1 & 0 & 3  \end{bmatrix}x(t)

Determine el valor de \dot{x}(t) para x_0=[0,1,0]'

Ecuación caracteristica:

\Delta(\lambda)=det(\lambda\mathbf{I}-\mathbf{A})^{-1}=(\lambda-1)^2(\lambda-2)=0

Para determinar e^{At} se puede utilizar el polinomio minimo de h(A).

f(\lambda)=e^{\lambda t}=h(\lambda)=\beta_0+\beta_1\lambda+\beta_2\lambda^2

f(1)=e^{ t}=h(1)=\beta_0+\beta_1+\beta_2

\dfrac{df(1)}{d\lambda}=\dfrac{dh(1)}{d\lambda}=te^{ t}=\beta_1+2\beta_2

f(2)=e^{ 2t}=h(2)=\beta_0+2\beta_1+4\beta_2

Coeficientes de \beta_i

\beta_0=-2te^{t}+e^{2t}, \beta_1=3te^{t}+2e^{t}-2e^{2t}, \beta_2=e^{2t}-e^{t}-te^{t}

Calculando \Phi(t)=e^{At}

\Phi(t)=e^{A t}=\beta_0\mathbf(I)+\beta_1\mathbf(A)+\beta_2\mathbf(A)^2

=\begin{bmatrix}  2e^{t}-e^{2t} & 0 & 2e^{t}-2e^{2t}\\  0 & e^{t} & 0\\  e^{2t}-e^{t} & 0 & 2e^{2t}-e^{t}  \end{bmatrix}

Solución:

x(t)=\Phi(t)x(0)=[0,e^{t},0]

Matriz de jordan 1

Matriz de Jordan 2

Matriz de Jordan 3

Forma de Jordan

Cuando tenemos autovalores repetidos, podemos obtener una representación muy próxima de la diagonal conocida como la Matriz de Jordan donde la construcción de esta matriz se basa en el concepto de autovectores Generalizados.

Representación de Jordan

 

Formas de Jordan

Veamos un video con el funcionamiento de las formas de Jordan:

 

EJEMPLO

Determinar la respuesta temporal del siguiente sistema por la forma diagonal, por cayley hamilton y por respuesta frecuuencial.

\dot{x}(t)=\begin{bmatrix}  0 & 1\\  -2 & -3  \end{bmatrix}x(t)+\begin{bmatrix}  0 \\  1  \end{bmatrix}u(t)
y(t)=\begin{bmatrix}  0 & 1  \end{bmatrix}x(t)

para x_1(0)=1, x_2(0)=1 y u(t)=0.

En la entrada pasada habíamos visto que existen diferentes formas de obtener la matriz de transición de estados para poder determinar la respuesta temporal de un sistema dinámico expresado en variables de estado. Y vimos que una de las formas más faciles de determinar esta matriz se daba cuando nuestra matriz dinámica A es una matriz diagonal.

¿Que sucede si mi matriz A no es diagonal?

Cuando mi matriz no es diagonal, se puede utilizar el concepto de transformación de similitud.

Para esto se parte del sistema original dado por:

\dot{x}(t)=Ax(t)+Bu(t), x(0)

y(t)=Cx(t)+Du(t)

Para poder obtener el A diagonal se utilizan dos conceptos:

1. Se utiliza la matriz de autovectores.

2. En el caso de autovalores repetidos se utiliza el concepto de autovectores generalizados que termina dando la forma de Jordan.

La representación de la matriz de Jordan, es una representación casi diagonal, cuando la matriz tiene autovalores repetidos.

Si consideremos que A_{diagonal}=\hat{A} y V como la matriz de autovectores, procedemos de la siguiente forma:

AV=V\hat{A}

\hat{A}=V^{-1}AV

Deseo obtener una nueva representación de mi sistema donde la matriz dinámica sea \hat{A}. Para eso multiplico por la izquierda el sistema por V^{-1}.

V^{-1}\dot{x}=V^{-1}Ax+V^{-1}Bu

De esa forma puedo renombrar los estados como:

\hat{x}=V^{-1}x

\dot{\hat{x}}=V^{-1}\dot{x}

x=V\hat{x}

Asi mi nuevo estado se expresa como:

\dot{\hat{x}}=V^{-1}AV\hat{x}+V^{-1}Bu

Mi nueva representación donde la matriz A es diagonal es:

\dot{\hat{x}}=\hat{A}\hat{x}+\hat{B}u

y(t)=CV\hat{x}+Du

donde

\hat{B}=V^{-1}B y \hat{C}=CV

Es necesario encontrar el vector de autovectores para poder calcular B.

\dot{\hat{x}}=\hat{A}\hat{x}+\hat{B}u, \hat{x}(0)=V^{-1}x(0)

y(t)=\hat{C}\hat{x}+Du

Resolviendo ese sistema consigo llegar en

\hat{x}(t)\rightarrow x(t)=V\hat{x}(t)

Si por alguna razón necesitara calcular la salida del sistema, no es necesario volver a la representación del sistema original. Note que la salida siempre es la misma en ambos sistemas.

 

Resolución en la Forma Diagonal

Procedimiento: Calcular ecuación característica, los autovalores, los autovectores, la forma diagonal y luego resolver.

La entradas y salidas del sistema en la representacion diagonal es la misma que el sistema original.

G(s)=C(sI-A)^{-1}B+D

 =\hat{C}(sI-A)^{-1}\hat{B}+D

La ecuación caracteristica es:

\Delta(\lambda)=(\lambda I-A)=\lambda^2+3\lambda+2=(\lambda+1)(\lambda+2)=0

Para conseguir los autovectores de una matriz se debe cumplir que:

Suponiendo que los autovalores de una matriz son distintos, entonces van a existir n autovectores de forma que:

Av_i=\lambda_iv_i,i=1,2,...,n

Y en ese caso los autovectores son todos distintos y linealmente independientes, asi pueden ser utilizados para construir una base con la siguiente matriz de transformación:

V=[v_1,v_2,...,v_n]

Los vectores que satisfacen la relacion:

(A-\lambda_i)v_i=0, i=1,2,...n

Son conocidos como autovectores de A.

Retomando el ejemplo

Primer autovalor

(A-(-1)I)v_1=0

(A+I)\begin{bmatrix}v_{11} \\v_{12}\end{bmatrix}=\begin{bmatrix}0 \\0\end{bmatrix}

Segundo autovalor

(A-(-2)I)v_2=0

(A+2I)\begin{bmatrix}v_{21} \\v_{22}\end{bmatrix}=\begin{bmatrix}0 \\0\end{bmatrix}

De estas dos representaciones obtengo una ecuación que es una recta.

Solucionando para el primer autovector

\begin{bmatrix}1&1 \\-2&-2\end{bmatrix}\begin{bmatrix}v_{11} \\v_{12}\end{bmatrix}=\begin{bmatrix}0 \\0\end{bmatrix}

Donde se obtiene:

v_{11}+v_{12}=0

Puedo escoger dos valores que cumplan esa igualdad (soluciones infinitas), entonces:

v_1=\begin{bmatrix}1 \\-1\end{bmatrix}

Solucionando para el segundo autovector

\begin{bmatrix}2&1 \\-2&-1\end{bmatrix}\begin{bmatrix}v_{21} \\v_{22}\end{bmatrix}=\begin{bmatrix}0 \\0\end{bmatrix}

Donde se obtiene:

2v_{21}+v_{22}=0

Puedo escoger dos valores que cumplan esa igualdad (soluciones infinitas), entonces:

v_2=\begin{bmatrix}1 \\-2\end{bmatrix}

El auto vector de A es:

v=\begin{bmatrix}1 &1\\-1&-2\end{bmatrix}

La nueva matriz A es

\hat{A}=v^{-1}Av

\hat{A}=\begin{bmatrix}-1&0 \\0&-2\end{bmatrix}

El nuevo vector B seria:

\hat{B}=v^{-1}B

\hat{B}=\begin{bmatrix}1 \\-1\end{bmatrix}

El nuevo vector C seria:

\hat{C}=Cv

\hat{C}=\begin{bmatrix}-1 &-2\end{bmatrix}

La nueva condicion inicial es

\hat{x}(0)=v^{-1}x(0)

\hat{x}(0)=\begin{bmatrix}3 \\-2\end{bmatrix}

El sistema diagonal queda dado por:

\dot{\hat{x}}(t)=\begin{bmatrix}-1&0 \\0&-2\end{bmatrix}x(t)+\begin{bmatrix}1 \\-1\end{bmatrix}u(t)

y(t)=\begin{bmatrix}-1 &-2\end{bmatrix}x(t)

La salida del sistema es:

y(t)=\hat{C}\hat{x}=\hat{C}\begin{bmatrix}e^{-1t}&0 \\0&e^{-2t}\end{bmatrix}\begin{bmatrix} 3 \\-2\end{bmatrix}= 4e^{-2t}-3e^{-t}

RESOLUCIÓN POR CAYLEY HAMILTON

La ecuación caracteristica es:

\Delta(\lambda)=(\lambda I-A)=\lambda^2+3\lambda+2=(\lambda+1)(\lambda+2)=0

Se determina h(\lambda) reemplazando los autovalores y hallando el sistema de ecuaciones

f(\lambda)=e^{\lambda t}=h(\lambda)=\beta_0+\beta_1\lambda

1. e^{-t}=\beta_0-\beta_1

2. e^{-2t}=\beta_0-2\beta_1

soluciono el sistema de ecuaciones

\beta_1=e^{-t}-e^{-2t}

\beta_0=2e^{-t}-e^{-2t}

De esa forma

h(\lambda)=(2e^{-t}-e^{-2t}) + (e^{-t}-e^{-2t})\lambda

h(A)=e^{At}=(2e^{-t}-e^{-2t})I + (e^{-t}-e^{-2t})A

e^{At}=\begin{bmatrix}2e^{-t}-e^{-2t}&0 \\0&2e^{-t}-e^{-2t}\end{bmatrix}+  \begin{bmatrix}0&e^{-t}-e^{-2t} \\-2(e^{-t}-e^{-2t})&-3(e^{-t}-e^{-2t})\end{bmatrix}

e^{At}=\begin{bmatrix}2e^{-t}-e^{-2t}&e^{-t}-e^{-2t} \\-2e^{-t}+2e^{-2t}&-e^{-t}+2e^{-2t}\end{bmatrix}

Los estados son:

x(t)=e^{At}\begin{bmatrix}1 \\1\end{bmatrix}=\begin{bmatrix}3e^{-t}-2e^{-2t} \\-3e^{-t}+4e^{-2t}\end{bmatrix}

y=-3e^{-t}+4e^{-2t}

RESOLUCIÓN DE LA FORMA FRECUENCIAL

\Phi(s)=(sI-A)^{-1}

\Phi(s)=\begin{bmatrix}s&-1 \\2&s+3\end{bmatrix}^{-1}

\Phi(s)=\begin{bmatrix}\dfrac{s+3}{s^2+3s+2}&\dfrac{1}{s^2+3s+2} \\ \dfrac{2}{s^2+3s+2}&\dfrac{5}{s^2+3s+2}\end{bmatrix}

\Phi(t)=L^{-1}(\Phi(s))=\begin{bmatrix}2e^{-t} - e^{-2t}& e^{-t} - e^{-2t} \\ 2e^{-2t} - 2e^{-t} & 2e^{-2t} - e^{-t}\end{bmatrix}

Multiplico por el vector de estados iniciales

x(t)=\Phi(t)x(0)

x(t)=\begin{bmatrix}3e^{-t} - 2e^{-2t} \\ 4e^{-2t} - 3e^{-t}\end{bmatrix}

La salida del proceso con la condicion inicial es:

y(t)=4e^{-2t} - 3e^{-t}

 

 

 

 

Si te sirvió de algo la información de la pagina, podrías invitarme a un café y ayudarme a seguir manteniendo en pie el sitio WEB. Solo cuesta $2USD y me ayudarías enormemente a seguir con mi trabajo. Muchas Gracias por tu visita.