Trabajo Final_cmv.docx

UNIVERSIDAD DE SANTIAGO DE CHILE FACULTAD DE INGENIERIA Departamento de Ingeniería Eléctrica

CONTRO MULTIVARIABLE LMI TRABAJO FINAL “Generador de Vapor Planta Refinería”

Nombre: Francisco Ramos Ccalla Profesor(a): Karina Acosta Barbosa

Santiago – Chile 2019

INDICE 1. RESUMEN 2. INTRODUCCION 3. DESCRIPCION DE LA PLANTA Y MODELO DEL SISTEMA 4. ANALISIS DE ESTABILIDAD 5. DISEÑO DE CONTROLADORES 5.1 POR REALIMENTACIÓN DE ESTADOS 5.2 POR REGIÓN LMI 5.3 OBSERVADOR DE ESTADO ORDEN COMPLETO 5.4 CONTROLADOR POR REALIMENTACION DE ESTADOS OBSERVADOS 6. NORMA INFINITA: PERTURBACIÓN SISTEMA 7. CONCLUSIONES 8. REFERENCIAS 9.

2

1. RESUMEN

En el presente informe se resuelven los problemas de diseños de varios tipos de controladores aplicados a sistemas físicos representados en sus espacios de estado. Para abordar el problema de diseño, se toma el enfoque de optimización usando LMIs como restricciones del problema de factibilidad u minimización, según corresponda.

3

2. INTRODUCCIÓN

El enfoque tradicional de sistemas de control se basa en el estudio de los sistemas en el dominio de la frecuencia. Para sistemas de mayor complejidad que involucran mayor cantidad de variables de entrada y salida, técnicas de diseño multivariables son más efectivas y fáciles de implementar con ayuda de programación en software. El diseño de controladores de sistemas MIMO involucra el algebra matricial para lograr la ubicación de polos, estabilización, retroalimentación de estados, etc. Estos métodos son de fácil implementación si se logra encontrar un conjunto de desigualdades matriciales lineales (LMI), que transforman el diseño en un problema de minimización o factibilidad. En el presente informe se aplican condiciones LMI para mostrar la conveniencia de estos métodos.

4

3. OBJETIVOS



Aplicar correctamente técnicas LMIs para el diseño de controladores por realimentación de estados.



Identificar las principales características de los controladores por realimentación de estados.

.

5

4. DESCRIPCIÓN DE LA PLANTA Y MODELO DEL SISTEMA

La planta en estudio es un generador de vapor de la planta de potencia integrada en la refinería de Syncrude Canada Limited (SCL) ubicada en Mildred Lake, Alberta. La planta suministra electricidad a la red de distribución eléctrica canadiense como al propio proces de refinería. El generador analizado es llamado “Generador de Servicio”.

Un liquido precalentado en un economizador se alimenta al colector superior que debe mantenerse a 6,306 MPa de presión y a cierto nivel de agua. (Figura 1). El flujo de agua de alimentación cae al colector de sedimentos, donde se calienta hasta condiciones de saturación. El horno genera calor con el flujo de combustible y aire de entrada controlado por un extractor. La mezcla de líquido y vapor sube por la tubería de elevación (riser), de modo que se tiene en el colector superior vapor y agua. El agua no evaporada regresa al colector de sedimentos por los tubos descendentes (downcomers), mientras que el vapor es separado y dirigido a los recalentadores primario y secundario donde se distribuye a los colectores de menor presión (no mostrados en la figura), con la temperatura de vapor necesaria. Entre los recalentadores se ubica un atemperador, el cual regula la temperatura de vapor que sale del recalentador secundario por medio de inyección de agua de menor temperatura.

Modelo del sistema

Las variables de entrada al sistema son las siguientes: 

𝑢1: Caudal de agua de alimentación [kg/s]



𝑢2: Caudal de combustible [kg/s]



𝑢3: Caudal de spray de atemperador [kg/s] Se tiene el vector de entradas 𝑢1 = (𝑢1, 𝑢2, 𝑢3)𝑇

Las variables de salida del sistema son las siguientes: 

𝑦1: Nivel del colector [m]



