Control predictivo no lineal basado en modelos de Volterra. Aplicación a una planta piloto

CONTROL PREDICTIVO NO LINEAL BASADO EN MODELOS DE VOLTERRA. APLICACIÓN A UNA PLANTA PILOTO. Jorn K. Gruber 1,Carlos Bordons Dpto. Ingeniería de Sistemas y Automática, Escuela Técnica Superior de Ingenieros, Universidad de Sevilla, Camino Descubrimientos, s/n., 41092 Sevilla Email: [email protected], bor[email protected] Resumen: El presente artículo muestra la aplicación de una estrategia de control predictivo no lineal basado en modelos (Nonlinear Model Predictive Control, NMPC) a una planta piloto en la que se controla la temperatura de un reactor. El controlador está basado en un modelo diagonal de Volterra de segundo orden para considerar efectos no lineales. Para calcular la acción de control se utiliza un método iterativo que requiere pocos recursos computacionales. El comportamiento del proceso y del controlador se presenta mediante resultados experimentales. Finalmente se comparan los resultados experimentales del NMPC con los resultados de un controlador predictivo lineal (MPC ). Copyright c 2007 CEA-IFAC Palabras clave: control predictivo no lineal, modelos de Volterra, secuencia pseudo-aleatoria, identificación, planta piloto. 1. INTRODUCCIÓN Los modelos de Volterra representan la extensión simple y lógica de los modelos de convolución que han sido aplicados con gran éxito en control predictivo lineal basado en modelos (Linear Model Predictive Control, MPC ). Estos modelos presentan genéricamente un buen comportamiento y se puede aprovechar su estructura en el diseño de controladores. En el caso de un modelo de segundo orden, el modelo tiene la forma 1Este trabajo ha sido financiado parcialmente por el Ministerio de Educación y Ciencia a través del proyecto de investigación DPI-2004-07444-C04-01. y(k)=h0+ ∞  i=1 aiu(k−i) + ∞  i=1 ∞  j=i bi,j u(k−i)u(k−j)(1) que corresponde al modelo lineal de convolución con la no linealidad como término adicional y aditivo, siendo y(k)la salida y u(k)la entrada del sistema en el instante actual. Los parámetros h0,aiybi,j representan el offset, los parámetros lineales y los parámetros de segundo orden, respectivamente. El control predictivo no lineal basado en modelos (Nonlinear Model Predictive Control, NMPC )es un campo joven en el que todavía existen pocas referencias de aplicaciones industriales (Camacho y Bordons 2004). Pero la capacidad de considerar http://riai.isa.upv.es ISSN: 1697-7912. Vol. 4, Núm. 3, Julio 2007, pp. 34-45 dinámicas más complejas y no lineales convierte al NMPC en una opción prometedora en el momento de controlar procesos reales. Aunque los procesos industriales muestran en general un comportamiento no lineal, la mayoría de los controladores predictivos están basados en modelos lineales. En los procesos fuertemente no lineales la aplicación 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ámica no lineal. Frente a esta gran ventaja existen también inconvenientes como: la dificultad en la obtención de un modelo no lineal a partir de datos experimentales y el problema del cálculo de la acción de control. Este artículo presenta la aplicación de un control predictivo no lineal a un proceso de laboratorio. El modelo dinámico está basado en un modelo diagonal de Volterra y el cálculo de la acción de control se realiza mediante un método iterativo basado en (Maner et al. 1996, Doyle et al. 2002). El artículo está 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ón 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ón del controlador no lineal al proceso real y los resultados de un controlador predictivo lineal como comparación. Finalmente, en el apartado 6 se presentan las conclusiones derivadas de la aplicación de este tipo de controlador. 2. DESCRIPCIÓN DEL PROBLEMA Para la aplicación de un control no lineal se ha elegido un proceso real representado por una planta piloto. El proceso analizado ha sido estudiado anteriormente por varios autores (Hugo et al. 1988, Szeifert et al. 1995, Cueli 2005). 2.1 Proceso de laboratorio El proceso de laboratorio, representado por una planta piloto, puede simular reacciones químicas exotérmicas a base de cambios en la temperatura. La planta piloto (ver la figura 1) está situada en la Escuela Superior de Ingenieros de la Universidad de Sevilla y ha sido utilizada como planta de referencia para propósitos de control por varios autores (Ramírez et al. 1999). Los elementos principales de la planta piloto son el reactor, el intercambiador de calor, la camisa de refrigeración y la válvula de recirculación para manipular el flujo en la camisa de refrigeración (ver la figura 2). Figura 1. Planta piloto utilizada para aplicar el NMPC. Intercambiador de calor v8 Fj,e Fj,s Ff,e Ff,s Figura 2. Diagrama de la planta piloto con los cuatro elementos: reactor, intercambiador de calor, camisa de refrigeración y válvula de recirculación. La reducción de la temperatura del contenido del reactor se realiza mediante una camisa de refrigeración. La disipación de calor provocado por la camisa de refrigeración puede ser regulada mediante la válvula v8que manipula el flujo en la camisa. La relación experimental entre el flujo Fjque circula por la camisa de refrigeración y la apertura de la válvula v8se puede ver en la figura 3. El líquido de refrigeración, agua, entra en la camisa con una temperatura constante. Con el flujo de alimentación Ff,e se suministra el reactivo de la reacción al reactor con el fin de mantener la reacción activa. Antes de entrar al reactor, el flujo de alimentación pasa por el intercambiador de calor para reducir la diferencia de temperatura entre el flujo de alimentación y el contenido del reactor. Para mantener el volumen del contenido del reactor constante se utiliza el flujo saliente Ff,s. Por consiguiente, como el flujo entrante y el flujo saliente tienen el mismo caudal y casi la misma temperatura, los dos flujos juntos apenas provocan cambios en la energía interna del reactor. J. K. Gruber, C. Bordons 35 Fjl s v8[%] 020 40 60 80 100 0 0.06 0.12 0.18 Figura 3. Relación entre el flujo Fjcirculando por la camisa de refrigeración y la apertura de la válvula v8. El calor de la reacción, que es función de la concentración del reactivo, puede ser emulado empleando para tal efecto la resistencia eléctrica según la ley que rige esta reacción. Este procedimiento de simulación de reacciones ha sido empleado por diversos autores, como (Santos et al. 2001). Para información más detallada sobre el proceso ver las publicaciones (Cueli 2005, Cueli y Bordons 2006). 2.2 Modelo matemático Aunque no es necesario disponer de un modelo matemático para el diseño del controlador basado en modelos de Volterra, se presenta en esta sección el modelo del proceso para poner de manifiesto las características no lineales del mismo. El modelo matemático también justifica la forma de emular el calor generado en la reacción química mediante la resistencia. En el modelo, representando un proceso de refino y anteriormente utilizado en (Lee et al. 2000), se ha considerado que el flujo de alimentación y el flujo saliente del reactor son equivalentes, así que Ff=Ff,e =Ff,s. Como consecuencia de esta suposición el volumen VylamasaMson constantes. Los cambios de la temperatura en el reactor de la planta piloto vienen definidos como: dT dt=−Fj V(Tj,e −Tj,s) +(−ΔH)·V MCp k0e−E/(RT)C2 A(2) donde el primer término corresponde a la disipación de calor por la camisa de refrigeración y el segundo término denota la energía generada por la reacción química exotérmica. Las variables Fj, Tj,e yTj,s representan el flujo que circula por la camisa de recirculación así como la temperatura del agua entrando y saliendo de la camisa, respectivamente. CAes la concentración del reactivo en el reactor. Como el flujo de alimentación pasa por el intercambiador de calor y entra en el reactor casi con la misma temperatura que tiene el contenido del reactor, se ha supuesto que el flujo de alimentación ni suministra ni extrae calor del reactor. Tabla 1. Parámetros del modelo y variables constantes Parámetro Valor Unidad k01.2650 ·1017 l mol·s Cp4.18 kJ K·kg ΔH−105.57 kJ mol E/R 13550 K Variable Valor Unidad V25 l M25 kg CA,e 1.2mol l Ff0.05 l s Tj,e 291.15 K La concentración CAdel reactivo en el reactor de la planta se calcula como: dCA dt=Ff V(CA,e −CA)−k0e−E/(RT)C2 A(3) donde CA,e representa la concentración del reactivo en el flujo de alimentación. Los parámetros del modelo y las variables con valores constantes están mostrados en la tabla 1. Como la temperatura del fluido de refrigeración es variable y no corresponde a la temperatura en el reactor, se ha utilizado la siguiente aproximación, obtenida experimentalmente, del modelo matemático para calcular el intercambio de calor: Fj·(Tj,s −Tj,e)=T−α β(1 −e−γFj)(4) con los parámetros α= 292.19 K,β=14.94 s/ly γ=13.18 s/l. La mayor parte de la no linealidad del modelo matemático es estática, pero debido a los términos exponenciales y cuadráticos en las ecuaciones diferenciales de la temperatura (2) y de la concentración (3) el modelo posee también una dinámica no lineal. La figura 4 muestra la no linealidad estática del modelo (2)-(4) a la que hay que añadir la relación no lineal entre el flujo de refrigeración Fjylaaperturadelaválvulav8, según se mostró en la figura 3. 3. IDENTIFICACIÓN DEL SISTEMA La identificación de parámetros de modelos consta en la mayoría de los casos de cinco partes: elección de la estructura del modelo, elección de señales de entrada adecuadas, ensayos con aplicación de las señales elegidas para conseguir datos de tipo entrada–salida, la propia identificación de los parámetros del modelo y la validación del modelo identificado. 36 Control Predictivo No Lineal Basado en Modelos de Volterra. Aplicación a una Planta Piloto T[K] CAmol l Fjl s Fjl s 0.02 0.02 0.05 0.05 0.08 0.08 0.11 0.11 0.14 0.14 0.17 0.17 0 0.2 0.4 0.6 300 325 350 375 Figura 4. Temperatura Ty concentración CAdel reactor en régimen estacionario en función del flujo de refrigeración Fjcalculados mediante el modelo matemático. 3.1 Estructura del modelo A pesar del hecho de 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 tk. Por estas razones existe la necesidad de describir sistemas continuos mediante modelos en tiempo discreto de la forma: y(k+1)=f(y(k),y(k−1), ..., u(k),u(k−1), ...)(5) Finalmente, para representar la dinámica no lineal del proceso mediante un modelo discreto en la forma de la ecuación (5), se ha elegido un modelo de Volterra de segundo orden finito, truncando las sumas infinitas de la expresión (1). Para sistemas con olvido (Boyd y Chua 1985) establecen que pueden ser aproximados por modelos de Volterra finitos. En el caso de un sistema estable, los parámetros lineales aitienden a 0para i→∞y los de segundo orden bi,j a0para i→∞yj→∞. Entonces se utiliza un orden de truncamiento de N1yN2con ai≈0∀i>N 1ybi,j ≈0∀i> N2,j>N 2, respectivamente. De esta manera el modelo puede ser definido como: y(k)=h0+ N1  i=1 aiu(k−i) + N2  i=1 N2  j=i bi,j u(k−i)u(k−j)(6) Otra posibilidad para identificar la dinámica del proceso sería el uso de un modelo con respuesta infinita al impulso (Infinite Impulse Response, IIR). El uso de este tipo de modelo resultaría en un control más complejo por lo que se ha renunciadoaestetipodemodeloparanoperder sencillez en el momento de controlar el proceso. 3.2 Señal de entrada Para la identificación de los parámetros del modelo diagonal de Volterra hay que utilizar datos de tipo entrada–salida adecuados. La obtención de estos datos requiere la elección de una señal de entrada que excita suficientemente el sistema. Los datos adquiridos deben permitir la posterior identificación de la dinámica compleja del sistema. En muchas identificaciones de sistemas se utiliza como señal de entrada el ruido blanco de Gauss (Nowak y Van Veen 1994). Esta señal es interesante en teoría por contener todas las frecuencias con la misma potencia y poseer una autocorrelación de cero (Parker et al. 2001). Pero la aplicación de la misma puede dañar el proceso por el movimiento constante de la válvula. Debido a este problema hay que encontrar una señal de entrada más adecuada que no provoca daños en el sistema pero permita la identificación de los parámetros del modelo. Para la identificación de sistemas lineales se usa frecuentemente una secuencia binaria pseudoaleatoria (Pseudo Random Binary Sequence, PRBS). Desafortunadamente, las PRBS tienen una gran deficiencia en el momento de identificar sistemas no lineales debido a la excitación insuficiente del sistema (Nowak y Van Veen 1994). Por los problemas con el ruido blanco de Gauss y la PRBS se ha elegido finalmente como señal de entrada una secuencia multinivel pseudo-aleatoria (Pseudo Random Multilevel Sequence,PRMS). Las PRMS tienen, en contraste a las PRBS,no solamente dos niveles, sino varios. La ampliación del número de diferentes niveles en la señal de entrada permite una mejor excitación del sistema. Las PRMS con N+1 niveles pueden excitar suficientemente sistemas de Volterra de orden N (Nowak y Van Veen 1994). Para el orden N=2 del modelo de Volterra que se utiliza en este artículo se usa una PRMS con tres niveles u(k)= {u−,u 0,u +}(ver la figura 5). Para los tres niveles se han elegido los valores u−=40%,u0=60%yu+=80%que resultan en flujos de refrigeración de Fj(u−)=0.03 l/s, Fj(u0)=0.06 l/syFj(u+)=0.12 l/s. Examinando los resultados del modelo matemático presentados en la figura 4 se puede observar para los valores del flujo de refrigeración una fuerte no linealidad en la temperatura en régimen estacionario. J. K. Gruber, C. Bordons 37 u(k) u+ u0 u− k Figura 5. Secuencia (PRMS) de entrada para la identificación de modelos de Volterra hasta segundo orden. T[K] 310 320 330 340 350 v8[%] 40 50 60 70 80 t[h] t[h] 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 Figura 6. Resultados del ensayo de la planta piloto con la entrada calculada mediante el método de PRMS con tres niveles de la apertura de la válvula de recirculación v8={40,60,80}%. 3.3 Ensayo Con el objetivo de conseguir datos para la identificación de parámetros del modelo de Volterra de segundo orden se han realizado varios ensayos de la planta piloto emulando la reacción química. Durante los ensayos se ha aplicado una PRMS en la válvula con períodos variables, pero bastante largos para observar la reacción de la planta piloto a cambios en la entrada (ver la figura 6). 3.4 Identificación de parámetros Finalmente, con los datos entrada–salida obtenidos en los ensayos, se ha realizado la identificación de los parámetros del modelo de Volterra de segundo orden mediante el método de mínimos cuadrados. Antes de identificar los parámetros, se ha normalizado la entrada de la siguiente manera: ˜u(k)=u(k)−u0 Δu(7) con u0=60yΔu=20así que la entrada normalizada tiene los niveles ˜u={−1,0,1}. ak[−] 0.2 0 −0.2 −0.4 −0.6 −0.8 k[−] 020 40 60 80 100 bj,k [−] 0 0.1 0.2 0.3 j[−] 0510 15 20 k[−] 0 5 10 15 20 Figura 7. Parámetros identificados del modelo no lineal con N1= 100 yN2=20. En primer lugar se ha identificado un modelo de Volterra completo con N1= 100,N2=20yun tiempo de muestro de tm=40s, dando lugar a los parámetros que se muestran en la figura 7. Se puede ver en la gráfica que el orden de truncamiento de N2=20no permite suficientemente la identificación de la dinámica no lineal sino corta partedelamisma.Comolamayorpartedela dinámica no lineal se encuentra en la diagonal principal (j=k) se ha decidido de utilizar un modelo diagonal de Volterra de segundo orden. Con los parámetros no diagonales de segundo orden bi,j =0∀i=jel modelo se puede representar de la siguiente manera: y(k)=h0+ N1  i=1 aiu(k−i)+ N2  i=1 biu(k−i)2(8) y permite un aumento de N2sin un incremento excesivo del número de parámetros. Con el modelo diagonal de Volterra de segundo orden se ha realizado una identificación con N1= 100 yN2=40resultando en los parámetros que muestra la figura 8. En este caso se puede observar que los últimos parámetros de segundo orden tienden a cero y, como consecuencia, el orden de truncamiento de N2=40permite una identificación adecuada de la dinámica no lineal. El offset del modelo tiene un valor de h0= 321.07 K. Para la validación del modelo diagonal identificado se ha utilizado datos de un segundo ensayo con una PRMS como señal de entrada, con los mismos niveles de la válvula de v8={40,60,80}%pero con períodos más cortos. La figura 9 muestra en la gráfica superior la identificación y en la gráfica 38 Control Predictivo No Lineal Basado en Modelos de Volterra. Aplicación a una Planta Piloto ak[−] 0.2 0.2 0 0 −0.2 −0.2 −0.4 −0.6 −0.8 bk[−] −0.1 0.1 0.3 k[−] k[−] 0 0 10 20 20 30 40 40 60 80 100 Figura 8. Parámetros identificados del modelo diagonal con N1= 100 yN2=40. T,Tid [K]T,Tvd [K] 310 310 315 320 320 325 330 330 335 340 350 t[h] t[h] 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 Figura 9. Comparación de los datos experimentales (—) y del modelo (---). Arriba: resultados de la identificación, abajo: resultados de la validación. Tabla 2. Comparación de errores cuadráticos medios Error Identificación Validación ECM 0.49972 0.55807 inferior la validación del modelo. La tabla 2 recoge los errores cuadráticos medios (ECM )dela identificación y de la validación del modelo. La figura 10 compara el modelo de Volterra,el modelo matemático de primeros principios y valores experimentales de la planta piloto en régimen estacionario. Para evitar problemas en la posterior implementación del controlador predictivo se normalizan la salida del sistema y la referencia, de tal manera que desaparece el offset del modelo: ˜y(k)=y(k)−h0(9) ˜r(k)=r(k)−h0(10) T[K] v8[%] 20 40 60 80 100 290 320 350 380 Figura 10. Temperatura en régimen estacionario delmodelodeVolterra (—), del modelo de primeros principios (---) y de la planta real (◦)enfuncióndelaaperturadelaválvulav8. 4. LEY DE CONTROL Para predecir el futuro comportamiento del proceso se utiliza el modelo no lineal de Volterra que se puede considerar como una extensión natural de los modelos usados en los controladores predictivos más extendidos en la industria, como el Dynamic Matrix Control (DMC ) (Cutler y Ramaker 1980). Para el modelo de segundo orden los futuros valores de la salida se pueden expresar como (Doyle et al. 2002): ˜y =G˜u+c+f(11) con c=H˜upas +d+g(12) donde las matrices GyHcontienen los parámetros lineales para considerar los términos lineales del futuro y del pasado, respectivamente. El vector ˜u contiene los futuros valores de la acción de control y, análogamente, el vector ˜upas representa los pasados valores de la acción de control. El vector d=[d(k+1),d(k+1), ..., d(k+1)] T contiene la diferencia entre la salida del proceso y la salida del modelo en el instante k. Los términos de segundo orden de tipo futuro–futuro y futuro– pasado se incluyen mediante el vector f, el vector grepresenta los elementos de segundo orden de tipo pasado–pasado. La definición de las matrices GyHasí como los vectores fygse explican en detalle en el Apéndice A. Para considerar los incrementos de la acción de control en la función de coste de la ley de control se transforma la ecuación (11) de tal manera que dependa de los incrementos de la acción de control ˜y =˜y(˜u)−→ ˜y =˜y(Δ˜u). Con la transformación (ver el Apéndice B), el cálculo de los futuros valores de la salida queda como: ˜y =G∗Δ˜u+c∗+f(13) con c∗=H∗˜upas +d+g(14) Con la transformación del sistema no lineal se puede incluir la ponderación de la acción de control en la función de coste que entonces viene definida como: J. K. Gruber, C. Bordons 39 J= P  i=1 (˜y(k+i)−˜r(k+i))2 + M−1  i=0 λΔ˜u(k+i)2(15) donde λrepresenta el parámetro de ponderación de la acción de control. Con la anterior representación del modelo no lineal y la función de coste se ha elegido un enfoque iterativo para calcular la acción de control. Este enfoque, está basado en el trabajo de (Maner et al. 1996, Doyle et al. 2002), al que se le ha añadido la consideración de la acción de control en la función de coste es el siguiente: Paso 1: hacer i=1yΔ˜u(0)(k)=0 Paso 2: resolver las ecuaciones sin restricciones a=(˜r −c∗−f)TG∗T(16) Δ˜u=G∗TG∗+λI−1a(17) Paso 3: comprobar si el primer elemento de la Δ˜ucalculada cumple la tolerancia δen la condición Δ˜u(i)(k)−Δ˜u(i−1)(k)<δ (18) Paso 4: si la anterior condición se cumple, hacer Δ˜u(k)=Δ˜u(i)(k)y la nueva acción de control en el instante k, con ˜usr(k)=˜u(k− 1) + Δ˜u(k)y la saturación física será: ˜u(k)=⎧ ⎨ ⎩ ˜umin si ˜usr(k)<˜umin ˜usr(k)enotrocaso ˜umax si ˜usr(k)>˜umax (19) Si la anterior condición no se cumple, recalcular fusando ˜u =⎡ ⎢ ⎢ ⎣ ˜u(k−1) ˜u(k−1) ... ˜u(k−1) ⎤ ⎥ ⎥ ⎦+ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Δ˜u(i)(k) 1  j=0 Δ˜u(i)(k+j) ... M−1  j=0 Δ˜u(i)(k+j) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (20) hacer i=i+1y volver a paso 2. Como se puede ver en el paso 4, en caso de no cumplir la condición de convergencia, en cada iteración el término cuadrático futuro fse calcula con el vector de acciones de control ˜u calculada en el paso 2 de la misma iteración. Con el nuevo vector fse calcula en la siguiente iteración un nuevo vector de incrementos de la acción de control. Después de cada iteración se calcula la diferencia entre el primer elemento del vector de incrementos de la actual iteración y de la anterior iteración. Se repite este procedimiento hasta que T[K] 320 323 326 329 332 v8[%] 35 45 55 65 75 t[min] t[min] 0 0 20 20 40 40 60 60 80 80 100 100 Figura 11. Cambios en la referencia del proceso simulado mediante el modelo de Volterra. Comparación del controlador propuesto (—) y un controlador basado en SQP (---). la diferencia cumpla la condición de convergencia. Para detalles sobre la convergencia del algoritmo ver el trabajo de (Economou 1986). Como se considera el término fdurante el paso 2 como constante, este término puede seguir dependiendo de los valores absolutos de la acción de control y no requiere una transformación para depender de los incrementos de la acción de control. En el trabajo de (Dorado 2006) se presenta cómo una adecuada elección de la ponderación de la acción de control en la función de coste puede garantizar la convergencia hacia el óptimo global. La aplicación posterior de la saturación se puede evitar integrando las restricciones en el paso 2 y resolver el problema de optimización mediante programación cuadrática (Quadratic Programming,QP) (Doyle et al. 2002). La nueva acción de control u(k)se calcula mediante denormalización de la acción de control ˜u(k)y se deduce directamente de la ecuación (7): u(k)=˜u(k)·Δu+u0(21) 5. RESULTADOS EXPERIMENTALES En primer lugar se han ajustado los parámetros del controlador propuesto simulando la reacción química mediante el modelo de Volterra.Pormedio de las simulaciones se ha fijado la ponderación de la acción de control a λ=0.8, los ordenes de truncamiento a N1= 100 yN2=20, el horizonte de predicción a P= 100 y el horizonte de control aM=15. Además se ha utilizado un controlador basado en el algoritmo SQP (Sequential Nonlinear Programming) de Programación No Lineal con el fin de verificar los resultados del controlador propuesto. La figura 11 muestra una comparación de los resultados casi idénticos de los dos controladores en el caso de cambios en la referencia. El proceso en cuestión está controlado mediante el sistema de control distribuido Simatic-IT que 40 Control Predictivo No Lineal Basado en Modelos de Volterra. Aplicación a una Planta Piloto u ˜u rr∗˜ry ˜y T ˜r=r∗−h0 1 Kts+1 ˜y=y−h0 u=˜uΔu+u0Planta piloto NMPC Filtro de primer orden Normalización Normalización Algoritmo iterativo Denormalización Sistema Figura 12. Diagrama de bloques de la estructura de control con el NMPC para controlar la temperatura de la planta piloto. incluye un Automata Programable denominado PMC 10. El algoritmo de control se ejecuta en Matlab e intercambia información con el sistema de control distribuido a través de OPC (OLE for Process Control). El algoritmo de control predictivo no lineal ha sido aplicado a la planta piloto usando la estructura de control que muestra la figura 12. En un primer paso se filtra la referencia ry se normalizan la referencia filtrada r∗y la salida del proceso y. Con los valores normalizados el NMPC calcula la nueva acción de control ˜umediante el algoritmo iterativo. Finalmente, se aplica la acción de control denormalizada ua la planta piloto. El uso del filtro de la referencia permite la elección de un valor bajo para el parámetro λde la ley de control. De esta manera, la combinación de filtro de referencia y valor bajo de λresultanenuncomportamiento suave del proceso en el momento de cambiar el punto de operación y una reacción más agresiva a divergencias en régimen estacionario. Con el controlador predictivo no lineal se han realizado dos series de ensayos, una con cambios en la referencia y otra con perturbaciones. Para poder comparar los datos obtenidos, se ha aplicado también un controlador MPC lineal mediante la misma estructura de control (ver la figura 12) al proceso. Para los ensayos con el NMPC se han utilizado los parámetros N1= 100 yN2=40 (orden de truncamiento), P= 100 (horizonte de predicción), M=15(horizonte de control) y λ=0.8(ponderación). El filtro de la referencia tiene una constante de tiempo de Kt= 100 s. Los mismos parámetros se han utilizado para el controlador predictivo lineal. Primero se ha realizado con el proceso controlado por el NMPC un ensayo con varios cambios en la referencia (ver la figura 13). En la gráfica se pueden ver desde arriba a abajo: la temperatura y la referencia, la concentración del reactivo, la apertura de la válvula y las iteraciones necesarias para cumplir la condición 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 diT[K] 320 323 326 329 332 CA[mol/l] 0.1 0.15 0.2 v8[%] 30 40 50 60 70 80 Iteraciones [−] 0 2 4 6 t[min] t[min] t[min] t[min] 0 0 0 0 20 20 20 20 40 40 40 40 60 60 60 60 80 80 80 80 100 100 100 100 Figura 13. Resultado del ensayo de la planta piloto con cambios en la referencia controlado por el NMPC. De arriba a abajo: temperatura y referencia, concentración del reactivo, acción de control, iteraciones para convergencia. vergencias entre temperatura y referencia, pero el controlador no lineal compensa estas divergencias en un tiempo razonable. Las acciones de control muestran una reacción rápida a las divergencias entre la temperatura y la referencia y presentan un comportamiento suave que no ha sido limitado por la saturación. Con respecto a las iteraciones para cumplir la condición de convergencia se puede ver que el algoritmo necesitaba entre 2 y 6 iteraciones, que significa un cálculo rápido con los ordenadores actuales. La figura 14 muestra los resultados del ensayo con cambios en la referencia controlado por el MPC lineal. Como se puede ver, la temperatura J. K. Gruber, C. Bordons 41 T[K] 320 323 326 329 332 CA[mol/l] 0.1 0.15 0.2 v8[%] 30 40 50 60 70 80 t[min] t[min] t[min] 0 0 0 20 20 20 40 40 40 60 60 60 80 80 80 100 100 100 Figura 14. Resultado del ensayo de la planta piloto con cambios en la referencia controlado por el MPC. De arriba a abajo: temperatura y referencia, concentración del reactivo, acción de control. muestra oscilaciones después de los cambios en la referencia. En régimen estacionario, temperatura y concentración del reactivo presentan más oscilaciones que el proceso controlado por el NMPC. Con respecto a la acción de control, la apertura de la válvula presenta oscilaciones considerables en régimen estacionario. En segundo lugar se ha comprobado el comportamiento del proceso controlado por el NMPC ante perturbaciones. Como perturbaciones se han elegido el flujo de alimentación Ff, ver la ecuación (3), ylaaperturadelaválvulav8. La figura 15 muestra los resultados de un ensayo con las perturbaciones mencionadas. Al proceso en régimen estacionario se ha aplicado en t=20min una perturbación aditiva en el flujo de alimentación de ΔFf=0.005 l/s que corresponde a un 10 % de la gama. Después de 60 minutos, en t=80min, se ha quitado la perturbación para aplicar en t= 140 min una perturbación en la válvula de Δv8=5%.Comose puede ver en los resultados, el controlador NMPC consigue compensar los errores provocados por las perturbaciones, pero tarda en todos los casos entre 30 y 40 minutos para quitar la divergencia. Después de la compensación de los errores, el proceso muestra oscilaciones significativas en la temperatura, la concentración y la acción de control. El número de iteraciones necesarias para calcular la acción de control no se ha visto influido por las perturbaciones. Para poder comparar los resultados del NMPC, se ha realizado un ensayo de la planta piloto con las mismas perturbaciones y el controlador preT[K] 320 321 322 CA[mol/l] 0.18 0.19 0.2 v8[%] 55 60 65 70 Iteraciones [−] 0 1 2 3 4 t[min] t[min] t[min] t[min] 0 0 0 0 50 50 50 50 100 100 100 100 150 150 150 150 200 200 200 200 Figura 15. Resultado del ensayo de la planta piloto con perturbaciones controlado por el NMPC. De arriba a abajo: temperatura y referencia, concentración del reactivo, acción de control, iteraciones para convergencia. dictivo lineal (ver la figura 16). El MPC también logra compensar los errores provocados por las perturbaciones, pero la temperatura muestra oscilaciones muy grandes alrededor de la referencia. La concentración del reactivo y la apertura de la válvula oscilan fuertemente durante todo el ensayo de tal manera que el MPC parece incapaz de reducir las oscilaciones presentes. Las divergencias entre temperatura y referencia durante los dos ensayos con perturbaciones (NMPC yMPC ) son similares, pero las oscilaciones fuertes en la temperatura, la concentración y la apertura de la válvula muestran que el MPC obtiene resultados inferiores en comparación con el NMPC. Para cuantificar la diferencia se ha utilizando la suma de errores cuadráticos = k=1 (r(k)−y(k))2(22) La tabla 3 recoge los resultados de los ensayos realizados. Puede observarse cómo en ambos casos se consigue una mejora, sobre todo teniendo en cuenta que la solución se ha implementado de forma sencilla sin necesidad de haber complicado la estrategia de control. En la tabla 4 se muestran los tiempos totales necesarios para calcular las acciones de control en 42 Control Predictivo No Lineal Basado en Modelos de Volterra. Aplicación a una Planta Piloto