8. SISTEMAS NO LINEALES 8.1 LINEARIZACION DE UNA FUNCION La mayoría de los procesos químicos son no lineales porque las ecuaciones diferenciales que los modelan contienen términos, en función de las variables de estado, que no son lineales. Es posible desarrollar la solución matemática de un modelo no lineal en el dominio del tiempo mediante métodos analíticos o numéricos y mucho más con la ayuda de lenguajes de programación matemática como Matlab. Para algunos propósitos, principalmente de simplificación, es una práctica muy usual aproximar el modelo a una forma lineal que facilita su solución y permite simularlo en el dominio de Laplace o compactarlo en la forma del espacio de los estados. La linearización es una técnica utilizada para aproximar la respuesta de un sistema no lineal a ecuaciones diferenciales lineales. La aproximación lineal es válida en una región próxima a un punto base alrededor del cual se haga la linearización. Para facilitar la manipulación de las ecuaciones linearizadas, se seleccionará el estado estacionario inicial como el punto base para la linearización. La linearización de una función no lineal con una sola variable, f (x) , consiste en expandir dicha función en una Serie de Taylor como se muestra en (8.1) alrededor de un punto base, x , y despreciar todos los términos de grado superior a uno. La ecuación (8.2) corresponde a la expansión simplificada de la función, f (x) , a una expresión lineal en términos de x.

1  ∂ 2 f ( x)   ∂f ( x)   ( x − x )2 + ... f ( x) = f ( x ) +   ( x − x ) +  2 2!  ∂x  x  ∂x  x  ∂f ( x )  f ( x) = f ( x ) +   ( x − x ) + ...  ∂x  x

(8.1)

(8.2)

La fórmula (8.2) es la básica para la linearización de una función de una sola variable. Como x es una constante, el miembro derecho de la ecuación (8.2) es lineal con respecto a la variable x.

130

Una función no lineal con dos o más variables se puede expandir en una Serie de Taylor en una forma similar a la ecuación (8.1), de tal manera que la aproximación linearizada es de la forma:

f (x1, x2 ,...xn ) = f (x1 , x2 ,...,xn ) +

n

∑ i=1

 ∂f (x1, x2 ,...,xn )    (xi − xi ) (8.3) x ∂ i ( x1,x2 ,...,xn ) 

En la expansión (8.3) se incluye un término lineal por cada una de las variables en la función

8.2 REACCION DE VAN DE VUSSE En este módulo se considera el sistema de reacciones de Van de Vusse para producir un compuesto B a partir de un componente A porque es un ejemplo de un esquema de reacciones que muestran una máxima concentración de producto para un determinado flujo a través del reactor. El esquema es el siguiente: K1

K2

A→ B →C K3

2A→ D Este sistema fue presentado por Van de Vusse en 1964 y muestra diferencias en su respuesta según sean las condiciones de operación seleccionadas. Se ha demostrado que varios procesos químicos se comportan de acuerdo al sistema propuesto por Van de Vusse. El proceso se desarrolla en un reactor de mezcla completa como el mostrado en la figura 8.1, que opera isotérmicamente, es decir, que no se necesita un balance de energía; y que los parámetros de velocidad de reacción son constantes. Se considera que la densidad es constante a través del sistema, como también el volumen de masa reaccionante en el tanque. Además, se asume que la corriente de alimentación solo contiene el componente A. Las reacciones 1 y 2 son de primer orden y la reacción 3 es de segundo orden. Inicialmente se analiza la respuesta del sistema en condiciones estacionarias para observar la variación de las concentraciones de A y B con el cambio en la velocidad espacio, Fs / V , y asignando los siguientes valores para los parámetros de las constantes de velocidades de reacción y la concentración de A en el alimento Mach

131

C Af = 10

gmol litro

K1 =

5 min −1 6

F CAf(t)

K2 =

5 min −1 3

V

K3 =

1 litro / mol − min 6

F CA(t)