𝑦2: Presión del colector [MPa]



𝑦3: Temperatura de vapor [grados C]

6

Se tiene el vector de salidas 𝑦 = (𝑦1, 𝑦2, 𝑦3)𝑇

Figura 1: Generador de Vapor. Planta bajo estudio

Al ser la planta de naturaleza no lineal, se espera que durante la operación normal las variables de salida se mantengan a valores predeterminados. Bajo este supuesto, se esperan desviaciones relativamente pequeñas con respecto a los valores nominales y un modelo LTI de segundo orden para la planta puede ser planteado. La matriz de funciones de transferencias es de orden 3x3: 𝐺11 𝐺𝑝 = (𝐺21 𝐺31

𝐺12 𝐺22 𝐺32

𝐺13 𝐺23 ) 𝐺33

donde los 𝐺𝑖𝑗 son las siguientes:

7

𝐺11 =

(−1.6𝑠 2 + 5.2 ∗ 10−1 𝑠 + 1.4 ∗ 10−2 ) ∗ 10−4 ) 𝑠(𝑠 + 1.68 ∗ 10−2 )

𝐺12 =

(3.1𝑠 − 3.2 ∗ 10−2 ) ∗ 10−3 𝑠(𝑠 + 2.15 ∗ 10−2 )

𝐺12 =

(3.1𝑠 − 3.2 ∗ 10−2 ) ∗ 10−3 𝑠(𝑠 + 2.15 ∗ 10−2 ) 𝐺13 = 0

𝐺23 =

𝐺21 = −

3.95 ∗ 10−5 (𝑠 + 1.8 ∗ 10−2 )

𝐺22 = −

2.51 ∗ 10−3 𝑠 + 1.57 ∗ 10−2

(5.88𝑠 2 + 2.015𝑠 + 9 ∗ 10−3 ) ∗ 10−4 (𝑠 2 + 3.52 ∗ 10−2 𝑠 + 1.42 ∗ 10−4 )

𝐺31 =

(−1.18𝑠 + 1.39 ∗ 10−1 ) ∗ 10−3 ) (𝑠 2 + 1.852 ∗ 10−2 𝑠 + 9.1 ∗ 10−5 )

𝐺32 =

𝐺33 =

(4.48𝑠 + 1.1 ∗ 10−2 ) ∗ 10−1 ) (𝑠 2 + 1.27 ∗ 10−2 𝑠 + 9.5 ∗ 10−5 ) (5.82𝑠 − 2.43 ∗ 10−1 ) ∗ 10−1 ) + 1.076 ∗ 10−1 𝑠 + 1.04 ∗ 10−3 )

(𝑠 2

Para el presente trabajo, se aborda el sistema con su modelo de espacio de estados. Utilizando las herramientas en MATLAB, se obtiene el siguiente modelo de espacio de estados

𝑥̇ = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) 𝑦 = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡)

8

Ya que existen valores de las matrices que son muy pequeños (del orden de 10−30 ), existen variables de estado cuyo aporte al comportamiento del sistema es despreciable. En base al modelo reducido (reducción de Hankel, por valores singulares), el conjunto de matrices es: −0.0112 −0.0057 −0.0027 𝐴= −0.0055 5.64 ∗ 10−5 ( 0

0.0112 −0.0020 0.0023 −0.0045 6.2877 ∗ 10−5 0

0.0043 −0.0025 −0.0039 0.0191 −3.2633 ∗ 10−4 0

0.0054 0.0191 −0.0038 𝐵= 0.0065 0.0365 (−0.0048

0.0033 𝐶 = (0.0191 0.6933

−0.0034 −0.0026 −0.2344

0.0044 0.0025 −0.1250

0.0116 −0.0059 −0.0282 −0.1029 0.0013 0

0.6132 −0.1464 0.1117 −0.0071 0.0017 0.0857

0 0 0 0 0 −1.9407 ∗ 10−18 )

−0.3236 −0.1822 0.0562 −0.3546 0.0074 1.737 ∗ 10−16 )

−5.5581 ∗ 10−4 −6.5321 ∗ 10−4 −0.3547

−1.6000 ∗ 10−14 𝐷=( 0 0

−0.0012 6.37 ∗ 10−4 0.0023 0.0214 −0.0092 0

0 0 0

0.0018 0.0010 0.0372

−0.0174 4.5416 ∗ 10−21 ) 7.0703 ∗ 10−35

0 5.88 ∗ 10−4 ) 0.5820

