SciELO - Scientific Electronic Library Online

 
vol.19 número6Ajuste del Equilibrio Químico del Agua Potable con Tendencia Corrosiva por Dióxido de CarbonoCombustión sin Llama de Mezclas Pobres Metano-Aire sobre Óxido de Magnesio Adicionado con Óxido de Calcio índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Información tecnológica

versión On-line ISSN 0718-0764

Inf. tecnol. v.19 n.6 La Serena  2008

http://dx.doi.org/10.4067/S0718-07642008000600011 

 

Información Tecnológica-Vol. 19 N°6-2008, pág.: 103-110
doi:10.1612/inf.tecnol.3971it.07

ARTÍCULOS VARIOS

Estimación Óptima de Parámetros y Determinación del Estado Inicial en Modelos Matemáticos No Lineales: Aplicación a un Reactor Químico

Optimal Parameter Estimation and Determination of the Initial Conditions in Non Linear Mathematical Models:  Application to a Chemical Reactor  

Francisco E. Aban, Raquel L. Michel y Orlando J. Domínguez
Universidad Nacional de Salta, Facultad de Ingeniería, Avda. Bolivia Nº 5150, (4400) Salta - Argentina (e-mail: rmichel@unsa.edu.ar, orlando@unsa.edu.ar)


Resumen

Se presenta un nuevo algoritmo de cálculo para la determinación de las condiciones iniciales y el ajuste de parámetros en modelos matemáticos no lineales. El trabajo fue desarrollado sobre la base de modelos derivados de las ecuaciones de balance de propiedades. Se estableció una relación matricial entre las variaciones en los grados de libertad del proceso y las variaciones en las variables de entrada para formular la función objetivo, la cual fue optimizada empleando métodos de proyección ortogonal. La técnica elaborada fue empleada en un reactor químico del tipo agitado continuo obteniéndose resultados consistentes con su comportamiento dinámico del reactor.

Palabras clave: modelos no lineales, determinación de parámetros, proyección ortogonal, dynamic reactor


Abstract

A novel algorithm is presented for the optimal estimation of parameters and for determining the initial conditions in non linear mathematical models. The work was developed based on models derived from property balance equations. A matrix relation between the variation of degrees of freedom of the process and the variations of the input variables was used to formulate the objective function, which was optimized via orthogonal projection methods. The technique was applied to a continuous stirred tank reactor obtaining results that were consistent with the dynamic behavior of the reactor.

Keywords: non linear models, parameters determination, orthogonal projection, dynamic reactor


INTRODUCCIÓN

Tanto en los problemas globales de optimización como en la síntesis de estructuras de control, la determinación óptima de las variables de estado es una cuestión que demanda una especial atención por la incidencia que tiene en la validez de los resultados a obtener.

A la fecha existe una profusa producción, tanto para sistemas determinísticos como para estocásticos, donde se pone de manifiesto de una manera muy clara, la importancia de una adecuada selección de los procedimientos de cálculo involucrados en la resolución del problema global de la determinación del estado. Así se encuentran aplicaciones del método de colocación ortogonal en combinación con diversas estrategias de optimización, al respecto se puede citar, la búsqueda de parámetros de transferencia de calor en lechos empacados (Tsang et al.,1976); la evaluación del estado y los parámetros en un reactor de lecho fijo combinado con programación no lineal empleando una estrategia secuencial de optimización (Cheng y Yuan, 1996); el ajuste de las suposiciones iniciales de parámetros cinéticos y de transferencia de calor en reactores de lecho fijo, desacoplando las ecuaciones de balance del modelo, encontrándose pocas ventajas en esta estrategia (Cheng et al., 1996); el control auto regulado de columnas de destilación (Karacan et al., 1997); en combinación con la técnica de optimización de evolución diferencial en la estimación efectiva de parámetros de transferencia de calor en reactores de lecho empacado (Babu y Sastry, 1999); para la evaluación de parámetros y estado en procesos biológicos (Damak, 2007). Recientemente se pudo comprobar que la solución obtenida mediante está vía puede no ser única (Goyal et al., 2007) y que incluso puede presentar problemas de estabilidad en la determinación de parámetros (Valdes-Parada et al., 2007). Otras líneas de trabajo apuntan a la elaboración de distintas estructuras matemáticas que puedan conducir al empleo de algoritmos de cálculo recursivos basados en el clásico mínimos cuadrados o esquemas de minimización por proyección ortogonal con ponderación de errores como las diversas variantes de filtros de Kalman (Baker et al., 2005; Yuceer et al., 2005; Bogaerts y Vande Wouwer, 2004; Cao y Gertler, 2004; Cao y Schwartz, 2004).