Figura 8.1 Reactor de Mezcla Completa Luego se analiza la respuesta dinámica del sistema y finalmente se plantean tres casos para observar los efectos que tiene en la respuesta del sistema el cambio en el valor de la velocidad espacio

Modelo No Lineal – Dominio Tiempo Un balance de materia del componente A en el reactor, a volumen constante, es: d (C A ) F = (C Af − C A ) − K 1C A − K 3 C A2 dt V

(8.4)

Un balance de materia del componente B en el reactor, a volumen constante, es: d (C B ) F = C B + K 1C A − K 2 C B dt V

(8.5)

Un balance de materia del componente C en el reactor, a volumen constante, es:

Mach

132 d (C C ) F = − CC + K 2 C B dt V

(8.6)

Un balance de materia del componente A en el reactor, a volumen constante, es: d (C D ) 1 F = − C D + K 3 C A2 2 dt V

(8.7)

Las ecuaciones (8.4), (8.5) (8.6) y (8.7) son no lineales porque incluyen términos no lineales como el producto de las variables flujo y concentración y otro de grado dos con respecto a la concentración de A. Por lo tanto, el modelo es no lineal y se desarrollará su solución como tal y se comparará con la del correspondiente modelo linearizado

Análisis en estado estacionario Las ecuaciones (8.4) a (8.7) se pueden combinar y deducir las expresiones matemáticas para calcular las concentraciones de los componentes del sistema de reacciones en el estado estacionario, así:

C As

F  −  K1 + s V =  2K 3

C Bs =

C Cs =

C Ds

Mach

K 1C As Fs + K2 V K 2 C Bs Fs V

2 K 3 C As = F 2 s V

2

F   F    K 1 + s  + 4 K 3  s C Af  V  V  +  2K 3

(8.8)

(8.9)

(8.10)

(8.11)

133

En las Figuras 8.2 y 8.3 se grafican las variaciones de las concentraciones de A y B con el cambio en la velocidad espacio y se observa un máximo en la concentración de B.

Figura 8.2 Concentración de A en función de la Velocidad Espacio

Figura 8.3 Concentración de B en función de la Velocidad Espacio Mach

134

Es decir, si el objetivo del proceso es la maximización de la producción de B, entonces existe una velocidad espacio óptima de 1.2921 min-1 para alcanzar dicho máximo con un valor de 1.2660 mol/litro (Ver archivos conmax.m y vandevusse.m).

Análisis en estado dinámico – Modelo No Lineal Las figura (8.4) muestran las respuestas dinámicas del sistema de reacciones para los componentes A y B asignando un valor para la velocidad espacio de 4/7 min-1 en estado estacionario y un cambio paso de 0.01 min-1.

Figura 8.4 Concentración de A y B para un cambio paso en la velocidad espacio Se observa un aumento monotónico para el componente A y una Respuesta Inversa para el componente B, es decir, que durante un periodo inicial de tiempo la variación de la concentración de B es inversa a la tendencia posterior que toma dicha variación. Un aumento en la velocidad espacio con una determinada concentración en A mayor que la que se tiene en el reactor hace que la concentración de A aumente pero con respecto a B ocurre una dilución que hace que durante un tiempo determinado la concentración disminuya y con ello su velocidad de formación de B.

Mach

135

Después de un tiempo el cambio en la concentración de B toma el perfil que se espera debido al aumento en la concentración del reaccionante A que lo produce. La figura (8.5) muestra las variaciones de las concentraciones de C y D. Se observa una disminución en la concentración de C y una Respuesta Inversa en la variación de la concentración del componente D

Figura 8.5 Concentración de C y D para un cambio paso en la velocidad espacio

Modelo Lineal – Dominio Tiempo – Espacio de los Estados Los miembros derechos de las ecuaciones (8.4) a (8.7) son funciones no lineales de las concentraciones y la velocidad espacio y se pueden considerar de la siguiente manea:. dC A F F = f 1 (C A , C B , C C , C D , ) = (C Af − C A ) − K 1C A − K 3 C A2 dt V V dC B F F = f1 (C A , C B , C C , C D , ) = − C B + K 1C A − K 2 C B dt V V