Ya que el modelo se obtuvo considerando su funcionamiento en torno al punto de operación nominal, el modelo se basa en variables de desviación de entradas y salidas en torno a los siguientes valores nominales de entrada y salida: 40.68 𝑢1∗ 𝑢 = (𝑢2∗ ) = (2.102) 0 𝑢3∗ ∗

𝑦1∗ 1,0 𝑦 = (𝑦2∗ ) = (6,306) 𝑦3∗ 466.7 ∗

9

5. ANALISIS DE ESTABILIDAD Considere el sistema lineal en lazo abierto 𝑥̇ = 𝐴 ∗ 𝑥(𝑡)

La matriz A es asintóticamente estable (Hurwitz-estable) si todos sus valores propios tienen parte real negativa. Si el sistema cumple con esta condición, el vector de estados tenderá al estado cero al trascurrir el tiempo. Se verifica la estabilidad con las siguientes restricciones LMI, basada en el teorema de estabilidad de Lyapunov

El sistema es Hurwitz-estable si y sólo si existe una matriz P simétrica tal que: 𝑃>0 𝑇

𝐴 𝑃 + 𝑃𝐴 < 0

El script ejecuta en MATLAB (vea Anexo) encuentra una solución P al problema de factibilidad, por tanto, se puede concluir que es asintóticamente estable. Se tabulan los autovalores de A:

𝜆1 = −0.0969 𝜆2 = −0.0064 + 0.0074𝑖 𝜆3 = −0.006 − 0.0074𝑖 𝜆4 = −0.0107 𝜆5 = −0.0089 𝜆6 ≈ 0

A pesar de haber encontrado solución al problema, existe un autovalor muy cercano al origen, esperable debido a que la matriz de transferencias presentaba integradores que implican problemas con la estabilidad. Sin embargo, el diseño de todos los controladores tratados aquí se cumplen la función de estabilizar el sistema, por lo cual no representará mayor problema.

10

DISEÑO DE CONTROLADORES 6. POR REALIMENTACIÓN DE ESTADOS.

La matriz A es asintóticamente estable (Hurwitz-estable) si todos sus valores propios tienen parte real negativa. Si el sistema cumple con esta condición, el vector de estados tenderá al estado cero al trascurrir el tiempo. Se verifica la estabilidad con las siguientes restricciones LMI, basada en el teorema de estabilidad de Lyapunov

𝑥̇ = (𝐴 + 𝐵 ∗ 𝐾)𝑥(𝑡)

son las siguientes: 𝐴𝑃 + 𝑃𝐴𝑇 + 𝐵𝑊 + 𝑊 𝑇 𝐵𝑇 < 0

Si si puede encontrar un matriz simétrica definida positiva P y una matriz W que satisfagan la desigualdad anterior, entonces existe una ley de control 𝑢(𝑡) = 𝐾𝑥(𝑡) tal que el lazo de control cerrado es asintóticamente estable, y la ganancia de retroalimentación tiene la forma:

𝐾 = 𝑊𝐵 −1

%Lazo realimentacion estabilizacion. u=W*P^-1 %Dimensiones matrices A, B n=size(A,1); [n1,m]=size(B); %Definicion de variables de decision 11

P=sdpvar(n,n); W=sdpvar(n1,m); %Seleccion de solver para optimizacion options = sdpsettings('solver','sdpt3'); % para usar el sdpt3 %Definicion de restricciones LMI y solucion del problema LMIs= [A*P+P*A’+B*W+W’*B’<=0,P>=0]; Solution = optimize(LMIs,[],options) check_LMI=check(LMIs) %Bloque para verificar la validez de soluciones encontradas. if check_LMI > 0, P = double(P); W = double(W); disp('las matrices soluciones son') P; W; disp('la matriz ganancia de realimentacion es: ') K=W*(P^-1) else disp('non feasible') end

