Saltar al contenido

Modelo de 2 tanques esféricos

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 h_1 y h_2 y dos variables manipuladas F_{i_1} y F_{i_2}.

Modelo tanques esfericos

La dinámica de los tanques puede ser modelada usando el principio de conservación de Masa:

\dfrac{dm}{dt}=\dfrac{dm_{in}}{dt}-\dfrac{dm_{out}}{dt}

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) \dfrac{dV_1}{dt}=F_{i_1}-F_{o_1}

(2) \dfrac{dV_2}{dt}=F_{i_2}+F_{o_1}-F_{o_2}

Donde V_i es el volumen acumulado en cada tanque y F_{i} y F_{o} son los caudales volumetricos de entrada y salida de cada tanque respectivamente. Dabiendo que dV=Adh:

(3) \dfrac{A_1(h_1)dh_1}{dt}=F_{i_1}-F_{o_1}

(4) \dfrac{A_2(h_2)dh_2}{dt}=F_{i_2}+F_{o_1}-F_{o_2}

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 h_1 y h_2, asi tomando en cuenta que dV=Adh procedemos a derivar el volumen de un tanque esferico con relacion a la altura, sin embargo para esto deberemos conocer la ecuación matemática del volumen de un casquete esférico:

V=\dfrac{\pi h^2}{3}(3R-h)

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:

Casquete esférico Tanques


Por lo tanto, los límites de integración serán de r - h a r.

Partiendo de la ecuación de la circunferencia con centro en el origen:

x=\sqrt{r^2+y^2}

Supongamos que el casquete esférico de altura h es formado por cilindros infinitos de alturas infinitesimales dy y radios x, donde x es variable para cada punto de h:

modelado casquete esférico


Sabemos que el volumen de un cilindro es:

V_c=A_bh=\pi r^2 h

Y en este caso:

V_c=\pi x^2 dy

La suma de estos cilindros infinitos de alturas infinitesimales forma el volumen del casquete esférico en los límites r - h y r. Entonces su volumen estará dado por:

V=\int_{r-h}^{r}\pi x^2 dy

V=\pi \int_{r-h}^{r} x^2 dy

Sustituyendo la ecuación del circulo tenemos


V=\pi \int_{r-h}^{r} (\sqrt{r^2+y^2})^2 dy

V=\pi \int_{r-h}^{r} (r^2+y^2) dy

Resolviendo

V=\pi\left[r^2y - \dfrac{y^3}{3} \right]_{r-h}^{r}

aplicamos limites

V=\pi\left[\left( r^2r - \dfrac{r^3}{3}\right) - \left( r^2(r-h) - \dfrac{(r-h)^3}{3} \right) \right]

V=\pi\left[\left( \dfrac{3r^3-r^3}{3}\right) - \left( \dfrac{3r^3-3hr^2-r^3-3hr^2+3rh^2-h^3}{3} \right) \right]

V=\dfrac{\pi h^2}{3}(3r-h)

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) A_1=\pi(2R_1h_1-h_1^2)

(6) A_2=\pi(2R_2h_2-h_2^2)

Los caudales de salida que pasan por una sección transversal y a través de una válvula pueden ser modelados como

(7) F_{o_1}=k_1a_1\sqrt{h_1}

(8) F_{o_2}=k_2a_2\sqrt{h_2}

De esa forma la ecuación diferencial que representa los tanques acoplados provenientes de la ecuación (3) y (4) vienen dado por:

(9) \pi(2R_1h_1-h_1^2)\dfrac{dh_1}{dt}=F_{i_1}-k_1a_1\sqrt{h_1}

(10) \pi(2R_2h_2-h_2^2)\dfrac{dh_2}{dt}=F_{i_2}+k_1a_1\sqrt{h_1}-k_2a_2\sqrt{h_2}

despejando el sistema tenemos:

(11) \dfrac{dh_1}{dt}=\dfrac{1}{\pi(2R_1h_1-h_1^2)}\left(F_{i_1}-k_1a_1\sqrt{h_1}\right)

(12) \dfrac{dh_2}{dt}=\dfrac{1}{\pi(2R_2h_2-h_2^2)}\left(F_{i_2}+k_1a_1\sqrt{h_1}-k_2a_2\sqrt{h_2}\right)

En el caso de los flujos de entrada, que provienen de las bombas, vamos a dejarlos indicados como F_{i_1} y F_{i_2}, sin embargo, por causa de usar bombas aquí, podemos representar esa entrada como:

F_{i_1}=K_pu_1(t) y F_{i_2}=K_pu_2(t) donde u_1(t) y u_2(t) son el voltaje aplicado a las bombas y el K_p es la constante de la bomba dada en unidades de caudal por voltaje, por ejemplo (L/v).

Estado Estacionario