Mach

136

d (C C ) F F = f 3 (C A , C B , C C , C D , ) = − C C + K 2 C B dt V V d (C D ) F F 1 = f 4 (C A , C B , C C , C D , ) = − C D + K 3 C A2 dt V V 2

La linearización de estas funciones con respecto a cada una de ellas permite expresar a dichas ecuaciones en términos lineales de las correspondientes variables desviación de la siguiente forma: dC A F = A11C A + A12 C B + A13 C C + A14 C D + B11 dt V

(8.12)

dC B F = A21C A + A22 C B + A23 C C + A24 C D + B21 dt V

(8.13)

dC C F = A31C A + A32 C B + A33 C C + A34 C D + B31 dt V

(8.14)

dC A F = A41C A + A42 C B + A43 C C + A44 C D + B41 dt V

(8.15)

Siendo, A11 =

F ∂f1 = − s − K 1 − 2 K 3 C As ; ∂C A V

A21 =

∂f 2 = K1 ; ∂C A

A22 =

F ∂f 2 = − s − K2 ; ∂C B V

A31 =

∂f 3 =0; ∂C A

A32 =

∂f 3 = K2 ; ∂C B

A41 =

∂f 4 ∂f = K 3 C As ; A42 = 4 = 0 ; ∂C A ∂C B

Mach

A12 =

∂f 1 =0; ∂C B

A13 =

A23 =

∂f 1 =0 ∂C C

∂f 2 = 0; ∂C C

A14 =

∂f 1 =0 ∂C D

A24 =

∂f 2 =0 ∂C D

A33 =

F ∂f 3 =− s ; V ∂C C

A34 =

∂f 3 =0 ∂C D

A43 =

∂f 3 = 0; ∂C C

A44 =

F ∂f 3 =− s V ∂C D

137

B11 =

∂f1 = C Af − C As ; ∂( F / V )

B21 =

∂f 2 = −C Bs ; ∂( F / V )

B31 =

∂f 3 = −C Cs ∂( F / V )

B41 =

∂f 4 = −C Ds ∂( F / V )

Para la escritura del modelo en la forma del espacio de los estados, las matrices A, B, C y D son las siguientes:  Fs − V − K 1 − 2 K 3 C AS    K1   A=  0     K 3 C As   1   0 C= 0   0

0

−

0

Fs − K2 V K2

0

0

−

Fs V 0

0 0 0   1 0 0  0 1 0   0 0 1

 0     0     0    Fs  −  V 

C Af − C As       − C Bs    B=   − C Cs       − C Ds 

0      0  D=    0      0 

Análisis en el Espacio de los Estados Para el análisis dinámico en el espacio de los estados y la comparación entre el modelo no lineal y el lineal se consideran tres casos en los cuales se asigna una velocidad espacio menor, igual y mayor que la óptima obtenida en estado estacionario, es decir, 1.2921 min-1. Los valores de la velocidad espacio correspondientes a cada uno de los casos son 4/7, 1.2921 y 2.8744 min-1.

Mach

138

Caso 1 En la Figura (8.6) se muestra la comparación gráfica para la variación de la concentración de B que resulta mediante la solución de tanto el modelo no lineal como el linearizado para una velocidad espacio de 4/7 min-1 y un cambio paso en esta de 0.01 Solo se ha representado la respuesta con respecto a B para analizar su respuesta inversa

Figura 8.6 Concentración de B – Modelo Lineal y No Lineal (0.01) Se observa una buena aproximación del modelo lineal al no lineal cuando la variable de entrada es perturbada con un pequeño cambio paso. En el archivo solnolinealtime.m, codificado con Matlab, se incluye la determinación de la función de transferencia de la concentración de B con respecto a la velocidad espacio y la solución de dicho programa muestra que es:

