Los métodos matemáticos de optimización tienen diversas aplicaciones en la ingeniería de procesos, que abarcan desde el diseño óptimo a la optimización del rendimiento económico. En esta entrada, se ilustra el uso de estas metodologías para optimizar el comportamiento dinámico de un reactor durante su arranque. Para ello, se formula el modelo dinámico del reactor y la optimización de su comportamiento durante el arranque como un problema de optimización convexo, que se resuelve empleando la toolbox CVX de Matlab.
1. Introducción
Los reactores de tanque agitado, como el mostrado esquemáticamente en la Figura 1, son unos de los equipos más utilizados industrialmente para llevar a cabo reacciones químicas. Estos reactores son dispositivos sencillos y fáciles de operar, que pueden operar tanto en discontinuo (por cargas) como en continuo (con un flujo continuo de reactivos y productos), y por ello son muy utilizados en industrias pequeñas y medianas.

Figura 1: Diagrama esquemático de un reactor continuo de tanque agitado
En una configuración para operación en continuo, como la mostrada en la Figura 1, el reactor consta de al menos una tubería de entrada mediante la que se introduce un flujo de reactivos, y una tubería de salida mediante la que se extrae un flujo de productos. El reactor en sí es un tanque de un volumen determinado; su parámetro característico es el tiempo de residencia del fluido en su interior, τ:
Siendo V el volumen del tanque y F el flujo de entrada de reactivos y de salida de productos (que en régimen estacionario son iguales para que el nivel de fluido en el interior del tanque se mantenga constante), el tiempo de residencia τ debe ser el adecuado para que se complete la reacción química que se desee llevar a cabo. Otros elementos habituales del reactor son un agitador que mezcla y homogeneiza el fluido contenido en él, y otros elementos auxiliares como un encamisado o un serpentín interior para ajustar la temperatura del fluido o elementos mecánicos que favorezcan la mezcla del fluido.
Los reactores de tanque agitado que operan en continuo se diseñan para que proporcionen el rendimiento deseado en su operación en estado estacionario, es decir, una vez se alcanzan condiciones de operación que son invariables con el tiempo. Brevemente, el diseño consiste en calcular el tiempo de residencia τ apropiado para que la reacción se produzca del modo deseado [1]. Las fases no estacionarias del funcionamiento del reactor, como el arranque y la parada, se estudian por lo general en menor detalle. Sin embargo, estas fases también tienen gran importancia tanto desde el punto de vista económico (minimizar la cantidad de productos/reactivos que no se producen de acuerdo con las especificaciones y que por lo tanto hay que desechar, con el consiguiente coste económico), como por los aspectos de seguridad (evitar posibles riesgos como picos de concentración/temperatura que pueden llevar a reacciones descontroladas y posibles accidentes).
Con este planteamiento, el objetivo de este trabajo es obtener una estrategia de arranque óptima para un reactor continuo de tanque agitado. Esta estrategia se obtiene resolviendo un problema de optimización convexo.
2. Metodología
2.1 Modelo del reactor
Se considera un reactor continuo de tanque agitado de un determinado volumen V y con un flujo determinado de reactivos y productos F que opera isotérmicamente, es decir, a temperatura constante. Para analizar un caso de una cierta complejidad, se considera una reacción cuyo mecanismo consta de dos etapas en paralelo en las que un reactivo A se convierte en sendos productos B y C:
En este mecanismo, se supone que B es el producto deseado, mientras que C es un subproducto no deseado. Se considera que ambas reacciones son de primer orden respecto del reactivo, con constantes cinéticas k1 y k2, respectivamente, con lo que sus velocidades de reacción vienen dadas por:
donde cA es la concentración del reactivo A. En una primera aproximación, para realizar un modelo sencillo de un reactor como este, se suele suponer que la mezcla que se alcanza en su interior es perfecta, de modo que todos los puntos internos del reactor tienen las mismas condiciones; en particular, en todos ellos hay la misma concentración de reactivo, cA. Además, como la corriente de productos se extrae de este volumen contenido en el reactor, la concentración en la corriente de salida también coincide con la concentración en cualquier punto del interior del reactor (ver Figura 1). En estas condiciones, la dinámica del reactor viene descrita por un balance de materia, que para las tres sustancias implicadas en el mecanismo de reacción supuesto tiene la siguiente forma [1]:
donde cA,0 es la concentración de reactivo que hay inicialmente en el reactor y se supone que las concentraciones iniciales de los dos productos son cero.
Para facilitar la resolución del problema y su implementación como un problema de optimización convexa, la solución de este sistema de ecuaciones diferenciales se ha aproximado mediante un método numérico de Euler explícito [2]. Con este método, las soluciones aproximadas se obtienen mediante las siguientes ecuaciones:
donde Δt es el paso de tiempo empleado en el método numérico de Euler. Desarrollando estas fórmulas iterativas junto con las condiciones iniciales, las concentraciones de salida se pueden calcular en función de la concentración de entrada de reactivo cA,e mediante las siguientes expresiones en forma matricial:
En estas expresiones, si se define f como la agrupación de constantes:
las matrices para el cálculo de la concentración del reactivo cA vienen dadas por:
siendo n el número de pasos de integración realizados mediante el método de Euler. Del mismo modo, definiendo la agrupación de constantes g como:
Las matrices para el cálculo de la concentración del producto B son las siguientes:
Las expresiones para la concentración del producto C son equivalentes, reemplazando k1 por k2.
2.2 Planteamiento del problema de optimización
El cálculo de la estrategia de arranque se plantea como un problema de optimización con un doble objetivo: obtener una evolución de las concentraciones durante el arranque que se aproxime lo más posible a la evolución deseada, y minimizar el coste del arranque en el sentido de minimizar el gasto de reactivo empleado en este proceso. Para ello, la variable manipulada que se optimiza es la concentración de entrada suministrada al reactor a lo largo del tiempo, cA,e.
Respecto del primer objetivo, se supone que el objetivo que se busca alcanzar es que la concentración de salida del producto deseado B sea igual a una determinada concentración cB,set, lo que se puede efectuar interpretando este objetivo como un problema de búsqueda de la mejor aproximación que minimice la distancia entre las funciones deseada y obtenida de la concentración. Así, esta especificación se formula a través de la norma 2 de la diferencia entre la concentración de salida del reactor y esta concentración deseada:
Respecto del primer objetivo, el coste del arranque se minimiza a través de la minimización de la concentración de entrada suministrada cA,e, con el matiz de que esta concentración (o el coste) no se puede disminuir hasta cero puesto que para que el reactor opere produciendo B según se desea, es preciso suministrarle una cierta cantidad de reactivo A. Esta cantidad se puede calcular resolviendo los balances de materia en estado estacionario:
de las cuales se deduce que la concentración de entrada requerida en el estado estacionario es:
El coste total del reactivo empleado se obtiene a partir de la suma en cada intervalo de tiempo de la cantidad de reactivo empleado, obtenida como el producto de su concentración por el flujo, multiplicada por el conste unitario de ese reactivo, K:
Con esto, y considerando que en el modelo propuesto el flujo F es constante y que en el método de Euler empleado Δt es también constante, y que por lo tanto ambas constantes pueden agruparse con K, la segunda especificación se formula minimizando la norma 1 de la desviación de la concentración de entrada respecto de esta concentración requerida:
Se tiene así un problema de optimización con dos componentes, en el que ambas componentes son convexas, que se reducen a una única función convexa mediante una combinación lineal de estas dos componentes a través de un parámetro de peso arbitrario γ, obteniéndose así los puntos óptimos de Pareto mediante el método de la escalarización [3]. Para el planteamiento del problema debe tenerse en cuenta también que las componentes del vector cA,e no pueden tomar cualquier valor, sino que están limitadas, inferiormente por 0 por consideraciones físicas (la concentración no puede ser negativa), y superiormente por un cierto valor cA,e,max por consideraciones prácticas (en la práctica no se podrá incrementar la concentración indefinidamente, sino solo hasta un valor máximo dado). Con esto, el problema de optimización se formula como:
minimizar2.3 Resolución del problema de optimización
El problema de optimización se ha resuelto mediante la toolbox CVX de Matlab [4, 5], empleando para ello el código incluido como Anexo.
La toolbox CVX resuelve el problema de optimización empleando el algoritmo SPDT3 [6]. Este algoritmo es del tipo Primal-dual interior point, o método primal-dual de punto interior [3]. Brevemente, este método, adecuado para resolver problemas con restricciones de desigualdad, como el planteado, puede considerarse como una extensión del método de la barrera, que mejora el rendimiento y la rapidez de convergencia de este método [3].
Como en otros algoritmos para la resolución de problemas de optimización con restricciones de desigualdad, el punto de partida son los algoritmos para problemas con restricciones de igualdad, que pueden resolverse, entre otros métodos, mediante un método de Newton modificado que, partiendo de una estimación inicial factible de la solución, emplee las condiciones de Karun-Kush-Tucker KKT para determinar una dirección de optimización que resulte factible. Esto se lleva a cabo aproximando la función a optimizar por su extensión de Taylor de segundo orden, que junto con las restricciones de igualdad se puede resolver analíticamente a través de sus condiciones KKT [3].
Los métodos de punto interior parten de este resultado para resolver problemas con restricciones de desigualdad a partir de una secuencia de problemas con restricciones de igualdad, que se resuelven mediante el método de Newton [3]. Así, el método de la barrera consiste en hacer implícitas las restricciones de desigualdad, incorporándolas a la función objetivo mediante una función de barrera que penaliza el incumplimiento de las restricciones de desigualdad; el problema se resuelve a través de una secuencia de funciones barrera cada vez más restrictivas, hasta alcanzar una solución con la precisión deseada. Así, un problema de la forma:
minimizarse transforma en un problema sin restricciones empleando la función de barrera logarítmica:
minimizarde modo que a medida que se incrementa t, la solución de este segundo problema converge a la del problema original. Este problema se resuelve mediante el método de Newton aplicado a las condiciones KKT del problema.
En un método primal-dual de punto interior se modifica el método de la barrera en lo referente a la búsqueda de la dirección de optimización en la que se aplica el método de Newton, la cual se obtiene simultáneamente para las variables del problema original y las de su problema dual. Así, partiendo de las condiciones KKT para el problema de la barrera:
solución en la que las direcciones de búsqueda del problema primal y dual están acopladas y dependen una de la otra.
El algoritmo SPDT3 es una versión más elaborada de este método básico primal-dual de punto interior en el que se implementa una forma mejorada de obtener las direcciones de optimización para la resolución del problema por el método de Newton. Concretamente, el algoritmo contempla dos formas de obtener estas direcciones, denominadas dirección HKM y dirección NT [6]. Estos métodos se diferencian en las matrices de escalado empleadas para resolver el cálculo de la dirección de optimización. En el problema presentado en este ejemplo, se ha empleado el método NT.
3. Resultados
El problema de optimización se ha resuelto adoptando los siguientes valores numéricos para las constantes del modelo, que pueden considerarse de un orden de magnitud típico para un sistema de estas características [1]:
Por otra parte, se ha supuesto que el reactor está cargado con una concentración inicial de reactivo cA,0 = 2.5 mol/L y que la concentración máxima que puede tener la corriente de entrada es cA,e,max = 10 mol/L. La concentración deseada de producto se ha tomado como cB,set = 2 mol/L.
La Figura 2 muestra los resultados obtenidos para diversos valores del factor de peso γ. La figura 2a corresponde a γ = 0, es decir, a un problema en el que el objetivo es arrancar el reactor de modo que se alcance la concentración de salida deseada con la mayor rapidez posible. Como puede apreciarse, la estrategia proporcionada por la optimización consta de tres pasos: en un primer paso, se introduce la alimentación con la mayor concentración de entrada posible durante un determinado tiempo, y a continuación la concentración de la corriente de entrada se baja a cero durante un segundo intervalo de tiempo. La longitud de estos intervalos propicia que, tras dejar la concentración de entrada a cero, la inercia del sistema haga que la concentración de producto siga subiendo hasta alcanzar el valor deseado. Tras esto, la concentración de entrada se fija en el valor correspondiente al estado estacionario, cA,e,ss. Al comienzo de esta tercera etapa, puede apreciarse una pequeña oscilación en la concentración de la corriente de entrada.