El objetivo de este trabajo fue realizar el ajuste óptimo de parámetros de modelos no lineales y predecir las condiciones iniciales de equipos empleados en ingeniería. La función objetivo empleada se elaboró como una relación entre las variaciones de los parámetros y de las condiciones iniciales con las variaciones de las variables de estado. La información empleada, así como la evolución temporal de la misma, fue obtenida experimentalmente durante la operación real del equipo de proceso. Se desarrolló una estructura de cálculo, sobre la base de los métodos de proyección ortogonal, para realizar a minimización. Este algoritmo fue codificado en lenguaje Fortran y empleado para realizar la evaluación buscada en un reactor del tipo tanque agitado continuo. Los resultados obtenidos fueron consistentes con el comportamiento físico del equipo.

DESARROLLO MATEMÁTICO

La formulación de las ecuaciones de balance, de las propiedades involucradas en los fenómenos que gobiernan los procesos, en término de variables de estado permite plantear una estructura matemática muy adecuada para abordar el análisis y diseño de estrategias óptimas de control de sistemas multivariables. Este conjunto de ecuaciones puede plantearse formalmente como:

(1)

Donde:

:

Vector de funciones no lineales de n componentes.

:

Vector de funciones, por lo general no lineales, de r componentes.

:

Vector de variables de estado de n componentes.

:

Vector de variables de entrada de m componentes.

:

Vector de variables de salida de r componentes.

:

Vector de parámetros de q componentes.

El comportamiento dinámico del sistema para todo tiempo t > 0 podrá estimarse si se conoce el estado X(t) = F (t , X(0) ,V(t) ,P). Una dinámica X(t) quedará unívocamente determinada, si se fija un estado particular X(0) para el tiempo adoptado como t = 0, y si se establecen para todo t > 0 los vectores V(t) y P. De toda esta información necesaria para poder calcular X(t), solo se puede disponer, en la mayoría de los casos de V(t), obtenida experimentalmente, y sólo en forma parcial. Por otro lado es imposible, físicamente, fijar como dato las condiciones X(0) y el vector P; es decir que para la determinación del estado tendremos n + q grados de libertad. Respecto de la dinámica del sistema real solo podremos disponer de la información que está contenida en el vector de variables medidas o de salida Y(t).

El ajuste entre la dinámica Y(t) = G (X (t)) correspondiente al modelo propuesto y la exhibida por el sistema físico Y(t) medido experimentalmente cada kT instante de muestreo con k = 1,2,3,...,K y T intervalo de muestreo, se realiza por medio de un proceso iterativo de búsqueda de extremo de una función escalar B(X( 0 / i + 1 ),P( i + 1)), con i índice de iteración, definida como:

(2)

Con:

 :

Vector de medidas experimentales, de (r K) componentes.

 :

Estima que correspondería al vector de variables de salida con (r K) componentes.

El algoritmo de cálculo que se desarrolla permite establecer una relación entre las variaciones en el vector de grados de libertad dU(i), y las correspondientes variaciones calculadas o estimadas del vector de variables de salida dŶM (KT/ i). En cada etapa i, de iteración, la minimización de la función, da por resultado un par X(0/i+1) y P (i+1), el proceso continúa hasta que se consigue que la norma del vector dŶM (KT/ i) sea menor o igual que un error preestablecido.

Algoritmo de cálculo

Dado un número real arbitrario e > 0 se definen las variaciones en los grados de libertad son X(0) y P, con dX(0/i) = X(0/i+1) - X(0 / i) y dP( i) = P(i+1) - P(i) , tal que sus normas verifiquen que:

½ ½ < e; ½ ½< e

El vector estado integrando con un método tipo Runge Kutta con un paso fijo h, en términos del tiempo, para los instantes de tiempo: t =(j+1) h  y  t = j h, será:

    

(3)

  

(4)

Donde: R1, R2, R3 y R4 son las pendientes que se ponderan en el método de integración. A partir de la ecuaciones (3) y (4) se puede obtener, en términos de variaciones:

(5)

Para R1, teniendo en cuenta su definición, se plantea:

  

(6)

Donde:

FX (X(jh/i), V(jh), P(i)): Jacobiano del vector de funciones F con respecto al vector estado X evaluado en t = jh en las condiciones que se corresponden con X(0/i) y P(i) y que estará simbolizado con X(jh/i). A esta matriz, de n x n, la denotaremos con MRX 1 (jh/i).