C B ( s) =

Mach

0.5848(−0.3549s + 1) F − 1.1170s + 3.1472 F ( s) = ( s) 0.1858s 2 + 0.8627 s + 1 V s + 4.6429s + 5.3821 V 2

(8.16)

139

En la función de transferencia (8.16) se observa la inclusión de un zero positivo con un valor de 2.8175 (= 3.1472/1.1170) y una ganancia positiva (0.5848), lo que es característico de una respuesta inversa La Figura (8.7) muestra el resultado obtenido con un cambio paso mas grande de 0.1 observándose una diferencia mayor entre los modelos, con una ganancia significativamente inferior para el modelo no lineal

Figura 8.7 Concentración de B – Modelo Lineal y No Lineal (4/7) La aproximación entre los modelos se limita a las cercanías próximas al estado estacionario de la condición inicial. Las figuras (8.6) y (8.7) son características de un sistema con ganancia positiva y respuesta inversa (zero a la derecha del plano)

Caso 2 En la Figura (8.8) se muestra la comparación gráfica para la variación de la concentración de B que resulta mediante la solución de tanto el modelo no lineal como el linearizado para una velocidad espacio de 2.8744 min-1 y un cambio paso en esta de 0.1.

Mach

140

Figura 8.8 Concentración de B – Modelo Lineal y No Lineal (2.8744 – 0.1) La figura 8.8 es la respuesta característica de un sistema con una ganancia negativa y muestra que para un pequeño cambio paso en la variable de entrada, el modelo lineal es una buena aproximación al modelo no lineal. La función de transferencia de la concentración de B con respecto a la velocidad espacio obtenida con el mismo código en Matlab es:

C B ( s) =

− 1.1170s − 3.1472 F − 0.01362(0.3549s + 1) F (s) = ( s) 0.0384s 2 + 0.3945s + 1 V s + 10.2778s + 26.0508 V 2

(8.17)

En la función de transferencia (8.17) se observa la inclusión de un zero negativo con un valor de -2.8175 (= -3.1472/1.1170) y una ganancia negativa (-0.01362). La Figura (8.9) muestra el resultado obtenido con un cambio paso mayor de 2.8744 min-1 observándose una diferencia mayor entre los modelos, con una ganancia significativamente inferior para el modelo no lineal. Las figuras (8.8) y (8.9) son características de un sistema con ganancia negativa y zero a la izquierda del plano. Para un pequeño cambio paso en la variable de entrada el modelo lineal es una buena aproximación al modelo no lineal.

Mach

141

Figura 8.8 Concentración de B – Modelo Lineal y No Lineal (2.8744 – 2.8744) Es interesante anotar que todos los puntos de operación a la izquierda del óptimo en la curva de condiciones estacionarias tienen zeros en el plano derecho de representación numérica (respuesta inversa), mientras que todos aquellos a la derecha del óptimo no lo tienen.

Caso 3 En la Figura (8.9) se muestra la comparación gráfica para la variación de la concentración de B que resulta mediante la solución de tanto el modelo no lineal como el linearizado para el valor óptimo de la velocidad espacio, es decir, 1.2921 min-1 y para un pequeño cambio paso en esta de 0.04. La función de transferencia de la concentración de B con respecto a la velocidad espacio obtenida con el mismo código en Matlab es:

C B ( s) =

Mach

− 1.2660s F ( s) s + 6.5825s + 10.7217 V 2

(8.18)

142

En la función de transferencia (8.18) se observa un zero en el origen y se nota muy poco cambio con respecto a la condición en estado estacionario para un pequeño cambio paso en la variable de entrada. Esto es consistente con la nula ganancia que se observa en la función (8.18)

Figura 8.9 Concentración de B – Modelo Lineal y No Lineal (1.2921 – 0.04)

Modelo Lineal - Espacio de los Estados – Dominio Laplace El modelo en la forma del espacio de los estados expresado en términos de las variables desviación se puede transformar según Laplace de tal manera que la ecuación matricial es de la forma:

