Aplicación de MPC no lineal no lineal basado en modelos de Volterra a una planta piloto

APLICACI ´ ON DE MPC NO LINEAL BASADO EN MODELOS DE VOLTERRA A UNA PLANTA PILOTO J. K. Gruber y C. Bordons Dpto. Ingenier´ıa de Sistemas y Autom´atica, Escuela T´ecnica Superior de Ingenieros Universidad de Sevilla, Camino Descubrimientos, s/n., 41092 Sevilla [email protected], [email protected] Resumen El presente documento muestra la aplicaci´on de una estrategia de control predictivo no lineal basado en modelo (NMPC) a un proceso de laboratorio: una planta piloto con una v´alvula de recirculaci´on manipulable para controlar la temperatura en el reactor. El controlador est´a basado en un modelo diagonal de Volterra de segundo orden para considerar efectos no lineales. Para calcular la acci´on de control se utiliza un m´etodo iterativo que requiere pocos recursos computacionales. El comportamiento del proceso y del controlador se presentan mediante resultados experimentales Palabras clave: control predictivo no lineal, modelos de Volterra, secuencia pseudo-aleatoria. 1 INTRODUCCI ´ ON Los modelos de Volterra representan la extension simple y l´ogica de los modelos de convoluci´on que han sido aplicados con gran ´exito en control predictivo lineal basado en modelos (linear model predictive control, MPC). Estos modelos presentan gen´ericamente un buen comportamiento y se puede aprovechar su estructura en el dise˜no de controladores, especialmente en el caso de modelos de segundo orden. En el caso de un modelo de segundo orden y truncamiento de t´erminos, el modelo tiene la forma y(k) = h0+ N1 X i=1 a(i)u(k−i) + N2 X i=1 N2 X i=j b(i, j)u(k−i)u(k−j) (1) que corresponde al modelo lineal de convoluci´on con la no linealidad como t´ermino adicional y aditivo. El control predictivo no lineal basado en modelos (nonlinear model predictive control, NMPC) es un campo joven y todav´ıa existen pocas referencias de aplicaciones industriales [1]. Pero la capacidad de considerar din´amicas m´as complejas y no lineales convierte el NMPC en una opci´on prometedora en el momento de controlar procesos reales. Aunque los procesos industriales demuestran en general un comportamiento no lineal, la mayor´ıa de los controladores predictivos est´an basados en modelos lineales. Sin embargo, los procesos fuertemente no lineales justifican el uso de NMPC. En los procesos fuertemente no lineales la aplicaci´on de MPC lineal puede ser poco efectiva y requiere el empleo de controladores no lineales. La principal ventaja de NMPC representa la posibilidad de considerar la din´amica no lineal. Enfrente a esta gran ventaja existen tambi´en desventajas que son: la obtenci´on de un modelo no lineal a partir de datos experimentales y el problema de optimizaci´on del c´alculo de la acci´on de control. Este documento presenta la aplicaci´on de control predictivo no lineal a un proceso de laboratorio. El modelo din´amico est´a basado en un modelo de Volterra y la optimizaci´on del c´alculo de la acci´on de control se realiza mediante un m´etodo iterativo presentado por Doyle et al. [5]. El documento est´a distribuido de la siguiente manera: en el apartado 2 se explica el proceso utilizado y el problema que resulta del mismo, en el apartado 3 se presenta la identificaci´on del modelo de Volterra y en el 4 se muestra la ley de control predictivo no lineal. En el apartado 5 se presentan los resultados de la aplicaci´on del controlador no lineal al proceso real y en el 6 se presentan las conclusiones derivadas de la aplicaci´on de este tipo de controlador. 2 DESCRIPCI ´ ON DEL PROBLEMA Para la aplicaci´on de un control no lineal se ha elegido un proceso real representado por una planta piloto, estudiado anteriormente por varios autores [6, 10, 2]. La planta piloto puede simular reacciones qu´ımicas exot´ermicas a base de cambios en la temperatura. Los cuatro elementos principales de la planta son: el reactor, la resistencia, el intercambiador de calor y la v´alvula de recirculaci´on, ver la figura 1. XXVII Jornadas de Automática Almería 2006 - ISBN: 84-689-9417-0 1235 Intercambiador de calor V´alvula Resistencia Reactor Figura 1: Diagrama de la planta piloto con los cuatro elementos: reactor, v´alvula de recirculaci´on, resistencia y intercambiador de calor. Mediante la resistencia se puede suministrar energ´ıa cal´orica y en consecuencia aumentar la temperatura del contenido del reactor. Para la reducci´on de la temperatura del contenido del reactor se utiliza el intercambiador de calor. El transporte de energ´ıa cal´orica desde el reactor al intercambiador de calor se puede regular mediante la v´alvula. El contenido saliente del reactor que pasa por el intercambiador de calor entra otra vez en el reactor tal que la planta piloto representa un sistema cerrado con respecto a la masa. La planta piloto representada (ver la figura 1) se utiliza en general para simular experimentos de reacciones qu´ımicas exot´ermicas. Para informaci´on m´as detallada sobre el proceso utilizado ver [2, 3]. En este documento no se utiliza un modelo matem´atico para calcular la energ´ıa cal´orica generada por una presunta reacci´on qu´ımica, sino que se fija la resistencia a un valor constante. Con la resistencia fijada a un valor constante, eso significa que la resistencia no representa una entrada manipulable del proceso, la planta piloto representa un sistema SISO (single in single out) con una entrada (la apertura de la v´alvula) y una salida (la temperatura en el reactor). La idea de este articulo es la obtenci´on de un modelo matem´atico discreto para representar el cambio de la temperatura Ten dependencia de la apertura de la v´alvula de recirculaci´on v8con un suministro constante de energ´ıa (resistencia fijado a un valor constante) y su posterior uso para el desarrollo de un controlador predictivo no lineal. Tal modelo se puede escribir de la siguiente manera: y(k+ 1) = f (u(k), u(k−1), u(k−2), . . .) (2) con la salida y(temperatura T) y la entrada u (apertura de la v´alvula v8). T[˚C] 50 60 70 80 v8[%] 40 50 60 70 80 90 100 Figura 2: Relaci´on entre apertura de la v´alvula y temperatura en regimen estacionario y resistencia con un valor constante del 50 %. El modelo para identificar deber´ıa considerar la din´amica y los efectos no lineales del proceso de laboratorio. Una de las no linealidades es la relaci´on entre la apertura de la v´alvula v8y la temperatura Ten regimen estacionario y la resistencia con un valor constante, ver la figura 2. 3 IDENTIFICACI ´ ON DE MODELOS La identificaci´on de par´ametros de modelos consiste en la mayor´ıa de los casos de cuatro partes: elecci´on de la estructura del modelo, elecci´on de se˜nales de entrada adecuadas, ensayos con aplicaci´on de las se˜nales elegidas para conseguir datos de tipo entrada–salida, as´ı como la propia identificaci´on de los par´ametros del modelo. 3.1 ESTRUCTURA DE MODELO A pesar del hecho que los sistemas f´ısicos evolucionan en general en tiempo continuo, la mayor´ıa de esquemas de control basados en modelos utilizan ordenadores para recoger las medidas de los sensores en el instante de tiempo discreto tky realizan las acciones de control algunos instante m´as tarde en tk+n. Por estas razones existe la necesidad de describir sistemas continuos mediante modelos en tiempo discreto. Los modelos de Volterra permiten la representaci´on de sistemas complejos y no lineales mediante modelos discretos y por lo tanto su uso en el control de estos sistemas. Finalmente, se ha elegido un modelo de Volterra de segundo orden que representa un modelo no lineal, con el fin de considerar la din´amica compleja del proceso real. La estructura general de un modelo de Volterra de segundo orden es la siguiente: y(k) = h0+ ∞ X i=1 a(i)u(k−i) + ∞ X i=1 ∞ X j=1 b(i, j)u(k−i)u(k−j) (3) XXVII Jornadas de Automática 1236 Almería 2006 - ISBN: 84-689-9417-0 Con el fin de reducir considerablemente el n´umero de par´ametros para identificar, se ha elegido un modelo diagonal de Volterra que representa un caso especial del modelo de Volterra. En este tipo de modelo se consideran solamente los t´erminos cuadr´aticos de los pasados valores de la entrada y no se utilizan los t´erminos cruzados de los mismos. Esta reducci´on de par´ametros resulta en la siguiente estructura de modelo: y(k) = h0+ ∞ X i=1 a(i)u(k−i) + ∞ X i=1 b(i)u(k−i)2(4) Con el truncamiento a N1t´erminos lineales y N2t´erminos no lineales el modelo diagonal de Volterra tiene la siguiente forma: y(k) = h0+ N1 X i=1 a(i)u(k−i) + N2 X i=1 b(i)u(k−i)2(5) 3.2 SE ˜ NAL DE ENTRADA Para la identificaci´on de los par´ametros del modelo diagonal de Volterra hay que utilizar datos de tipo entrada–salida adecuados. La obtenci´on de estos datos requiere la elecci´on de una se˜nal de entrada adecuada. La se˜nal de entrada debe permitir la identificaci´on de la din´amica compleja del sistema. En muchas identificaciones de sistemas se utiliza como se˜nal de entrada el ruido blanco de Gauss [8]. No obstante, el ruido blanco de Gauss es bastante atractivo en teor´ıa, pero poco atractivo en el momento de aplicarlo al sistema [9]. El uso de ruido blanco como se˜nal de entrada para una v´alvula resulta en un movimiento constante de la misma y puede da˜narla. Por estas razones hay que encontrar una se˜nal de entrada aplicable que no da˜ne la v´alvula pero permita la identificaci´on de los par´ametros del modelo. Para la identificaci´on de sistemas lineales se usa en muchos casos una PRBS (Pseudo Random Binary Sequence). Desafortunadamente, las PRBS tienen una gran deficiencia en el momento de identificar sistemas no lineales debido a la excitaci´on insuficiente del sistema [8]. Por los problemas con el ruido blanco de Gauss y la PRBS se ha elegido finalmente como se˜nal de entrada una PRMS (Pseudo Random Multilevel Sequence). Las PRMS tienen, al contrario que las PRBS, no solamente dos niveles, sino varios u(k) u+ u0 u− k Figura 3: Secuencia (PRMS) de entrada para la identificaci´on de modelos de Volterra hasta segundo orden. niveles. La ampliaci´on del n´umero de diferentes niveles en la se˜nal de entrada permite una mejor excitaci´on del sistema. Las PRMS con N+ 1 niveles pueden excitar suficientemente sistemas de Volterra de orden N[8]. Por el orden N= 2 del modelo diagonal de Volterra que se utiliza en este documento se usa una PRMS con tres niveles u(k) = {u−, u0, u+}, ver la figura 3. 3.3 ENSAYO En un primer paso se han realizado varios ensayos de la planta piloto con el objetivo de conseguir datos para la identificaci´on de un modelo de tipo Volterra. Para los ensayos se ha fijado la resistencia a un valor constante del 50 % (→7500 W) y la apertura de la v´alvula se ha calculado mediante el m´etodo PRMS con tres niveles v8={30,60,90}%. Los ensayos se realizaron con periodos variables de la PRMS, pero bastante largos para observar la reacci´on de la planta piloto a cambios en la entrada, ver la figura 4. T[˚C] 55 65 75 v8[%] 30 40 50 50 60 60 70 70 80 80 90 t[h] t[h] 0 0 1 1 2 2 3 3 4 4 5 5 Figura 4: Resultados de ensayo de la planta piloto con la entrada calculada mediante el m´etodo de PRMS con tres niveles de la apertura de la v´alvula de recirculaci´on v8={30,60,90}%. XXVII Jornadas de Automática Almería 2006 - ISBN: 84-689-9417-0 1237 a(k) [−] −0.2 −0.4 −0.6 −0.8 1 b(k) [−] −0.1 0 0 0.1 0.2 0.3 k[−] k[−] 0 0 10 10 20 20 30 30 40 40 50 50 60 60 70 70 80 80 Figura 5: Par´ametros identificados del modelo no lineal con N1= 80 y N2= 45. 3.4 IDENTIFICACI ´ ON DE PAR ´ AMETROS Finalmente, con los datos entrada–salida obtenidos en los ensayos, se ha realizado la identificaci´on de los par´ametros del modelo diagonal de Volterra de segundo orden. Para la identificaci´on de los par´ametros se ha utilizado el m´etodo de m´ınimos cuadrados. Antes de identificar los par´ametros, se ha normalizado la entrada de la siguiente manera: ˜u=u−u0 ∆u(6) con u0= 60 y ∆u= 30 as´ı que la entrada normalizada tiene los niveles ˜u={−1,0,1}. Durante el proceso de identificaci´on se ha observado que la elecci´on de un tiempo de muestreo de tm= 30 s representa un buen compromiso entre la bondad de la identificaci´on y el n´umero de par´ametros. Adem´as se ha observado que, con el tiempo de muestreo de tm= 30 s, el sistema tiene un retardo de 1 periodo de muestreo. Finalmente, se han identificado N1= 80 par´ametros de primer orden y N2= 45 par´ametros de segundo orden. Los par´ametros identificados se pueden ver en la figura 5. El offset del modelo tiene un valor de h0= 59.1oC. Para evitar problemas en la posterior implementaci´on del controlador predictivo se normaliza la salida del sistema y la referencia, de tal manera que desaparece el offset del modelo: ˜y=y−h0(7) ˜r=r−h0(8) 4 LEY DE CONTROL Para predecir el futuro comportamiento del proceso se utiliza el modelo no lineal identificado en una representaci´on de DMC (Dynamic Matrix Control) [4] extendida al caso no lineal. Para el modelo de segundo orden los futuros valores de la salida vienen definidos como: ˜y =G˜u +c+f(9) con f=˜uTG2˜u (10) c=H˜upas +˜uT pasH2˜upas +d(11) y las matrices G(12) y H(13) que representan la parte lineal del modelo identificado. Las matrices G2yH2consideran la parte no lineal identificada y tienen la forma que se puede ver en (14) y (15), respectivamente. El vector d= [d(k+ 1), d(k+ 1), . . . , d(k+ 1),]Tcontiene la diferencia entre la salida del proceso y la salida del modelo en el instante k. G=               a10... 0 a2a1 ...0 . . .. . ....a1 . . .. . ....a1+a2 . . .. . ..... . . aPaP−1... P+M+1 P i=1 ai               (12) H=             a2a3... ... aN0 a3a4. . . aN−10 0 . . .. . .. . .. . .. . .a1 aP−1aP . . .. . .. . .a1+a2 aP0. . .. . .. . .. . . 0 0 . . . 0 0 0             (13) G2=               b10... 0 b2b1 ...0 . . .. . ....b1 . . .. . ....b1+b2 . . .. . ..... . . bPbP−1... P+M+1 P i=1 bi               (14) XXVII Jornadas de Automática 1238 Almería 2006 - ISBN: 84-689-9417-0 H2=             b2b3... ... bN0 b3b4. . . bN−10 0 . . .. . .. . .. . .. . .b1 bP−1bP . . .. . .. . .b1+b2 bP0. . .. . .. . .. . . 0 0 . . . 0 0 0             (15) Con la anterior representaci´on del modelo no lineal se ha elegido un enfoque iterativo para calcular la acci´on de control. Este enfoque, presentado por Doyle et al. [5, 7], no considera restricciones en la acci´on de control. El esquema para calcular la acci´on de control es el siguiente: •Paso 1: poner i= 1 •Paso 2: resolver las ecuaciones a=(r−c−f)TGT(16) ˜u =GTG−1a(17) •Paso 3: comprobar si la ucalculada cumple la condici´on ˜u(i)(k)−˜u(i)(k−1)< δ (18) •Paso 4: Si la anterior condici´on est´a cumplida, poner ˜u(k) = ˜u(i)(k) y aplica ˜u(k). Si la anterior condici´on no est´a cumplida, recalcular fusando ˜u(i)(k), poner i=i+ 1 y volver a paso 2. Como no se utilizan restricciones para calcular las acciones de control se usa la saturaci´on ˜u(k) =    ˜umin si ˜u(k)<˜umin ˜u(k) en otro caso ˜umax si ˜u(k)>˜umax (19) que limita la acci´on de control al intervalo ˜umin ≤ ˜u(k)≤˜umax. 5 RESULTADOS EXPERIMENTALES El proceso se controla mediante el sistema de control distribuido Simatic-IT que incluye un Automata Programable denominado PMC 10. El algoritmo de control se ejecuta en Matlab e intercambia informaci´on con el sistema de control distribuido a trav´es de OPC. La figura 6 muestra los resultados de un ensayo de la planta piloto controlado mediante el NMPC. Por los 80 elementos lineales del modelo no lineal identificado se ha elegido un orden de truncamiento de N= 80. Para el horizonte de predicci´on y el horizonte de control se han utilizado P= 80 y M= 15, respectivamente. En la gr´afica se pueden ver arriba la temperatura y la referencia, en el centro la acci´on de control y abajo las iteraciones necesarias para cumplir la condici´on de convergencia. La temperatura del proceso sigue bien la referencia y muestra pocas oscilaciones. Solamente en el momento de cambio de la referencia hay grandes divergencias entre temperatura y referencia, pero el controlador no lineal compensa estas divergencias en un tiempo razonable. Las acciones de control muestran considerables oscilaciones a lo largo del ensayo, incluso en instantes en los que no hay divergencias entre temperatura y referencia. En los momentos de cambio de referencia, la acci´on de control (la apertura de la v´alvula) ha sido limitada por la saturaci´on a los valores m´ınimos y m´aximos permitidos. Con respecto a las iteraciones para cumplir la condici´on de divergencia se puede ver que el algoritmo necesitaba entre 3 y 11 iteraciones, que significa un c´alculo r´apido con los ordenadores actuales. Adem´as se han realizado ensayos de control predictivo lineal para obtener datos de comparaci´on. Con el NMPC se ha podido reducir la suma del error cuadr´atico entre salida del proceso y referencia un 5.1 % frente al error del proceso controlado por un MPC. 6 CONCLUSIONES En este documento se ha mostrado la aplicaci´on de un NMPC basado en un modelo diagonal de Volterra a un proceso de laboratorio. El procedimiento de desarrollo del controlador es simple (por el uso de un modelo de Volterra) y r´apido. El c´alculo de la acci´on de control se ha realizado con un m´etodo iterativo que converge en pocas iteraciones y por consiguiente requiere pocos recursos. En los resultados se ha podido ver un buen comportamiento del proceso, pero tambi´en cambios permanentes en la entrada del sistema. Para reducir estos cambios permanentes en la apertura de la v´alvula se puede pensar en una extensi´on de la ley de control, de tal manera que se incluye una ponderaci´on de los cambios de control ∆uen el algoritmo de control. Para reducir el error entre temperatura y referencia se puede incluir restricciones para la acci´on de control, de tal manera que se puede eliminar la saturaci´on utilizada en la entrada del sistema. Esta medida obligar´ıa a cambiar del algoritmo iterativo a proXXVII Jornadas de Automática Almería 2006 - ISBN: 84-689-9417-0 1239 T[˚C] 55 60 65 v8[%] 50 100 Iteraciones [−] 0 0 5 10 t[min] t[min] t[min] 0 0 0 10 10 10 20 20 20 30 30 30 40 40 40 50 50 50 60 60 60 Figura 6: Resultado del ensayo de la planta piloto controlado por el MPC no lineal con c´alculo iterativo de la acci´on de control. Arriba: la temperatura y la referencia, centro: la acci´on de control, abajo: las iteraciones necesarias para cumplir la condici´on de convergencia. gramaci´on cuadr´atica (QP). Agradecimientos Este trabajo ha sido financiado parcialmente por el Ministerio de Educaci´on y Ciencia a trav´es del proyecto de investigaci´on DPI-2004-07444-C04-01. Referencias [1] Camacho, E.F., Bordons, C., (2004) “Model Predictive Control”, Second Edition, Springer, London. [2] Cueli, J.R., (2005) “Control predictivo para procesos repetitivos”, Tesis doctoral, Universidad de Sevilla. [3] Cueli, J.R., Bordons, C. (2005) “Iterative Nonlinear Control of a Semibatch Reactor. Stability Analysis”, In Proceedings of the 44th IEEE Conference on Decision and Control, and the European Control Conference 2005, Sevilla. [4] Cutler, C.R., Ramaker, B.C., (1980) “Dynamic Matrix Control – A Computer Control Algorithm”, In Automatic Control Conference, San Francisco. [5] Doyle, F.J., Pearson, R.K., Ogunnaike, B.A., (2002) “Identification and Control Using Volterra Models”, Springer, London. [6] Hugo, P., Steinbach, J., Stoessel, F., (1988) “Calculation of the maximum temperature in stirred tank reactors in case of breakdown of cooling”, Chemical Engineering Science, Vol. 43, No. 8, pp 2147-2152. [7] Maner, B.R., Doyle, F.J., Ogunnaike, B.A., Pearson, R.K., (1996) “Nonlinear Model Predictive Control of a Simulated Multivariable Polymerization Reactor Using Second-order Volterra Models”, Automatica, pp 1285-1301. [8] Nowak, R.D., Van Veen, B.D., (1994) “Random and Pseudorandom Inputs for Volterra Filter Identification”, IEEE Transactions on Signal Processing, Vol. 42, No. 8, pp 21242135. [9] Parker, R.S., Heemstra, D., Doyle, F.J., Pearson, R.K., Ogunnaike, B.A., (2001) “The identification of nonlinear models for process control using tailored “plant-friendly” input sequences”, Journal of Process Control, Vol. 11, No. 2, pp 237-250. [10] Szeifert, F., Chovan, T., Nagy, L., (1995) “Process dynamics and temperature control of fed-batch reactors”, Computers & and Chemical Engineering, Vol. 19, No. 1, pp 447452. XXVII Jornadas de Automática 1240 Almería 2006 - ISBN: 84-689-9417-0