Control robusto de un sistema de dos rotores en cuadratura

CONTROL ROBUSTO DE UN SISTEMA DE DOS ROTORES EN CUADRATURA M.L´opez-Mart´ınez, M.G.Ortega y F.R.Rubio Dpto.Ingenier´ıa de Sistemas y Autom´atica. Universidad de Sevilla. Escuela Superior de Ingenieros. Tel´efono: +34 95 4487361. Fax: +34 95 4487340 correo electr´onico: {mlm,ortega,rubio}@cartuja.us.es Resumen En este art´ıculo se presenta el dise˜no de un control H∞para un sistema de dos rotores en cuadratura. Este sistema se caracteriza por ser multivariable (de 2 entradas y 2 salidas) altamente no lineal y fuertemente acoplado. Se ha hecho uso de la teor´ıa de H∞para sistemas lineales, aplic´andola al sistema linealizado mediante realimentaci´on. Palabras clave: Control Robusto, Control H∞, Linealizaci´on por realimentaci´on, Control Multivariable, Doble Rotor. 1 INTRODUCCI ´ ON En este art´ıculo se presenta un controlador robus- to dise˜nado para un sistema no lineal multivariable y fuertemente acoplado, que est´a sometido a perturbaciones. Dicho sistema ha sido linealizado por realimentaci´on mediante la t´ecnica cl´asica del control por par calculado utilizado en al ´ambito de la rob´otica. Debido a las incertidumbre existentes en el modelo, se ha observado que el desacoplamiento del sistema no se produce en alta frecuencia, existiendo una incertidumbre en la din´amica del sistema linealizado. Este hecho justifica el dise˜no de un controlador robusto, eligiendo la teor´ıa de control H∞gracias a sus buenas propiedades de control tanto desde el punto de vista de robustez como de comportamiento. El art´ıculo est´a estructurado de la siguiente manera: en el apartado 2 se describe el sistema. Seguidamente, en la secci´on 3 se realizan ciertas consideraciones sobre la linealizaci´on por par calculado, describiendo la identificaci´on llevada a cabo en alta frecuencia del sistema linealizado. A continuaci´on, en el apartado 4 se muestra c´omo se ha realizado el dise˜no del controlador H∞. En la secci´on 5 se presentan los distintos resultados. Finalmente, en el apartado 6 se muestran las conclusiones. 2 DESCRIPCI ´ ON DEL SISTEMA El equipo empleado es un artilugio mec´anico de dos grados de libertad (ver figura 1) impulsado por el movimiento de dos rotores en cuadratura ([4]). Los grados de libertad son cabeceo y gui˜nada (elevaci´on y orientaci´on). Este equipo se caracteriza por ser multivariable, no lineal y fuertemente acoplado. Figura 1: Doble rotor en cuadratura. El modelado del sistema, as´ı como el valor de los par´ametros utilizados pueden ser encontrados en [5]. 3 LINEALIZACI ´ ON POR REALIMENTACI ´ ON Utilizando el modelo no lineal del sistema desarrollado [5], se ha realizado una linealizaci´on por realimentaci´on utilizando la t´ecnica del control por par calculado empleada com´unmente en el ´ambito de la rob´otica ([10]), obteni´endose a su vez un desacoplamiento parcial de las entradas y salidas del sistema linealizado. Debido a las caracter´ısticas de la linealizaci´on llevada a cabo y al ruido que aparece en los sensores a partir de los cuales realimentamos, se observa que el desacoplamiento es tanto peor conforme se Referencia Linealización por Par Calculado τ Adaptación ω Velocidades de las hélices Posiciones y Velocidades en Orientación y Elevación Control LQR Posiciones y Velocidades en Orientación y Elevación Sistema Doble Rotor u Control en Velocidad de hélices A B Figura 2: Linealizaci´on por Par calculado y sistema a identificar. va aumentando la frecuencia de la se˜nal de referencia del sistema. Por tanto, existe una incertidumbre en el grado de desacoplamiento del sistema, el cual depende de la frecuencia de las se˜nales, siendo mayor dicha incertidumbre conforme aumenta la frecuencia de trabajo del sistema. A fin de poder calcular un controlador robusto, se ha llevado a cabo una identificaci´on del sistema en alta frecuencia para acotar la incertidumbre sobre el modelo. El primer inconveniente que surge a la hora de identificar el nuevo sistema resultante de aplicar el par calculado (sistema de A a B en la figura 2) es que el ´este es inestable en bucle abierto, por lo que no se puede hacer una identificaci´on tradicional. Para poder realizar la identificaci´on, se ha cerrado el bucle de control mediante un controlador LQR (ver figura 2 previamente dise˜nado en [5]). Se ha llevado al sistema a un punto de trabajo y se le ha aplicado una secuencia binaria pseudo aleatoria (PRBS) en la referencia de amplitud constante y de per´ıodo variable entre 0.02 s. y 1s (frecuencia entre 1Hz y 50Hz). Esto se ha realizado tanto para la orientaci´on (gui˜nada) como para la elevaci´on (cabeceo) (ver figura 3). En el dise˜no del PRBS se ha tenido en cuenta que el tiempo de muestreo del sistema es de T=0.01 s, y que el tiempo de subida ante entradas en escal´on es del orden de segundos. Barriendo este rango de frecuencias se garantiza la riqueza de la excitaci´on del sistema para poder llevar a cabo la identificaci´on. Sin embargo, al estar el sistema en bucle cerrado, la entrada depende de la salida a trav´es del controlador situado en el lazo de realimentaci´on. Por esta raz´on, se ha inyectado un ruido blanco adicional en la entrada del sistema a identificar, disminuyendo as´ı la correlaci´on entre las se˜nales de entrada y salida del sistema. 10 15 20 25 30 35 40 −5 −4 −3 −2 −1 0 1 2 3 4 5 PRBSxi Time ξ ( Elevación - Cabeceo ) Figura 3: Secuencia binaria pseudo aleatoria. En la figura 4 se muestra la respuesta del sistema al ser sometido al PRBS junto con el ruido blanco, habiendo llevado previamente al sistema a un punto de trabajo. En este caso, el punto de trabajo corresponde a 25oen elevaci´on (cabeceo) y −90oen orientaci´on (gui˜nada). Tras realizar una serie de experimentos con el sistema linealizado se observa que la respuesta del mismo sigue estando acoplada cuando la frecuencia de la se˜nal de entrada es elevada, esto es, que sucesivos cambios bruscos en la entrada afectan a las dos salidas. Por lo tanto, la identificaci´on del sistema en dicho rango de frecuencias da lugar a una matriz no diagonal de cuatro funciones de transferencia. Usando un algoritmo de identificaci´on multivariable basado en variables de estado de orden 4 ([6]), el modelo obtenido, expresado como matriz de funciones de transferencia, es el siguiente: G(s) = 1 d(s)  n11(s)n12(s) n21(s)n22(s) (1) 12 14 16 18 20 22 24 26 28 −140 −130 −120 −110 −100 −90 −80 −70 −60 −50 −40 theta−Time Time θ ( Orientación - Guiñada ) 12 14 16 18 20 22 24 18 20 22 24 26 28 30 32 xi−Time Time ξ ( Elevación - Cabeceo ) Figura 4: Respuesta del sistema ante el PRBS. n11(s) = (−0.2537s4 −25.640s3 −15.640s2 −64.690s−0.2368) n21(s) = (−0.0001s4 −0.6141s3+43.160s2+125.700s+96.290) n12(s) = (−0.0014s4+ 0.350s3+ 0.5438s2 −6.1820s−3.0020) n22(s) = (−0.2959s4 −6.8510s3 −13.760s2+ 5.750s−4.070) d(s) = 180 π(1.000s4+ 3.8350s3+ 4.8540s2+ 5.5660s+ 1.1450) En la figura 5 se muestra una comparativa entre la respuesta del sistema real y el sistema simulado. Puede observarse como ambas respuestas (la del modelo y la obtenida experimentalmente) son cualitativamente parecidas. Este es el modelo que se va a usar para acotar la incertidumbre del sistema, y a partir de las cuales se obtendr´a un controlador H∞. 4 S´ INTESIS DEL CONTROLADOR H∞ Una configuraci´on general para expresar el problema de control H∞puede mostrarse en el diagrama de bloques representado en la figura 6, donde Prepresenta a una planta generalizada, K representa al controlador, ues la se˜nal de control, ves el vector de medidas, ωes el vector de perturbaciones externas y zes el vector objetivo cuya norma se quiere hacer peque˜na. 30 31 32 33 34 35 36 37 38 39 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 Time Elevación en Grados Salida Real y Simulada Real Simulada Figura 5: Respuesta del sistema real y del identificado ante el PRBS. K(s) z vu ωP(s) Figura 6: Configuraci´on general de control. Llamando Tzω a la matriz de funciones de transferencia desde ωhasta z(ver figura 7), esta vendr´a dada por la transformaci´on lineal fraccional inferior z=Tzωω=Fl(P, K)ω, donde Fl(P, K) = Pzω +PzuK(I−PvuK)−1Pvω. Seg´un el conocido Teorema de la peque˜na ganancia ([11]), una condici´on suficiente para la estabilidad robusta del sistema en bucle cerrado es que se cumpla kTzwk∞k4k∞<1, donde k·k∞representa a la norma infinito. Cuando esta es aplicada a una funci´on de transferencia G, su expresi´on es kGk∞≡supω¯σ(G(jw)) donde ¯σ(·) representa el valor singular m´aximo. ∆(s) Tzω(s) y∆ u∆ Figura 7: Configuraci´on para Teorema de la peque˜na ganancia. El algoritmo de s´ıntesis utilizado para el c´alculo del controlador de esta aplicaci´on esta basado en la soluci´on presentada en [3], y que esta implementado en diversos paquetes de programas comerciales, como por ejemplo [1, 2]. A continuaci´on se exponen los pasos que permitir´an llevar a cabo la s´ıntesis del controlador: 4.1 CONFIGURACI ´ ON DE LA PLANTA AUMENTADA En la figura 8 se presenta un esquema de un sistema de control en bucle cerrado, en el que est´an presentes perturbaciones, d, y ruidos en los sensores de medida. y d n u e r Referencia G Planta K Controlador Figura 8: Sistema en bucle cerrado. Teniendo en cuenta las definiciones: L(s) = G(s)·K(s) (2) S(s)=(I+L(s))−1(3) T(s) = L(s)·S(s) = I−S(s) (4) podemos expresar la salida, el error y la se˜nal de control como: y=T·r+S·d−T·n(5) e=S·(r+d)−T·n u=K·S·(r−n−d) Para conseguir que la salida siga la referencia, interesa que T(s) tenga sus valores singulares en torno a la unidad en baja frecuencia, y que ´estos sean peque˜nos en alta frecuencia. Por otro lado, interesa que S(s) tenga su valores singulares peque˜nos en baja frecuencia. De esta forma tambi´en se consigue que el error tienda a cero. En cuanto a la se˜nal de control, interesa que K(s)·S(s) sea tal que no produzca saturaci´on en los actuadores. Para moldear S(s) y T(s), y calcular un controlador K(s), se van a introducir unas funciones de ponderaci´on en el esquema de control, de forma que la planta queda aumentada de la siguiente forma: Figura 9: Planta aumentada P(s). La resoluci´on del problema H∞con esta configuraci´on recibe el nombre de problema de sensibilidad mixta S/KS/T. La expresi´on a minimizar en este caso queda de la siguiente manera: kTzw(S)k∞=       WS(s)·So(s) WU(s)·K(s)·So(s) WT(s)·To(s)      ∞ (6) 4.2 ESTIMACI ´ ON DE LA INCERTIDUMBRE MULTIPLICATIVA A LA SALIDA Para la estimaci´on de la incertidumbre del modelo respecto a la del sistema se ha hecho uso de la matriz de funciones de transferencia identificada a esas frecuencias y del modelo linealizado que da lugar a un integrador doble multiplicando a una matriz diagonal. Siguiendo el m´etodo explicado en [7], se ha realizado una normalizaci´on el sistema haciendo uso de las m´aximas se˜nales de control permitidas, Du, y de los m´aximos errores permitidos, De. La normalizaci´on se realiza mediante la expresi´on: GN=D−1 e·G·Du(7) A continuaci´on se estiman las incertidumbres normalizadas respecto al sistema nominal: Gnom(s) =   1 (s−ε)20 01 (s−ε)2  (8) donde el par´ametro εse ha introducido para solventar problemas num´ericos. Se le ha dado un valor de 10−5. As´ı, para cada frecuencia se calcular´a la incertidumbre multiplicativa a la salida del sistema respecto al anterior modelo nominal mediante la expresi´on: Ems(s) = (GN(s)−GNnom(s))·G−1 Nnom(s) (9) donde se ha establecido la hip´otesis de que: G(s)−Gnom(s) = 1 d(s)  n11(s)n12(s) n21(s)n22(s) (10) En la figura 10 se ha representado los valores singulares de Ems(jω) en la frecuencia. 10−5 10−4 10−3 10−2 10−1 100101102 −250 −200 −150 −100 −50 0 50 frecuencia (rad/s) V.S. (dB) σM(Ems) σm(Ems) 1 2 Figura 10: Incertidumbre normalizada. C´alculo de la matriz WT(s) La funci´on de ponderaci´on WT(s) es utilizada para robustificar al sistema frente a incertidumbres multiplicativas. Se ha escogido WT(s) como una matriz de funciones de transferencia cuadrada, diagonal y de dimensi´on la del n´umero de salidas del sistema: WT(s) =   WT1(s) 0 0WT2(s) (11) Las funciones de su diagonal WT i(s) deber´an ser estables, de fase m´ınima y de forma que para cada frecuencia se cumpla que: |WT1(jω)| ≥ σ1(Ems(jω)) (12) |WT2(jω)| ≥ σ2(Ems(jω)) Habr´a que asegurarse de que WT i(s) tenga en alta frecuencia un m´odulo grande para conseguir que T(jω) tenga un m´odulo peque˜no y pueda rechazar ruidos en los sensores. Se han dise˜nado las siguientes WT i(s): WT1(s) = s2 0.0001s2+ 0.02s+ 1 (13) WT2(s) = 0.1·s2 0.0001s2+ 0.02s+ 1 Como puede observarse en la figura 11, la elecci´on de estas funciones de ponderaci´on cumplen las restricciones impuestas en la expresi´on (12). 10−5 10−4 10−3 10−2 10−1 100101102 −250 −200 −150 −100 −50 0 50 100 frecuencia (rad/s) V.S. (dB) σ1 σ2 Wt1 Wt2 Figura 11: Matriz de ponderaci´on WT(s). C´alculo de la matriz WS(s) La matriz WS(s) va a ser utilizada para imponer condiciones de comportamiento sobre el sistema. En concreto, seg´un la ecuaci´on (5), interesa que S(s) tenga valores singulares peque˜nos en baja frecuencia para rechazar perturbaciones a la salida y para reducir el error. Por ello los valores singulares de WS(s) tendr´an que ser grandes en baja frecuencia. Se propone WS(s)como una matriz cuadrada y diagonal de dimensi´on el n´umero de salidas del sistema: WS(s) =   WS1(s) 0 0WS2(s) (14) donde las funciones de transferencia de su diagonal han de cumplir que: |W−1 S1(jw)| ≥ σ1(S(jw)) (15) |W−1 S2(jw)| ≥ σ2(S(jw)) Se han propuesto unas WSi(s)de la forma: WSi(s) = "N √αi·S+ 10(κi−1)·ωT S+N √βi·10(κi−1)·ωT#N (16) siendo ambas iguales (ver figura 12), con los siguientes valores en los par´ametros: •αi= 0.5 : Ganancia en alta frecuencia. •βi= 1e−4: Ganancia en baja frecuencia. •ωT: Frecuencia de corte de WT1(s). •κi= 2 : Determina la frecuencia de corte de la funci´on. •N= 1 : Determina la pendiente de los polos y ceros. 10−5 10−4 10−3 10−2 10−1 100101102 0 10 20 30 40 50 60 70 80 frecuencia (rad/s) dB Ws1 Ws2 Figura 12: Matriz de ponderaci´on WS(s). C´alculo de la matriz WU(s) El objetivo del uso de la funci´on WU(s)ser´a el de disminuir la sobreoscilaci´on de la respuesta temporal afectando poco a la rapidez del mismo. Asimismo, la inclusi´on de WUpermite evitar problemas num´ericos en el c´alculo del controlador. Se ha elegido como matriz WU(s), una matriz identidad de dimensi´on el n´umero de salidas. WU=  WU10 0WU2 =1 0 0 1 (17) 4.3 C ´ ALCULO DEL CONTROLADOR H∞ Siguiendo el m´etodo explicado en [7] se ha calculado el controlador H∞mediante el algoritmo implementado en [1], en el que se obtiene un controlador sub´optimo iterando con un par´ametro γ. El controlador ser´a sintetizado en continuo, discretizado mediante una transformaci´on bilineal de Tustin con T=0.01s y, finalmente, ser´a desescalado mediante la transformaci´on: K(z) = Du·b K(z)·D−1 e(18) donde DeyDuson las mismas matrices que se utilizaron para escalar el sistema seg´un la ecuaci´on (7). El controlador resultante viene dado por las siguientes expresiones: K(z) = 1 dk(z)  nk11(z)nk12(z) nk21(z)nk22(z) (19) nk11(z) = ( 0.0272 z8−0.0872 z5+0.0820 z2+ +0.0183 z7−0.0858 z4+0.0647 z+ −0.0231 z6+0.0042 z3−0.0003 ) nk12(z) = ( 4.4327e−7z8−1.4239e−6z5+1.3444e−6z2+ +2.9551e−7z7−1.4030e−6z4+1.0601e−6z+ −3.7975e−7z6+6.8262e−8z3−4.9787e−9) nk21(z) = (−4.2110e−7z8+8.5055e−7z5+1.0962e−9z2+ −8.5723e−7z7+3.7413e−7z4+5.9606e−8z+ +3.5838e−8z6−5.2926e−8z3+1.0035e−8) nk22(z) = ( 0.03621 z8−0.07218 z5−0.00108 z2+ +0.07273 z7−0.03119 z4−0.00500 z+ −0.00310 z6+0.00445 z3−0.00084 ) dk = ( 1.0000 z8−3.8635 z5+5.5515 z2+ −3.4970 z7+0.9092 z4−0.3206 z+ +0.3463 z6−0.1469 z3+0.0209 ) 5 RESULTADOS DE SIMULACI ´ ON En este apartado se van a presentar unas simulaciones con el controlador dise˜nado utilizando el modelo estimado del sistema real. La figura 13 muestra la respuesta dada por el sistema con incertidumbre ante varios escalones dados en las referencias de gui˜nada (orientaci´on) y cabeceo (elevaci´on). Por otro lado en la figura 14 se representan los valores de la se˜nal de control del sistema con incertidumbre en el acoplamiento. Un hecho importante es que en la planta real se produce saturaci´on 0 20 40 60 80 100 120 −1 0 1 2 3 4 5 6 7Sistema con Incertidumbre Ángulo en grados Tiempo ξ:Cabeceo−Elevación θ: Guiñada−Orientación Figura 13: Respuesta del sistema nominal y del sistema con incertidumbre. de los actuadores cuando esta se˜nal de control supera el valor de 0.5. Puede observarse que se est´a lejos de que se produzca este fen´omeno, por lo que se podr´ıa hacer el control algo m´as agresivo sin problemas de saturaci´on. 0 20 40 60 80 100 120 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 Sistema con Incertidumbre Señal de control Tiempo ξ:Cabeceo−Elevación θ: Guiñada−Orientación Figura 14: Se˜nal de control del sistema con incertidumbre en el acoplamiento. La figura 15 presenta un detalle en el que se observa el efecto de acoplamiento en los cambios bruscos de referencia. De hecho, para el cambio de referencia dado, el sistema se comporta como un sistema de fase no m´ınima. La explicaci´on f´ısica de este fen´omeno es la siguiente: si se proporciona como referencia un incremento de ∆ξhacia arriba (v´ease la figura 16), se traducir´a en un incremento de la velocidad del rotor principal, ∆γ, en sentido contrario a las agujas del reloj. Esto dar´a lugar a un par de reacci´on en sentido contrario en O alrededor de Z0, provo- 1 2 3 4 5 6 7 8 9 10 11 0 0.5 1 1.5 2 2.5 3 3.5 Sistema con Incertidumbre Ángulo en grados Tiempo ξ:Cabeceo−Elevación θ: Guiñada−Orientación Figura 15: Respuesta del sistema no nominal. cando un desplazamiento en orientaci´on ∆θ. Como consecuencia, el sistema aumentar´a la velocidad del rotor de cola ∆αen sentido contrario a las agujas del reloj para controlar la orientaci´on, provocando a su vez un par de reacci´on de sentido contrario en A alrededor de X2. Esto dar´a lugar, inicialmente, a un desplazamiento ∆ξhacia abajo. Figura 16: Sistemas de referencia y nomenclatura. En la figura 17 se muestran las se˜nales de control Al igual que suced´ıa anteriormente, se est´a bastante lejos de los valores que producen saturaci´on en los actuadores (0.6). 6 CONCLUSIONES Se ha realizado un modelo lineal de la planta en torno a un punto de equilibrio identificando din´amica de alta frecuencia, de la cual se ha extra´ıdo un comportamiento de acoplamiento del sistema, y relacionado con ´el, un comportamiento t´ıpico de los sistemas de fase no m´ınima. Se ha dise˜nado un controlador robusto lineal y multivariable para dicha planta. Se ha visto como el sistema identificado, con incertidumbre en el acoplamiento, es controlado y consigue mantener 1 2 3 4 5 6 7 8 9 10 11 0 0.01 0.02 0.03 0.04 0.05 Sistema con Incertidumbre Señal de control Tiempo ξ:Cabeceo−Elevación θ: Guiñada−Orientación Figura 17: Se˜nal de control del sistema con incertidumbre en el acoplamiento. las se˜nales de control dentro del rango de valores que aseguran la no saturaci´on del sistema. Se consigue estabilidad y comportamiento robustos. Por otro lado, se ha llegado a una soluci´on de compromiso entre rapidez y comportamiento. Referencias [1] Balas, G., Doyle, J., Glover, K., Packard, A. and Smith, R. (1998). “µAnalysis Toolbox for use with MatLab”, The MathWorks, Inc. [2] Chiang, Y. and Safonov, R.(1998). “Robust Control Toolbox for use with MatLab”, The MathWorks, Inc. [3] Doyle, J., Grover, K., Khargonekar, P. and Francis, B. (1989). “State-Space Solutions to Standard H2and H∞Control Problems”. IEEE Trans. on Automatic Control. Vol. 34, no. 8, pp. 831-846, August. [4] Feedback Instruments Ltd., (1996). “Twin Rotor Mimo System : Installation and Commissioning. Getting Started. Reference Manual. External Interface to Real Time Kernel. Advanced Teaching Manual I”. Park Road, Crowborough, E. Sussex, TN6 2QR, UK. [5] L´opez-Mart´ınez, M. (1998). “Modelo y control LQR por Par Calculado de un Sistema Mec´anico de 2 GDL Impulsado por Dos H´elices”, Proyecto Final de Carrera. Escuela Superior de Ingenieros. Universidad de Sevilla. [6] Ljung, L. (1995). “System Identification Toolbox for use with Matlab”, The Math- Works, Inc. [7] Ortega, M.G. (2001). “Aportaciones al Control H∞de Sistemas Multivariables”. Tesis Doctoral. Depto. Ingenier´ıa de Sistemas y Autom´atica. Universidad de Sevilla. [8] Rubio, F.R. y L´opez S´anchez, M.J. (1996). “Control Adaptativo y Robusto”. Secretariado de publicaciones de la Universidad de Sevilla. [9] Skogestad, S. and Postlethwaite, I. (1996). “Multivariable Feedback Control. Analysis and Design.” John Wiley & Sons. New York. [10] Spong,M.W. and Vidyasagar,M., (1989). “Robot Dynamics and Control”. John Wiley & Sons. New York. [11] Zames, G. (1966) “On the input-output stability of time-varying nonlinear feedback systems, Part 1.” IEEE Trans. on Automatic Control. Vol. 11.