Figura 2: Resultados obtenidos con diferentes valores de γ. (a) γ = 0, (b) γ = 10-5, (c) γ = 10-3, (d) γ = 0.03.
La Figura 2d corresponde a la situación opuesta de un valor muy alto de γ = 0.03, que hace que la optimización busque minimizar el coste de emplear valores extremos de la concentración de entrada sin atender a la rapidez del arranque. Como puede observarse, en este caso la solución óptima se reduce a aplicar desde un comienzo la concentración de entrada correspondiente al estado estacionario, cA,e,ss, de modo que el arranque es mucho más lento y de hecho el estado estacionario aún no se ha alcanzado en el tiempo de simulación considerado.
Las Figuras 2b y 2c corresponden a casos intermedios, con γ = 10-5 y γ = 10-3, respectivamente. Puede observarse que a medida que se incrementa γ y, por lo tanto, el peso del coste en la función objetivo, la duración de las dos etapas iniciales de concentración de entrada máxima y mínima se acorta, reemplazándolas por un periodo más largo en el que se aplica la concentración de estado estacionario. Con γ = 10-5 se observa una pequeño periodo en el que se aplica esta concentración de estado estacionario intercalado entre los periodos de concentración máxima y mínima, y con γ = 10-3 la etapa de concentración 0 se elimina por completo. Estas variaciones resultan en un arranque algo más lento.
Esta evolución del comportamiento puede apreciarse también en la Figura 3, que presenta la curva de compensación entre las dos componentes de la función objetivo. En esta figura se observa que la curva presenta un pequeño codo, que puede interpretarse como el punto que mejor balancea las dos componentes de la función objetivo [3]. Este punto está ubicado aproximadamente a un valor γ = 10-3, que se corresponde con el comportamiento representado en la Figura 2c.
Figura 3: Curva de compensación (trade-off curve) del problema de optimización.
Debe resaltarse por último que la metodología propuesta es general, en el sentido de que puede aplicarse a cualquier ruta de arranque deseada. Como ilustración de este hecho, en la figura 4 se muestra el resultado del algoritmo de optimización para una ruta de arranque en rampa y para un proceso en el que se produzcan varios cambios sucesivos en escalón en la concentración de salida deseada. Como puede apreciarse, el algoritmo de optimización proporciona el mejor ajuste (dentro de las limitaciones impuestas por la dinámica del sistema) a estas situaciones.