Ahora, se necesita un script para garantizar que los polos de lazo cerrado estén en la región LMI que se define como sigue: 𝑆(𝛼, 𝑟, 𝜃) = {(𝑥, 𝑦), 𝑥 < −𝛼 < 0, |𝑥 + 𝑗𝑦| < 𝑟, |𝑦| < −𝑥𝑡𝑎𝑛(𝜃)}

La figura 1 muestra la región en el plano complejo. Esta región es una generalización de las regiones D y las regiones H. La ventaja de diseñar un controlador que ubique los polos dentro de la región S, es el control sobre los límites de parámetros como el amortiguamiento, frecuencia natural amortiguada, tiempo de subida, entre otros.

12

Figura 1: Región LMI punto (b).

El script para garantizar la estabilidad del dentro de esta región es el siguiente:

%Lazo realimentacion estabilizacion. u=W*P^-1 n = size(A,1); [n1,m]=size(B); %Dimensiones matrices A, B P= sdpvar(n,n); W= sdpvar(n1,m); %Seleccion de solver para optimizacion alfa=(1/50); r=2; theta=((45*2*pi)/360); a=sin(theta); b=cos(theta); %Seleccion de solver para optimizacion options = sdpsettings('solver','sdpt3'); % para usar el sdpt3 %Definicion de restricciones LMI y solucion del problema const1= 2*alfa*P+A*P+P*A'+B*W+W'*B'; const2 = [-r*P A*P+B*W P*A'+W'*B' -r*P]; const3 = [(A*P+P*A+B*W+W'*B')*a (A*P-P*A+B*W-W'*B')*b (P*A'-A*P+W'*B'-B*W)*b (A*P+P*A'+B*W+W'*B')*a]; 13

LMIs = [const1<=0,const2<=0,const3<=0,P>=0]; %Bloque para verificar la validez de soluciones encontradas. prueba = optimize(LMIs,[],options) check_LMI=check(LMIs) if check_LMI > 0, P = double(P); W = double(W); disp('la matriz P,W solucion:') P W disp('la matriz ganancia de realimentacion es: ') K=W*(P^-1) else disp('non feasible') end

14

ACTIVIDADES 7. ACTIVIDAD 2

Considere el siguiente sistema dinámico multivariable utilizado en la experiencia anterior:

𝑥̇ = (

−0.06120 0

0 0.0187 ) 𝑥(𝑡) + ( −0.0408 0.0281

−0.0281 ) 𝑢(𝑡) −0.0281

Diseñe una ley de control de manera que el sistema en lazo cerrado sea: a) Estable b)

𝑆{1/50,2,45}

c)

𝑆{1/20,4,45}

P

W

Autovalores

K

(A+B*K) (a)

(b)

(c)

(

(

(

10.0009 −0.0004

−0.0004 ) 10.006

10.7494 −1.6863

−1.6863 ) 9.5867

7.9922 −1.0926

−1.0926 ) 5.3884

3.3198 2.3653 ( ) ∗ 10−3 2.3653 2.5287

(

47.4875 94.47894 (

55.4765 57.9174

331.9603 236.5219

236.5326 ) 252.8635

−0.4998

94. .47794 ) 166.3303

6.1329 ( 11.8376

10.9340 ) 19.4323

−0.0454

57.9174 ) 84.7603

(

8.6506 9.6652

12.5025 ) 17.6898

−0.0918