FP (X(jh/i), V(jh), P(i)): Jacobiano del vector de funciones F con respecto al vector de parámetros P evaluado en t = jh. Esta matriz será denotada por MRP 1(j h / i) y será de n x m.

O(e): Infinitésimo de orden superior al primero.

En forma análoga se puede obtener la expresión para las variaciones de las otras pendientes de tal manera de escribir:

(7)

Con:

(8)


(9)

En el algoritmo de cálculo empleado se hacen corresponder J pasos de integración h con cada período T de muestreo en tiempo real, generando así un conjunto de J ecuaciones. Por un procedimiento de sustitución progresiva se puede obtener, para el final de un período de muestreo T, una relación entre dX(Jh/i) y los grados de libertad dX(0/i) y dP(i):

 

(10)

Que se puede escribir como:

 

(11)

Donde:

    (12)


(13)

En virtud de la relación entre las variables de salida y las variables de estado, (1), se tendrá:

(14)

Donde GX (X(Jh/i)) es la matriz jacobiana de las funciones G con respecto al vector estado X evaluada en t =Jh, de r x n, esta relación reemplazando (11) será:

(15)


(16)

Con:

(17)

En la implementación computacional del algoritmo de cálculo, se debe asumir como dato al tiempo máximo de muestreo de los datos experimentales que exige la dinámica del sistema físico con el que se trabaja. Este tiempo, T = Jh se respeta ajustando los parámetros, paso de integración h y número de etapas de integración J que se eligen en función de la complejidad y extensión de todos los cálculos involucrados en el algoritmo propuesto. Considerando la ecuación (2) y haciendo uso de la relación (16), se puede plantear que la minimización de la función objetivo por elección de los grados de libertad, tomando k = 1,…, K conduce a:

(18)

Finalmente, a partir de esta condición de extremo, se puede encontrar

(19)

De esta manera termina una etapa del proceso iterativo. El mismo continúa con la incorporación de las medidas experimentales correspondientes a una nueva instancia de muestreo (K+1) y culminará cuando se alcance el error o tolerancia que se propone para la norma del vector dYM (KT/i).

RESULTADOS Y DISCUSIÓN

Para la implementación computacional del algoritmo propuesto se elaboró un programa en lenguaje Fortran y como aplicación se simuló la operación dinámica de un reactor tipo agitado continuo con un sistema de refrigeración. La reacción química considerada es del tipo A ® B, exotérmica, con cinética de primer orden e irreversible. Los datos adoptados son: volumen del reactor: 1000 litros; caudal de alimentación de reactivos: 10 litros/seg; densidad de reactivo: r =1 kilogramo/litro; calor específico de reactivo: Cp=1 Kcal/kgºK; concentración de reactivo en la alimentación: CA0 = 0.0065 mol/litro; temperatura de alimentación de reactivo: T0=340°K; temperatura media del fluido refrigerante: TR =430 °K; coeficiente de transferencia de calor multiplicado por el área total de intercambio térmico: UA=10 Kcal/seg°K; constante cinética: k=k0 exp(-DE/RT) con k0=7.86x1012 seg-1 y DE/R=14000 °K; calor de reacción: DH=-27000 Kcal/mol. Para estas condiciones de operación existe un solo estado estacionario: CA*=6.27396 x 10-5 mol/litro y T*=471.903 °K.

Los grados de libertad escogidos para su ajuste en la modelación propuesta son las condiciones iniciales: gl1=CA(0); gl2=T(0) y los parámetros gl3=DE/R y gl4=UA. Analizando la dinámica del sistema se observó que el tiempo de muestreo podía ser regulado en el intervalo de 1 a 10 seg. La primera etapa de minimización fue realizada tomando las primeras diez medidas de la temperatura de salida del reactor; mientras que las etapas posteriores acontecieron cada vez que se incorporó una nueva muestra. El control del proceso iterativo fue realizado por evaluación de la siguiente ecuación:

(20)

: vector de medidas de temperaturas en cada instancia de muestreo.

: vector valores temperaturas calculados en cada instancia de muestreo.

El criterio adoptado fue que los resultados alcanzasen un valor tal que se cumpliera que eK £ 10-5. Los distintos cálculos efectuados fueron programados de forma tal que se pudiera tener un conjunto lo más representativo posible de las situaciones dinámicas propias del equipo adoptado para la aplicación. Las iteraciones fueron ejecutadas para dos tipos de condiciones iniciales: condiciones iniciales de baja conversión, esto es gl1min=CA(0) > CA* y gl2min=T(0) < T* y condiciones iniciales de alta conversión, cuando gl1max=CA(0) < CA*  y  gl2max=T(0) > T*; para ambas condiciones se tomaron distintas combinaciones de gl3 y gl4. Los resultados obtenidos se muestran en la Tabla 1, para dos valores diferentes de partida en la búsqueda de los grados de libertad en la zona de baja conversión; y en la Tabla 2 para dos juegos de puntos iniciales diferentes en la zona de alta conversión.