Figura 4: Estrategias óptimas de arranque para un arranque en rampa y para cambios de objetivo en escalón
4. Conclusiones
Se ha resuelto el problema del cálculo de la estrategia óptima de arranque de un reactor continuo de tanque agitado, que proporcione el mejor ajuste a la respuesta deseada del reactor y el menor coste, planteando para ello un problema de optimización convexa vectorial y resolviéndolo mediante la toolbox CVX de Matlab. Se han analizado las condiciones que proporcionan el mejor balance entre las dos componentes de la función objetivo, seleccionando unas condiciones óptimas que proporcionan un arranque rápido sin oscilaciones prolongadas en la concentración de entrada.
Como propuestas de ampliación o mejora futura del trabajo, puede sugerirse en primer lugar la incorporación del balance de energía en el modelo del reactor, haciéndolo no-isotérmico. Las variaciones de temperatura complican notablemente el comportamiento dinámico del reactor, ya que si por ejemplo la reacción produce calor, que hace que aumente la temperatura en el reactor, el aumento de temperatura propicia que se aceleren las cinéticas de reacción, produciéndose a su vez aún más calor y entrándose así en un proceso de realimentación que es responsable de muchos de los problemas de operación y de seguridad que pueden darse en el arranque de uno de estos reactores. Una segunda posible ampliación consistiría en tener en consideración las perturbaciones que inevitablemente afectan a un equipo de estas características, como pueden ser las perturbaciones en las concentraciones, flujos o temperaturas de las corrientes o el reactor. Estas perturbaciones no son conocidas a priori, pero se pueden describir mediante variables aleatorias.
Anexo: Código de Matlab
% Controla la concentración de salida cA de un reactor de tanque agitado % isotérmico con una reacción de orden 1 A ->B, A->C variando la concentración % de entrada ce % Modelo dinámico del reactor: % dcA/dt = -k*cA + (ce-cA)/tr % Modelo discretizado (método de Euler explícito): % cAi = cA(i-1) + dt·(-k*cA(i-1) + (cei-1)-cA(i-1))/tr) clear tr = 300; % Tiempo de residencia en el reactor k1 = 0.05; % Constante cinética de la reacción 1 k2 = 0.025; % Constante cinética de la reacción 2 cA0 = 2; % Concentración inicial en el reactor cB0 = 0; cC0 = 0; cBset = 2; % Concentración de salida objetivo en el reactor dt = 1; % Paso de tiempo de la discretización n = 500; % Número de pasos de tiempo gamma = 2e-2; % Parámetro de peso cass = cBset/(k1*tr); % Concentración de A en estado estacionario cess = cass*(1+(k1+k2)*tr); % Concentración de entrada de A en estado estacionario % Construye matriz del sistema discretizado, cA = A*ce + b k = k1 + k2; AcA = zeros(n,n); bcA = zeros(n,1); f = (1 - dt*k - dt/tr); bcA(1) = cA0; for i = 2:n for j = 1:i-1 AcA(i,j) = dt/tr*f^ (i-j-1); end bcA(i) = f^ (i-1)*cA0; end % Construye matriz del sistema discretizado, cB = A*ce + b AcB = zeros(n,n); bcB = zeros(n,1); fB = (1 - dt/tr); bcB(1) = cB0; bcB(2) = fB*cB0 + k1*dt*cA0; for i = 3:n for j = 1:i-2 sum = 0; for indxsum = 2:i-1 sum = sum + fB^ (i-indxsum-1)*k1*dt* AA(indxsum,j); end AcB(i,j) = sum; end sum = 0; for indxsum = 1:i sum = sum + fB^ (i-indxsum)*k1*dt*bcA(i); end bcB(i) = f^ (i-1)*cB0 + sum; end % Construye matriz del sistema discretizado, cC = A*ce + b AcC = zeros(n,n); bcC = zeros(n,1); fB = (1 - dt/tr); bcC(1) = cC0; bcC(2) = fB*cC0 + k2*dt*cA0; for i = 3:n for j = 1:i-2 sum = 0; for indxsum = 2:i-1 sum = sum + fB^ (i-indxsum-1)*k2*dt* AcA(indxsum,j); end AcC(i,j) = sum; end sum = 0; for indxsum = 1:i sum = sum + fB^ (i-indxsum)*k2*dt*bcA(i); end bcC(i) = f^ (i-1)*cC0 + sum; end % Resuelve optimización cvx begin variable ce(n) minimize(norm(AcB*ce+bcB-cBset) + gamma*(norm(ce-cess,1))) subject to 0 <= ce <= 10 cvx end %Representa resultados cA = AcA*ce + bcA; cB = AcB*ce + bcB; cC = AcC*ce + bcC; hold off plot(cA(1:n-2),’-b’) hold on plot(cB(1:n-2),’-g’) plot(cC(1:n-2),’-r’) plot(ce(1:n-2),’-k’) xlabel (’t(s)’) ylabel (’c(mol/L)’); legend(’cA’,’cB’,’cC’,’cAe’);
Referencias
[1] O. Levenspiel. El omnilibro de los reactores químicos. Ed. Reverté, 1986.
[2] C. Moreno González. Introducción al Cálculo Numérico. Ed. UNED, 2011.
[3] S. Boyd y L. Vandenberghe. Convex Optimization. Ed. Cambridge University Press, 2004.
[4] CVX Research, Inc. CVX: Matlab software for disciplined convex programming, version 2.0.,http://cvxr.com/cvx, April 2011.
[5] M. Grant y S. Boyd. Graph implementations for nonsmooth convex programs, Recent Advances in Learning and Control (a tribute to M. Vidyasagar), V. Blondel, S. Boyd y H. Kimura, editores, pp. 95-110, Lecture Notes in Control and Information Sciences, Springer, 2008.
[6] R. H. Tütüncü, K. C. Toh y M. J. Todd. Solving semidefinite-qadratic-linear programs using SDPT3. Mathematical Programming, Series B, 95:189-217, 2003.