(

± 2.6821𝑖

−0.5134

−0.2658

15

Figura

2:

Respuesta

a

entradas

sistema

(b)

(p=50,

v=40)

Figura

3:

Respuesta

a

entradas

sistema

(c)

(p=50,

v=40)

16

Figura 4: Respuesta a entradas sistema (a) (p=50, v=40)

Al calcular los autovalores de lazo cerrado, se observa que los polos de (b) y (c) se encuentran dentro de la región LMI deseada. Todos los controladores estabilizan las salidas en dentro de los 20-30 segundos de haber aplicado las entradas. Para el controlador en (a), solo se garantiza la estabilidad y no la zona donde se ubiquen los polos. Esto se observa calculando el ángulo de los polos 𝜃 = 80 𝑔𝑟𝑎𝑑𝑜𝑠, aproximadamente, polos que no pertenecen a las zonas LMI de (b) y (c). Las respuestas (b) y (c) poseen amortiguamiento igual a 1, mientras que el controlador en (a) genera una respuesta subamortiguada, debido a la parte imaginaria de sus polos. Por último, se comprueba que la respuesta (c) se estabiliza más rápidamente que la respuesta (b), debido a que posee mayor frecuencia, por tanto velocidad (r=4, en comparación a r=2 en (b))

17

ACTIVIDADES 8. ACTIVIDAD 3

Para el sistema anterior diseñe un controlador del tipo LQR, donde

a) Q=I y R=[1 0; 0 15] b) Q=I y R=[15 0; 0 1] c) Q=[9 0; 0 0,1] y R= I d) z= 3𝑥1 + √15𝑢1



Para cada ítem, simule el sistema en lazo cerrado para una condición inicial no nula, y haga un análisis comparativo de las señales de salida y de control.



En cada ítem explique cuál es el objetivo del control.

Q (a) (b) (c)

1 0 ) 0 1 1 0 ( ) 0 1 90 ( ) 0 0.1 (

R (

1 0

0 ) 15

15 0 ) 0 1 1 0 ( ) 0 1

(

K 0.1369 −0.0141 0.0086 ( −0.2027 0.9287 ( −1.4025 (

0.3025 ) −0.0199 0.0201 ) −0.2951 0.0247 ) −0.0200

18

Figura

5:

Figura

6:

Respuesta

Respuesta

a

a

entradas

entradas

sistema

sistema

(a)

(b)

(p=50,

v=40)

(p=50,

v=40)

19

Figura 7: Respuesta a entradas sistema (c) (p=50, v=40)

Las respuestas (a) y (b) son similares, ya que la matriz Q es la misma, y los estados en esta planta son las salidas de las alturas de los depósitos. La matriz R en (a) pondera más la entrada de potencia que la potencia de balance, mientras que la matriz R en (b) pondera más la entrada de balance de depósitos. El objetivo de control en (a) es poder ahorra potencia de entrada y concentrar la acción de control en el balance que pueden hacer las válvulas. El objetivo de control en (b) busca la misma acción sobre cada salida, pero ahora el cambio se lograr con mayor potencia por parte de la bomba.

En la respuesta (c) se cambia el objetivo de control. Ahora, se cambias las ponderaciones en los estados, es decir, las salidas. Se conserva la misma ponderación para las entradas. El objetivo es entonces darle prioridad al cambio en la altura del depósito h4. Observe que el cambio en el h4 es mucho más pronunciado que la salida h1, el cual tiene un factor de 0.1, que es muy bajo respecto al factor de 9 de h4.

20

9. CONCLUSIONES Se puede diseñar una serie de controladores que garantizan la estabilidad del sistema en lazo cerrado. Para cada uno de los controladores diseñados se garantiza que el sistema será asintóticamente estable. Si se requiere ubicar los polos para obtener ciertas características deseadas en la respuesta, se opta por un estabilizador que ubique los polos en la región S vista en las actividades. Si bien este controlador no da una ubicación exacta, garantiza que los parámetros de la respuesta este acotado a ciertos límites según se ha vista. Cuando se quiere priorizar o penalizar las acciones de control y su efecto sobre los estados, el controlador LQR es el adecuado. Las matrices Q y R le otorgan una importancia relativa de cada entrada y estado con respecto a otra. La minimización del parámetro de desempeño J permite obtener un equilibrio entre la acción de control y la acción sobre los estados, según el objetivo de control que se este buscando en el problema Al igual que en un SISO, es posible realizar un seguimiento de la referencia con el controlador con un integrador incorporado, de modo de garantizar el error en estado estacionario cero. Cada uno de los controladores son fácilmente implementables con técnicas LMIs, y debe decidirse cuál es la más adecuada para cada aplicación.

21

22