Tabla 1: Zona de baja conversión. Punto de ajuste: gl1min =4.487 x10-4, gl2min=441.690 °K,
gl3min=14000 °K y gl4min=10 Kcal/seg °K.

Punto de partida de la iteración: gl1°=2.9x10-4 mol/l, gl2°=416.213 °K, gl3°=13000 °K, gl4°=12 Kcal/seg°K

Iteración

gl1

gl2

gl3

gl4

10

2.545857584177e-4

441.691023640174

13820.60094157880

32.016331385355

11

3.683340516425e-4

441.692282883219

14004.83221459412

-7.024746894916

12

4.814725534829e-4

441.683432873790

14035.37172233161

10.035305568008

13

4.523230283348e-4

441.689277367254

14003.51463901728

10.201600557591

14

4.486860362325e-4

441.689997494997

13989.99722443157

10.000621366976

15

4.487000048341e-4

441.690002442114

13999.99999437232

9.9999999748916

16

4.487000114669e-4

441.690002441404

13999.99999993420

9.9999999997497

Punto de partida de la iteración: gl1°=2.651x10-4 mol/l, gl2°=433.523 °K, gl3°=13500 °K, gl4°=8 Kcal/seg°K

Iteración

gl1

gl2

gl3

gl4

10

4.810836241233e-4

441.678419863300

14105.01237767562

24.512786175357

11

1.917417900092e-4

441.699579351681

13738.64170013089

16.331778193329

12

2.972560674153e-4

441.677342716786

13880.89406340208

1.150218758416

13

3.905949014983e-4

441.688448012662

13958.48267993342

8.040112306636

14

4.376721308023e-4

441.690139531258

13992.85803076330

9.605101868060

15

4.481954090461e-4

441.690020874886

13999.67344829840

9.982132464800

16

4.487988400917e-4

441.690002498734

13999.99923599795

9.999958664735

17

4.487000114576e-4

441.690002441406

13999.99999998461

9.999999998182

Tabla 2: Zona de alta conversión. Punto de ajuste: gl1max = 2.454 x10-5 mol/l, gl2max = 487.419 °K,
gl3max=14000 °K y gl4max=10 Kcal/seg °K.

Punto de partida de la iteración: gl1°=3.385x10-5 mol/l, gl2°=499.214 °K, gl3°=14500 °K, gl4°=8 Kcal/seg°K

Iteración

gl1

gl2

gl3

gl4

10

3.508657569994e-5

487.422415644243

14180.54746603966

9.531970610260

11

2.216248635500e-5

487.422490814547

13997.33793605492

10.024559084384

12

2.457443981390e-5

487.418440590160

14000.14284631703

9.999867450517

13

2.454602086192e-5

487.419005933996

14000.00025685834

9.999999810355

14

2.454599613268e-5

487.419006347595

13999.99999981399

10.000000000108

Punto de partida de la iteración: gl1°=5.533x10-5 mol/l, gl2°=507.781 °K, gl3°=15000 °K, gl4°=12 Kcal/seg°K

Iteración

gl1

gl2

gl3

gl4

10

6.251112282296e-5

487.435311277178

14486.86348702014

9.949266674404

11

7.341393035087e-5

487.425867986360

14556.54531081021

9.246729590537

12

-1.616091517275e-5

487.443352933337

14053.32988358801

10.380540074268

13

3.142063920302e-5

487.356583385624

14085.64380902378

9.924377924570

14

2.611160321941e-5

487.404125412053

14024.46753864051

9.984547798017

15

2.461013099511e-5

487.418052685434

14001.16388216993

9.999681762250

16

2.454606987239e-5

487.419005494171

14000.00222897587

9.999999316394

17

2.454599963891e-5

487.419006347611

14000.00000007900

9.999999999720