Basta con cerar las derivadas y encontrar los estacionarios para ambas alturas.

Comenzamos con la altura 1

0=\dfrac{F_{i_1}^*}{\pi(2R_1h_1^*-(h_1^*)^2)}-\dfrac{k_1a_1\sqrt{h_1^*}}{ \pi(2R_1h_1^*-(h_1^*)^2)}

(13) h_1^*=\left(\dfrac{F_{i_1}^*}{k_1a_1}\right)^2

Ahora la altura 2

0=\dfrac{F_{i_2}^*}{\pi(2R_2h_2^*-(h_2^*)^2)}+\dfrac{k_1a_1\sqrt{h_1^*}}{\pi(2R_2h_2^*-(h_2^*)^2)}-\dfrac{k_2a_2\sqrt{h_2^*}}{\pi(2R_2h_2^*-(h_2^*)^2)}

h_2^*=\left(\dfrac{F_{i_2}^*}{k_2a_2}+\dfrac{F_{i_1}^*}{k_2a_2}\right)^2

(14) h_2^*=\left(\dfrac{F_{i_2}^*}{k_2a_2}\right)^2+\dfrac{2F_{i_1}^*F_{i_2}^*}{(k_2a_2)^2}+\left(\dfrac{F_{i_1}^*}{k_2a_1}\right)^2

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) G_1(h_1,h_2,F_{i_1},F_{i_2})=\dfrac{1}{\pi(2R_1h_1-h_1^2)}\left(F_{i_1}-k_1a_1\sqrt{h_1}\right)

(16) G_2(h_1,h_2,F_{i_1},F_{i_2})=\dfrac{1}{\pi(2R_2h_2-h_2^2)}\left(F_{i_2}+k_1a_1\sqrt{h_1}-k_2a_2\sqrt{h_2}\right)

\begin{bmatrix}\dot{h_1}\\ \\\dot{h_2}\end{bmatrix}=\begin{bmatrix}\dfrac{\partial G_1}{\partial h_1} & \dfrac{\partial G_1}{\partial h_2}\\ \\\dfrac{\partial G_2}{\partial h_1} & \dfrac{\partial G_2}{\partial h_2}\end{bmatrix}\begin{bmatrix}{h_1}\\ \\ \\h_2\end{bmatrix}+\begin{bmatrix}\dfrac{\partial G_1}{\partial F_{i_1}} & \dfrac{\partial G_2}{\partial F_{i_2}}\\ \\\dfrac{\partial G_2}{\partial F_{i_1}} & \dfrac{\partial G_2}{\partial F_{i_2}}\end{bmatrix}\begin{bmatrix}{F_{i_1}}\\ \\ \\F_{i_2}\end{bmatrix}

donde:

\dfrac{\partial G_1}{\partial h_1}=\dfrac{-F_{i_1}(2\pi R_1-2\pi h_1^*)}{(2\pi R_1 h_1^2-\pi (h_1^*)^2)^2}+ \dfrac{ \dfrac{-k_1a_1}{2\sqrt{h1^*}}(2\pi R_1h_1^*-\pi (h_1^*)^2)+k_1a_1\sqrt{h_1^*}(2\pi R_1-2\pi h_1^*) } { (2\pi R_1 h_1^2-\pi (h_1^*)^2)^2 }

\dfrac{\partial G_2}{\partial F_{i_2}}=0

\dfrac{\partial G_2}{\partial F_{i_1}}= \dfrac{ \dfrac{k_1a_1}{2\sqrt{h_1^*}} }{ (2\pi R_1 h_2^2-\pi (h_2^*)^2) }

    \begin{align*}\dfrac{\partial G_2}{\partial F_{i_2}}= & \dfrac{1}{(2\pi R_2 h_2^2-\pi (h_2^*)^2)^2} \\& \left(-F_{i_2}(2\pi R_2-2\pi h_2^*)+ \dfrac{-k_2a_2}{2\sqrt{h_2^*}}(2\pi R_2h_2^*-\pi (h_2^*)^2)+k_2a_2\sqrt{h_2^*}(2\pi R_2-2\pi h_2^*) \right.\\& \left. -k_1a_1\sqrt{h_1^*}(2\pi R_2-2\pi h_2^*)\right)\end{align*}

\dfrac{\partial G_1}{\partial F_{i_1}}=\dfrac{1}{\pi(2R_1h_1^*-h_1^2)}

\dfrac{\partial G_2}{\partial F_{i_2}} = \dfrac{\partial G_2}{\partial F_{i_1}} = 0

\dfrac{\partial G_2}{\partial F_{i_2}}=\dfrac{1}{\pi(2R_2h_2^*-h_2^2)}

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.

Tanques Esfericos
Linealización Tanques Esfericos

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í

>>>> BAJAR ARCHIVOS <<<<<

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.