X ( s ) = ( sI − A) −1 Bu ( s )

(8.19)

Como usualmente, el vector de las variables de salida se expresa solamente en función de las variables de entrada, la matriz D es tal que anula el segundo término de la ecuación algebraica del modelo y la función de transferencia en forma matricial a partir del modelo del espacio de los estados se expresa como:

Mach

143

Siendo

Y ( s ) = C ( sI − A) −1 Bu ( s )

(8.20)

G ( s ) = C ( sI − A) −1 B

(8.21)

8.3 MATLAB: COMANDOS UTILIZADOS Comando fminbnd Para el análisis estacionario del sistema de reacciones de Van de Vusse se desarrolla el archivo vandevusse.m codificado con Matlab y que grafica las concentraciones de A y B en función de la velocidad espacio asignando valores entre 0 y 9 a intervalos de 0.2. En el mismo archivo se incluye la determinación de las condiciones óptimas para el componente B, para lo cual es necesario construir previamente otro archivo denominado conmax.m y utilizar el comando fminbnd cuya sintaxis es: [x,f] = fminbnd(@conmax,x1,x2) Para la obtención del valor máximo de la concentración de B se hace necesario plantear como función argumento en el comando fminbnd la función de signo opuesto porque con este se halla el mínimo de una función y que, por lo tanto, el máximo corresponde al mismo valor con signo opuesto. Se observa en la sintaxis que el nombre del archivo se precede de la letra @ y adicionalmente se requieren dos argumentos correspondientes al valor inicial y final del intervalo de búsqueda del valor óptimo.

Comando figure Para el desarrollo de la solución del modelo se aplican algunos comandos no utilizados como el figure y plotyy. Con el primero se enumera la ventana donde se construye una gráfica y puede invocarse posteriormente para colocar sobre ella una nueva gráfica y otra gráfica. La sintaxis del comando figure es: figure(h) Siendo h el número que se asigna a la ventana.

Mach

144

Comando plotyy El comando plotyy se aplica cuando se necesita construir dos gráficas sobre una misma ventana, especialmente cuando para un mismo intervalo de valores de la variable independiente los valores de las variables dependientes son muy diferentes. La sintaxis del comando plotyy es: plotyy(t,y1,t,y2) Los argumentos son la variable independiente y la dependiente correspondiente a cada una de las gráficas. Matlab muestra a la izquierda y a la derecha de la ventana de la figura las escalas correspondientes a cada una de ellas.

Comandos de conversión de funciones de transferencia Para la conversión del modelo expresado en la forma del espacio de los estados a las funciones de transferencias correspondientes, Matlab dispone del comando ss2tf cuya sintaxis es: [num,den] = ss2tf(A, B, C, D) En forma similar, Matlab dispone del comando tf2ss para la conversión de la función de transferencia a la forma del espacio de los estados y cuya sintaxis es: [A, B, C, D] = tf2ss(num,den)

8.4 MATLAB: SOLUCION DE LOS MODELOS El análisis estacionario se realiza con los archivos conmax.m y vandevusse.m. La respuesta paso del sistema del sistema de reacciones de Van de Vusse se desarrolla con los archivos nolinealtime.m y solnonolinealtime.m. En este último se incluye la solución linearizada del modelo utilizando su expresión en la forma del espacio de los estados. En la sección 8.5 se incluyen los códigos de los archivos mencionados y se deja como ejercicio para el estudiante la revisión de ellos y la construcción de un archivo unificado que muestre para cada uno de los casos los perfiles de A, C y D y otro archivo que muestre la comparación entre el modelo no lineal y el modelo lineal empleando la función de transferencia a partir de la forma del espacio de los estados Mach

145

8.5 MATLAB - PROGRAMAS CODIFICADOS Archivo conmax.m function f = conmax(x) global K1 K2 K3 K4 Cf f = - K1*(-(K1 + x) + sqrt((K1 + x)^2 + 4*K3*Cf*x))/(2*K3)/(K2 + x);