El hecho de considerar las leyes físicas que gobiernan los fenómenos que se producen en las unidades de proceso, garantiza que la dinámica obtenida por cálculo a partir de tales modelos sea consistente con la información que se puede adquirir durante el funcionamiento real de los equipos. Esto puede ser comprobado a partir de los resultados obtenidos. En todos los casos se alcanza el valor de ajuste pero no siempre de la misma forma. En la zona de alta conversión, Tabla 2, prácticamente no hay restricción sobre el valor que se tome de punto de partida , siempre se alcanzó el valor buscado, esto coincide con el hecho de que el ejemplo elegido solo posea un estado estacionario y éste sea el de alta conversión, es decir con valores de temperatura mayores al estacionario. Por otro lado la zona de baja conversión, Tabla 1, no presenta la misma flexibilidad en la selección de los valores de partida, el entorno con el cual se puede trabajar es mucho menor que en el caso anterior, esto coincide con la situación física de seleccionar valores que lleven a la situación de reactor apagado, es decir cuando la reacción aún no ha empezado a suceder. Los valores extraordinarios, como por ejemplo menores que cero, alcanzados en las primeras iteraciones para la concentración inicial y en el coeficiente global de transferencia de calor pueden ser atribuidos al método de integración numérica de las ecuaciones no lineales, de cualquier forma al avanzar en la incorporación de datos experimentales esto se corrige favorablemente.

En cualquier caso el algoritmo de ajuste desarrollado encontró una única solución con el grado de error que se le estableció como criterio de parada. No se presentaron casos de múltiple o ninguna solución, situación que si bien es posible matemáticamente, en este equipo no es factible desde un punto de vista físico.

CONCLUSIONES

Los resultados obtenidos permiten concluir que el algoritmo desarrollado arroja en todos los casos una solución única, en la determinación de condiciones iniciales y de parámetros, con total independencia de los puntos de partida adoptados en las iteraciones; Esto es plenamente consistente con el basamento físico y termodinámico del modelo matemático con el que se trabaja. En esta cuestión el algoritmo de cálculo que se propone ofrece una notoria ventaja frente a algunos otros presentados en trabajos donde se plantea la existencia de más de una solución para la estimación de parámetros.

REFERENCIAS

Babu, B.V. y K.N. Sastry; Estimation of heat transfer parameters in a trickle-bed reactor using differential evolution and orthogonal collocation, Computers & Chemical Engineering: 23(3), 327-339 (1999).        [ Links ]

Baker, C.T.H. y otros tres autores; Computational modelling with functional differential equations: Identification, selection, and sensitivity, Applied Numerical Mathematics: 53(2-4), 107-129 (2005).        [ Links ]

Bogaerts, P. y A.V.A. Vande Wouwer; Parameter identification for state estimation-application to bioprocess software sensors,Chemical Engineering Science: 59(12), 2465-2476 (2004).        [ Links ]

Cao, J. y J. Gertler; Noise-induced bias in last principal component modeling of linear system, Journal of Process Control: 14(4), 365-376 (2004).        [ Links ]

Cao. L. y H.M. Schwartz; Analysis of the Kalman filter based estimation algorithm: an orthogonal decomposition approach, Automatica: 40 (1), 5-19 (2004).        [ Links ]

Cheng, Z. y W. Yuan; Initial estimation of heat transfer and kinetic parameters of a wall-cooled fixed-bed reactor, Computers & Chemical Engineering: 21 (5), 511-519 (1996).        [ Links ]

Cheng, Y.S., C.F. Abi y L.S. Kershenbaum; On-line estimation for a fixed-bed reactor with catalyst deactivation using nonlinear programming techniques, Computers & Chemical Engineering: 20 Supp. 2, S793-S798 (1996).        [ Links ]

Damak T.; Procedure for asymptotic state and parameter estimation of nonlinear distributed parameter bioreactors, Applied Mathematical Modelling: 31 (7), 1293-1307 (2007).        [ Links ]

Goyal, A. y otros tres autores; Application of orthogonal collocation and regression techniques for recovering parameters of a two-pathway transdermal drug-delivery model, Computers & Chemical Engineering: 31( 3), 107-120 (2007)        [ Links ]

Karacan, S. y otros tres autores; Pole placement self tuning control for packed distillation column, Chemical Engineering and Processing: 36 (4) 309-315 (1997).        [ Links ]

Tsang, T.H., T.F. Edgar y J.O. Hougen; Estimation of heat transfer parameters in a packed bed, The Chemical Engineering Journal: 11 (1), 57-66 (1976).        [ Links ]

Valdes-Parada F.J. y otros tres autores; On Green’s function methods to solve nonlinear reaction–diffusion systems, Computers & Chemical Engineering: 32(3), 503-511 (2007).        [ Links ]

Yuceer, M., A. Ilknur y B. Ridvan; An integration based optimization approach for parameter estimation in dynamic models, Computer Aided Chemical Engineering: 20(1), 631-636 (2005).        [ Links ]