Hola controleros y controleras sean bienvenidos una vez más a otra entrada de la pagina de Control Automático Educación, donde vamos a continuar con nuestro curso de Análisis de Sistemas. Hoy vamos a aprender como obtener el modelo matemático de un sistema de dos tanques esféricos en cascada a través de ecuaciones diferenciales, vamos a linealizar el sistema en un punto de operación para obtener sus representaciones en función de transferencia y en espacio de estados y finalmente simularemos la planta en MATLAB/SIMULINK.
Antes que nada te invito a VER EL CURSO GRATUITO DE ANALISIS DE SISTEMAS.
Modelo de 2 tanques esféricos
Continuemos con nuestros ejercicios de tanques de almacenamiento, vistos en nuestro curso de análisis de sistema.
A continuación se muestra un sistema de 2 tanques esféricos, el cual conforma un sistema MIMO (multi-Input Multi-Output) con dos variables controladas y
La dinámica de los tanques puede ser modelada usando el principio de conservación de Masa:
Se asume que el fluido dentro de los tanques es incompresible y que el fluido posee una densidad constante, de esa forma las ecuaciones diferenciales de los dos tanques son formulados como la variación del volumen en ambos tanques
(1)
(2)
Donde
(3)
(4)
Además, el área de la superficie libre de fluido en los tanques 1 y 2 podrían expresarse en función de las de alturas de los líquidos
Si no sabes obtener esta ecuación, te voy a dejar el procedimiento para el cálculo del casquete esférico a través del calculo integral:
Volumen de un Tanque Esférico (Casquete Esférico)
A continuación mostramos el procedimiento para el cálculo del volumen de una cascara esférica.
El casquete esférico es el sólido generado a partir de una esfera cuando está seccionada por un plano:
Por lo tanto, los límites de integración serán de
Partiendo de la ecuación de la circunferencia con centro en el origen:
Supongamos que el casquete esférico de altura
Sabemos que el volumen de un cilindro es:
Y en este caso:
La suma de estos cilindros infinitos de alturas infinitesimales forma el volumen del casquete esférico en los límites
Sustituyendo la ecuación del circulo tenemos
Resolviendo
aplicamos limites
Modelado de dos tanques de Nivel en Cascada
Control PID con Raspberry Pi Pico
Modelado de un Reactor CSTR de Van de Vusse
Continuación del Modelo de los 2 Tanques Esféricos
Volviendo al ecuacionamiento dinámico de los dos tanques esféricos acoplados, podremos entender como se comporta el llenado y vaciado de un tanque esférico:
(5)
(6)
Los caudales de salida que pasan por una sección transversal y a través de una válvula pueden ser modelados como
(7)
(8)
De esa forma la ecuación diferencial que representa los tanques acoplados provenientes de la ecuación (3) y (4) vienen dado por:
(9)
(10)
despejando el sistema tenemos:
(11)
(12)
En el caso de los flujos de entrada, que provienen de las bombas, vamos a dejarlos indicados como
Estado Estacionario
Basta con cerar las derivadas y encontrar los estacionarios para ambas alturas.
Comenzamos con la altura 1
(13)
Ahora la altura 2
(14)
Espacio de Estados
A partir de las ecuaciones (11) y (12) es posible obtener la representacón en espacio de estados del sistema linealizado usando la Jacobiana y empleando los puntos de equilibrio de las ecuaciones (13) y (14). Para eso reescribimos las ecuaciones (11) y (12) como:
(15)
(16)
donde:
Comparación Modelo No Lineal y Modelo Lineal
La siguiente grafica muestra la comparación entre los dos modelos, en un punto muy cercano al punto de equilibrio donde fue linealizado el sistema, con una pequeña desviación de la entrada del tanque 1 de magnitud 0.1. La gráfica nos sirve para entender el comportamiento dinámico de este sistema, como por ejemplo determinar su ganancia estática o saber el tiempo de llenado de un tanque.
Simulación de 2 Tanques Esféricos en Matlab
A continuación esta disponible el código que te permitirá modelar y simular los tanques esfericos acoplados que vimos en esta entrada. Este código tiene la capacidad de obtener el perfil dinámico de los tanques, obtener el modelo lineal en un punto de operación y finalmente hacer una comparación entre el modelo lineal y no lineal.
Para tener acceso al código basta con que compartas esta información para ayudar que este sitio web siga llegando a mas personas y continue compartiendo información de interés.
[sociallocker id=»948″]
Los dos Archivos están disponibles aquí
Código Principal
%% Sergio Castaño % Rio de Janeiro - 2019 % Dos tanques Esféricos proceso MIMO %% clc close all clear all FS=14; %Tamanho da Fonte das Figuras nit=300; Ts=1; tsim = 0:Ts:(nit-1)*Ts; %Tiempo de Simulación %Parametros R1 = 4; %Radio 1 [m] R2 = 6; %Radio 2 [m] k1 = 3.2; %Coeficiente de la Válvula Tanque 1 [m^2/s] k2 = 6.8; %Coeficiente de la Válvula Tanque 1 [m^2/s] a1 = 0.4; a2 = 0.5; par.R1 = R1; %Radio 1 [m] par.R2 = R2; %Radio 2 [m] par.k1 = k1; %Coeficiente de la Válvula Tanque 1 [m^2/s] par.k2 = k2; %Coeficiente de la Válvula Tanque 1 [m^2/s] par.a1 = a1; par.a2 = a2; %Entrada F1=1.75; F2=1; u.F1=F1; u.F2=F2; %Condição inicial x0=[0.1 0.1]; %Simulação do Sistema [t,x] = ode45(@(t,x)TwoSphericalT(t,x,u,par), tsim , x0, u); %Grafica do Sistema figure plot(t,x(:,1),[t(1) t(end)],[F1 F1],'linewidth',2) title('Altura h_1'); ylabel('h_1'); xlabel('tempo (s)'); figure plot(t,x(:,2),[t(1) t(end)],[F2 F2],'linewidth',2) title('Altura h_2'); ylabel('h_2'); xlabel('tempo (s)'); %Grafico estacionario F1x=0:0.1:3; F2x=1; H1x=(F1x/(k1*a1)).^2; H2x=(F2x/(k2*a2)).^2+(F1x/(k2*a2)).^2+((2*F1x*F2x)/(a2*k2)^2); figure subplot(2,1,1); plot(F1x,H1x),grid title('Variación de nivel h1'); ylabel('Altura (h)'); xlabel('Variação F1 com F2=1'); subplot(2,1,2); plot(F1x,H2x),grid title('Variación de nivel h2'); ylabel('Variação F1 com F2=1'); xlabel('Abertura de Valvula a1'); %Encuentra el estado Estacionario X = fsolve(@(x)TwoSphericalT(t,x,u,par),x0); %X=[h1 h2] %Establesco los Estados en el Estado Estacionario h1s = X(1); h2s = X(2); Eq=[h1s h2s F1 F2]; %% Linealizacion del sistema por Jacobiana (Libreria SYMBOLICA) syms x1 x2 u1 u2 A1 = pi*(2*R1*x1-x1^2); A2 = pi*(2*R2*x2-x2^2); % Funciones para realizar la jacobiana fx1 = 1/A1 * (u1 -a1*k1*sqrt(x1)); fx2 = 1/A2 * (u2 + a1*k1*sqrt(x1) - a2*k2*sqrt(x2)); %Conforma vectores fx = [fx1;fx2]; x = [x1;x2]; uj = [u1;u2]; %Linealiza por medio de la Jacobiana A = jacobian(fx,x); B = jacobian(fx,uj); C = eye(2); %Matriz de salida (Cb) %Reemplazo los puntos de equilibrio en la Jacobiana display('Representación en Espacio de Estado Caso 1') Ad = double(subs(A,{x1,x2,u1,u2},{Eq})) Bd = double(subs(B,{x1,x2,u1,u2},{Eq})) % Linealización usando Jacobiana %Cálculo da matriz A a11 = (-F1*(2*pi*R1-2*pi*h1s))/(2*pi*R1*h1s-pi*h1s^2)^2 ... + ( (-k1*a1)/(2*sqrt(h1s))*(2*pi*R1*h1s-pi*h1s^2)+k1*a1*sqrt(h1s)*(2*pi*R1-2*pi*h1s) ) / (2*pi*R1*h1s-pi*h1s^2)^2; a12 = 0; a21 = ( (k1*a1)/(2*sqrt(h1s)) ) / (2*pi*R2*h2s-pi*h2s^2); a22 = (-F2*(2*pi*R2-2*pi*h2s))/(2*pi*R2*h2s-pi*h2s^2)^2 ... + ( (-k2*a2)/(2*sqrt(h2s))*(2*pi*R2*h2s-pi*h2s^2)+k2*a2*sqrt(h2s)*(2*pi*R2-2*pi*h2s) ) / (2*pi*R2*h2s-pi*h2s^2)^2 ... +( -k1*a1*sqrt(h1s)*(2*pi*R2-2*pi*h2s) ) / (2*pi*R2*h2s-pi*h2s^2)^2; A = [a11 a12 a21 a22]; %Cálculo da matriz B B = [1/(pi*(2*R1*h1s-h1s^2)) 0 0 1/(pi*(2*R2*h2s-h2s^2))]; %Cálculo da matriz C C = [1 0 0 1]; %Cálculo da matriz D D = [0 0 0 0]; % %Determino la Funcion de Transferencia para la entrada 1 Cb % display('Representación en Función de Transferencia para Cb') % [num1,den1]=ss2tf(A,B,C,D,1); % ft1=tf(num1,den1) H = ss(A,B,C,D); [z,P,k] = zpkdata(H); K = zpk(z,P,k) %Comprobar Modelo Lineal VS Modelo NO Lineal du=0.1; u.F1=u.F1+du; [t,x] = ode45(@(t,x)TwoSphericalT(t,x,u,par), tsim , X, u); in=zeros(length(tsim),2); in(:,1)=du*ones(length(tsim),1); in(:,2)=0*ones(length(tsim),1); [ylin]=lsim(H,in,tsim); %Grafica figure plot(t,x(:,1),'-b',tsim,ylin(:,1)+X(1),'--r','LineWidth',2) ylabel('h_1','fontsize',FS); titulo='Comparação modelo Linear vs Não Linear (h_1) para u= ','fontsize',FS; titulo=strcat(titulo,num2str(F1),'h^{-1}'); title(titulo); xlabel('tempo (s)','fontsize',FS); legend('Não Linear','Linear','Location','SouthEast') set(gca,'fontsize',FS); figure plot(t,x(:,2),'-b',tsim,ylin(:,2)+X(2),'--r','LineWidth',2) ylabel('h_2','fontsize',FS); titulo="Comparação modelo Linear vs Não Linear (h_2) para u= "; titulo=strcat(titulo,num2str(F2),"h^{-1}"); title(titulo); xlabel('tempo (s)','fontsize',FS); legend('Não Linear','Linear','Location','SouthEast') set(gca,'fontsize',FS);
Modelo de los 2 Tanques Esféricos
function dhdt = TwoSphericalT(t,x,u,par) %% Datos % Entradas Fi1 = u.F1; %h^-1 (Flujo) Fi2 = u.F2; %h^-1 (Flujo) %Parametros R1 = par.R1; %Radio 1 [m] R2 = par.R2; %Radio 2 [m] k1 = par.k1; %Coeficiente de la Válvula Tanque 1 [m^2/s] k2 = par.k2; %Coeficiente de la Válvula Tanque 1 [m^2/s] a1 = par.a1; a2 = par.a2; % % Notación para las Variables de Estado % h1 = x(1); %height 1 [m] h2 = x(2); %height 2 [m] % Areas A1 = pi*(2*R1*h1-h1^2); A2 = pi*(2*R2*h2-h2^2); % % Ecuaciones Diferenciales del Sistema % dhdt = zeros(2,1); dhdt(1) = 1/A1 * (Fi1 -a1*k1*sqrt(h1)); dhdt(2) = 1/A2 * (Fi2 + a1*k1*sqrt(h1) - a2*k2*sqrt(h2));
[/sociallocker]
Bibliografia
- Tavakolpour-Saleh, A. ; Setoodeh, A. ; Ansari, E. (2016), ‘Iterative Learning Control of Two Coupled Nonlinear Spherical Tanks ‘, World Academy of Science, Engineering and Technology, Open Science Index 119, International Journal of Mechanical and Mechatronics Engineering, 10(11), 1862 – 1869.
Eso es todo por la entrada del dia de hoy, espero les haya gustado y hayan aprendido algo nuevo. Si te ha servido el contenido de esta entrada, de los videos y los códigos de implementación y deseas apoyar mi trabajo invitandome a un café super barato, puedes hacerlo en el siguiente link:
👉 Invitar a Sergio a un Café ☕️
Que esten muy bien, nos vemos en la siguiente entrada.
Mi nombre es Sergio Andres Castaño Giraldo, y en este sitio web voy a compartir una de las cosas que mas me gusta en la vida y es sobre la Ingeniería de Control y Automatización. El sitio web estará en constante crecimiento, voy a ir publicando material sobre el asunto desde temas básicos hasta temas un poco más complejos. Suscríbete al sitio web, dale me gusta a la página en Facebook y únete al canal de youtube. Espero de corazón que la información que comparto en este sitio, te pueda ser de utilidad. Y nuevamente te doy las gracias y la bienvenida a control automático educación.