Archivo solconmax.m global K1 K2 K3 K4 Cf K1 = input('Valor de la constante K1 = '); K2 = input('Valor de la constante K2 = '); K3 = input('Valor de la constante K3 = '); Cf= input('Concentracion de A en el alimento = '); x1 = input('Valor inicial de busqueda = '); x2 = input('Valor final de busqueda = '); [x,f] = fminbnd(@conmax,0,2)

Archivo vandevusse.m % Reaccion de Van de Vusse-Condiciones Estacionarias % Cas, concentracion de A en estado estacionario % Cbs, concentracion de B en estado estacionario % x, Velocidad Espacio en el reactor K1 = input('Valor de la constante K1 = '); K2 = input('Valor de la constante K2 = '); K3 = input('Valor de la constante K3 = '); Cf= input('Concentracion de A en el alimento = '); x1 = input('Valor inicial del intervalo de busqueda = '); x2 = input('Valor final del intervalo de busqueda = '); for i = 1:46 x(i) = 0.2*(i - 1); Cas(i) = (-(K1 + x(i))/(2*K3)) + sqrt(((K1 + x(i))^2) + 4*K3*Cf*x(i))/(2*K3); Cbs(i) = K1*Cas(i)/(x(i) + K2); end plot(x,Cas) figure; plot(x,Cbs)

Archivo nolinealtime.m function dy = nolinealtime(t,y) Mach

146

global Fv Caf K1 K2 K3 dy = [Fv*(Caf - y(1)) - K1*y(1) - K3*y(1).^2; -Fv*y(2) + K1*y(1) - K2*y(2); Fv*y(3) + K2*y(2); -Fv*y(4) + 0.5*K3*y(1).^2];

Archivo solnolinealtime.m function f = solnolinealtime(t,y) global Fvs Fv DFv Caf Cas Cbs Ccs Cds K1 K2 K3 Inicio Rango K1 = input('Valor de la constante K1, 1/min = '); K2 = input('Valor de la constante K2, 1/min = '); K3 = input('Valor de la constante K3, l/mol-min = '); Caf = input('Concentracion de A en el alimento, mol/litro = '); Fvs = input('Valor de la velocidad espacio en estado estacionario, 1/min = '); DFv = input('Valor del cambio en la velocidad espacio, 1/min = '); Rango = input('Intervalo de tiempo = '); Cas = -(K1 + Fvs)/(2*K3) + (1/(2*K3))*sqrt((K1 + Fvs)^2 + 4*K3*Fvs*Caf) Cbs = K1*Cas/(K2 + Fvs) Ccs = K2*Cbs/Fvs Cds = 0.5*K3*Cas^2/Fvs Fv = Fvs + DFv; Inicio = [Cas Cbs Ccs Cds]; [t,y] = ode45('nolinealtime', Rango, Inicio); figure(1); plotyy(t,y(:,1),t,y(:,2)) figure(2); plotyy(t,y(:,3),t,y(:,4)) figure(3); plot(t,y(:,2),'k') %Espacio de los Estados % A = [-Fvs-K1-2*K3*Cas 0 0 0; K1 -Fvs-K2 0 0; 0 K2 -Fvs 0; K3*Cas 0 0 -Fvs]; B = [Caf-Cas; -Cbs; -Ccs; -Cds]; C = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]; D = [0; 0; 0; 0]; H = ss(A,B,C,D); H1 = DFv*H; [y,t] = step(H1,5); figure(3) hold on plot(t,y(:,2)+Cbs,'r') % Funcion de Transferencia entre A y B A1 = [A(1,1) A(1,2); A(2,1) A(2,2)]; B1 = [B(1,1); B(2,1)]; C1 = [C(1,1) C(1,2); C(2,1) C(2,2)]; D1 = [0; 0]; [num,den] = ss2tf(A1,B1,C1,D1) Mach