Modelado y Linealización en Sistemas de Tanques

En esta entrada vamos a resolver paso a paso un ejercicio completo de modelado, linealización y análisis de estabilidad de un sistema de dos tanques interactuantes, donde uno de ellos es un tanque cónico. Este ejercicio integra prácticamente todos los conceptos centrales del curso de Control de Procesos: balance de masas, linealización por matriz Jacobiana, espacio de estados, función de transferencia, criterio de Routh-Hurwitz, aproximación de Padé, error en estado estacionario y principio del modelo interno.

Si lo deseas, puedes ver nuestro curso de Análisis de Sistemas para aprender a modelar otros tipos de procesos, y suscribirte al canal de YouTube para seguir aprendiendo más sobre el modelado y control de procesos industriales.

Si te interesa modelar otros sistemas hidráulicos, puedes darle un vistazo a las entradas relacionadas:

• Modelo Matemático de un Tanque de Nivel (un único tanque cilíndrico).
• Modelado matemático de Dos Tanques Interconectados en Serie.
• Modelo de 2 Tanques Esféricos.

Video en Español

Enunciado del Sistema

Consideremos el siguiente sistema conformado por dos tanques interactuantes en paralelo: un tanque cónico seguido de un tanque cilíndrico, comunicados por la parte inferior.

Donde:

• : Radio máximo del tanque cónico (m).
• : Altura total del tanque cónico (m).
• : Caudal de entrada al tanque cónico (con presión constante).
• : Caudal de interacción entre el tanque cónico y el cilíndrico.
• : Caudal de salida del sistema.
• : Altura del líquido en el tanque cónico (m).
• : Altura del líquido en el tanque cilíndrico (m).
• : Abertura de la válvula lineal de entrada (entre 0 y 1).
• : Área transversal del tanque cilíndrico.

Los caudales obedecen a las siguientes ecuaciones:

El caudal de entrada corresponde a una válvula con comportamiento estático lineal (abertura entre 0 y 1), mientras que y provienen del principio de Bernoulli (consulta el modelado del tanque de nivel donde se demuestra cómo se llega a esa expresión).

Los parámetros del modelo son:

Recordemos que el volumen de un cono viene dado por:

Modelado Matemático: Balance de Masa

Aplicamos el balance de masas a cada tanque. Lo que entra menos lo que sale es igual a la variación del volumen acumulado.

Tanque Cónico

Sustituyendo el volumen del cono:

Aquí aparece una complicación adicional: en el tanque cónico, el radio del líquido también cambia con la altura. Para eliminar esa dependencia usamos la proporción lineal del triángulo:

Sustituyendo en el volumen:

Para simplificar la notación definimos:

Aplicando la regla de la cadena :

Despejando la derivada llegamos a la primera ecuación diferencial no lineal:

Tanque Cilíndrico

El balance es más sencillo porque el área es constante:

A_2 \frac{dh_2}{dt} = f_1 - f_2 = k_1\sqrt{2g}\sqrt{h_1 - h_2} - k_2\sqrt{2g}\sqrt{h_2}

\boxed{\frac{dh_2}{dt} = \frac{k_1\sqrt{2g}\sqrt{h_1 - h_2} - k_2\sqrt{2g}\sqrt{h_2}}{A_2}}

Ya tenemos las dos EDOs no lineales que describen el comportamiento del sistema. Note la no linealidad: las alturas aparecen elevadas al cuadrado (en el denominador) y dentro de raíces cuadradas.

Curva Estática del Proceso

Antes de linealizar, conviene visualizar el comportamiento estático del sistema: cómo varía la altura del tanque cilíndrico (nuestra salida de interés) en función de la abertura de la válvula de entrada cuando el sistema alcanza el equilibrio.

En equilibrio, las derivadas se hacen cero, lo que conduce a un sistema lineal en términos de y que se puede resolver explícitamente. La curva estática resultante es:

Vemos que la respuesta estática es claramente no lineal (forma parabólica). Esto significa que la ganancia del proceso varía según el punto de operación, lo cual es típico en sistemas hidráulicos donde el caudal depende de .

Linealización por Matriz Jacobiana

Como tenemos un sistema multivariable con dos ecuaciones diferenciales acopladas, la forma más limpia de linealizar es usando la matriz Jacobiana, que es simplemente la extensión matricial del desarrollo de Taylor de primer orden para sistemas con varias variables.

Definimos las funciones:

Y trabajamos en variables de desviación:

El sistema linealizado en formato matricial queda:

Como solo nos interesa la altura del tanque cilíndrico:

Cálculo de las Derivadas Parciales

El término más laborioso es , ya que aparece en varios sitios. Aplicando regla del producto sobre :

Combinando los dos últimos términos sobre denominador común :

Para aplicamos regla de la cadena con y :

Las derivadas para son más directas:

Y para la entrada:

Punto de Equilibrio

Igualando las derivadas a cero llegamos al siguiente sistema lineal en y :

Que se resuelve de forma directa con MATLAB:

Ao = [1, -1;
      1, -((k2/k1)^2+1)];
bo = [(a0*k0/(k1*sqrt(2*g)))^2; 0];

Ho  = Ao\bo;
h10 = Ho(1);
h20 = Ho(2);

Para nuestro punto de operación , los valores de equilibrio son exactamente:

Función de Transferencia

Una vez obtenidas las matrices , y del sistema linealizado en espacio de estados, podemos calcular la función de transferencia aplicando:

En MATLAB usamos directamente ss2tf, obteniendo una función de transferencia de segundo orden:

Con polos reales en y , lo cual indica un sistema sobreamortiguado con dos constantes de tiempo dominantes (aproximadamente s y s). La ganancia estática resulta ser .

Comparación: Modelo Lineal vs No Lineal

Implementamos ambos modelos en Simulink y le aplicamos un pequeño escalón de a partir del punto de operación. La gráfica resultante muestra:

Ambas respuestas son bastante parecidas para este pequeño escalón, lo que valida la linealización. Existe una leve diferencia en estado estacionario, lo cual es esperado: el modelo lineal es solamente una aproximación local del sistema real alrededor del punto de operación. Si nos alejamos demasiado de ese punto, la diferencia crecerá.

Es importante notar que la función de transferencia trabaja desde el origen , mientras que el sistema no lineal arranca desde el punto de equilibrio m. Por eso, en Simulink se le suma este offset al sistema lineal para que ambas respuestas se puedan comparar gráficamente.

Comparación: Lazo Abierto vs Lazo Cerrado

Otro análisis interesante consiste en comparar la respuesta del sistema en lazo abierto (planta sola) versus lazo cerrado (con realimentación unitaria, equivalente a un control proporcional con ).

Vemos dos efectos claros del feedback:

• La constante de tiempo cambia: el lazo cerrado se estabiliza más rápido que el lazo abierto.
• La ganancia cambia: el valor final del lazo cerrado es mucho menor que el de lazo abierto (existe un error en estado estacionario).

Esta es la magia del feedback: la realimentación acelera la dinámica del sistema, a costa de modificar la ganancia estacionaria. Esto motiva el siguiente análisis.

El Reto del Tiempo Muerto

En sistemas reales casi siempre existe un tiempo muerto o retardo de transporte (por ejemplo, retardo en la comunicación, transporte de fluido, etc.). Vamos a incluir un retardo de 1 segundo y analizar cómo afecta la estabilidad bajo un control proporcional .

El problema es que un retardo puro no es una función racional, por lo que no se puede aplicar Routh-Hurwitz directamente. Para sortear esto usamos la aproximación de Padé de primer orden:

En MATLAB:

Gd = G;
Gd.iodelay = 1;            % Función de transferencia con retardo puro
[num, den] = pade(1, 1);
Gp = G * tf(num, den);     % Función de transferencia con retardo aproximado

Multiplicando por la aproximación, el sistema pasa de orden 2 a orden 3:

Margen de Ganancia con Routh-Hurwitz

La ecuación característica del sistema en lazo cerrado con un control proporcional es:

Que al desarrollar resulta en:

Construyendo la tabla de Routh:

Donde:

Para que el sistema sea estable, todos los términos de la columna pivote deben ser positivos. Esto nos da dos condiciones simultáneas:

Por lo tanto, el rango de estabilidad estimado por Routh-Hurwitz con la aproximación de Padé es:

Verificación con el comando margin

Como MATLAB permite trabajar con el retardo puro mediante la propiedad iodelay, podemos verificar el margen de ganancia real (sin aproximar) usando el comando margin:

[Gm, Pm] = margin(Gd);

El resultado es . Es decir, la aproximación de Padé(1,1) sobreestimó el margen en aproximadamente un 9% (4,62 vs 4,25). Esta es una lección pedagógica importante: las aproximaciones racionales del retardo (Padé de primer orden) son útiles para análisis algebraicos como Routh-Hurwitz, pero no son exactas. Si el análisis es crítico, conviene complementar con el comando margin sobre el sistema con retardo puro, o usar Padé de orden superior.

Análisis No Lineal

Aplicando al sistema no lineal con retardo (en Simulink), observamos que el sistema aún sigue siendo estable, e incluso podemos subir hasta sin que se inestabilice. Esto ocurre porque:

1. El análisis de margen de ganancia es solo válido localmente alrededor del punto de operación.
2. Las no linealidades del sistema real pueden ampliar (o reducir) el margen efectivo, dependiendo del proceso.

Este es un buen ejemplo de por qué los análisis con función de transferencia son guías de diseño útiles, pero siempre deben validarse en la planta real (o en una simulación no lineal).

Error en Estado Estacionario y Principio del Modelo Interno

Con un estable, la planta lineal con retardo es de tipo 0 (no tiene integradores en el lazo). Por lo tanto, ante una entrada escalón presentará un error en estado estacionario finito, y ante una rampa el error será infinito. Si quieres profundizar en este tema, te dejo la entrada de Error en Estado Estacionario.

Para llevar el error a cero ante una entrada escalón, aplicamos el Principio del Modelo Interno: el lazo de control debe contener una copia del modelo de la señal que se quiere seguir (o rechazar). Como un escalón se transforma como por Laplace, basta con añadir un integrador al lazo (es decir, pasar de un controlador puramente proporcional a uno PI).

Hay que tener en cuenta que añadir un integrador agrega dinámica adicional al sistema, por lo que la ganancia debe re-sintonizarse para mantener la estabilidad. Esto lo profundizaremos en próximas entradas cuando veamos el cálculo de controladores PI.

Análisis de Robustez: ¿Y si el retardo crece?

¿Qué pasa si el retardo aumenta de 1 a 10 segundos? Si dejamos el mismo controlador proporcional sintonizado para retardo de 1s, el sistema se vuelve completamente inestable.

¿Es posible «arreglar» el problema simplemente bajando ? Técnicamente sí, pero a costa de hacer el sistema extremadamente lento. Y aún así no es la mejor estrategia.

Cuando el retardo es dominante (mayor que la constante de tiempo del proceso), las estructuras de control clásicas (P, PI, PID) tienen un margen de fase muy pequeño y el desempeño se degrada drásticamente. La estrategia adecuada es usar predictores, como el Predictor de Smith, que separa el problema de control de la planta sin retardo del manejo del retardo, recuperando el desempeño nominal.

Este tema lo trataremos en detalle en una próxima entrada.

Código MATLAB y Simulink

A continuación les dejo los archivos completos en MATLAB y Simulink para que puedan reproducir todos los resultados del video. El paquete contiene el script taller_1_tanques.m y el modelo de Simulink taller1.slx con los bloques no lineal y lineal del sistema. Recuerden ejecutarlos en la misma carpeta.

Descargar Archivos en MATLAB y Simulink

%% Taller 1 de Sistemas de Control Continuo
% Linealización de un sistema de tanques interactuantes (cónico + cilíndrico)
% Prof. Sergio A. Castaño Giraldo
% 2025 - 2

clc
clear all
close all

%% Parámetros del modelo
A2 = 0.5;  %m^2  - Área del tanque cilíndrico
k0 = 0.25; %m^3/s - Constante de la válvula de entrada
k1 = 0.05; %m^2  - Constante de interacción
k2 = 0.05; %m^2  - Constante de la válvula de salida
g  = 10;   %m/s^2
H  = 1.5;  %m   - Altura máxima del tanque cónico
R  = 0.4;  %m   - Radio máximo del tanque cónico

% Constante geométrica del tanque cónico
alpha = pi * (R/H)^2;

%% Curva Estática variando a0
a0_vec = 0.0:0.05:1;
h2_vec = 0*a0_vec;
for i = 1:length(a0_vec)
    Ao = [1, -1;
          1, -((k2/k1)^2 + 1)];
    bo = [(a0_vec(i)*k0/(k1*sqrt(2*g)))^2; 0];
    
    Ho = Ao\bo;
    h2_vec(i) = Ho(2);
end
figure
plot(a0_vec, h2_vec, 'LineWidth', 2.5), grid on
xlabel('Abertura de válvula: a_0')
ylabel('Altura del tanque 2: h_2 [m]')
title('Curva estática del proceso')

%% Punto de operación
a0       = 0.6;       % Abertura de válvula en el punto de equilibrio
Delta_a0 = 0.05;      % Pequeño escalón para validar la linealización

%% Punto de Equilibrio
Ao = [1, -1;
      1, -((k2/k1)^2 + 1)];
bo = [(a0*k0/(k1*sqrt(2*g)))^2; 0];

Ho  = Ao\bo
h10 = Ho(1);
h20 = Ho(2);

fprintf('Punto de equilibrio:\n');
fprintf('  h10 = %.4f m\n', h10);
fprintf('  h20 = %.4f m\n', h20);

%% Linealización por Matriz Jacobiana
A = [ -2*a0*k0/(alpha*h10^3) + k1*sqrt(2*g)*(4*(h10-h20)-h10)/(2*alpha*h10^3*sqrt(h10-h20)), ...
       k1*sqrt(2*g)/(2*alpha*h10^2*sqrt(h10-h20));
       k1*sqrt(2*g)/(2*A2*sqrt(h10-h20)), ...
      -(k1*sqrt(2*g)/(2*A2*sqrt(h10-h20)) + k2*sqrt(2*g)/(2*A2*sqrt(h20))) ];

B = [k0/(alpha*h10^2); 0];

C = [0 1];

D = 0;

%% Función de Transferencia
[NUM, DEN] = ss2tf(A, B, C, D);

G  = tf(NUM, DEN)        % FT lineal sin retardo

% FT con retardo puro (para usar margin)
Gd = G;
Gd.iodelay = 1

% FT con aproximación de Padé de primer orden (para Routh-Hurwitz)
[num_pade, den_pade] = pade(1, 1);
Gp = G * tf(num_pade, den_pade)

%% Análisis de estabilidad
% Margen de ganancia con retardo puro
[Gm, Pm] = margin(Gd);
fprintf('\nMargen de ganancia (retardo puro) = %.4f\n', Gm);
fprintf('Margen de fase = %.4f grados\n', Pm);

% Mapa de polos y ceros con la aproximación de Padé
figure
pzmap(Gp)
title('Mapa de polos y ceros con aproximación de Padé(1,1)')
grid on

%% Comparación Modelo Lineal vs No Lineal
% Requiere el archivo taller1.slx en la misma carpeta
sim('taller1.slx')
figure
plot(t, ynl, '-r', 'LineWidth', 2.5), hold on
plot(t, yl,  '--k', 'LineWidth', 2.5), grid on
legend('No Lineal', 'Lineal', 'Location', 'best')
xlabel('Tiempo (s)')
ylabel('Altura h_2 (m)')
title('Dinámica de la altura del Tanque 2: Lineal vs No Lineal')

Conclusiones

En esta entrada hemos cubierto un flujo completo de análisis de un sistema de control:

1. Modelado no lineal de un sistema de tanques interactuantes con un tanque cónico, usando balance de masas y la proporción del triángulo para tratar la geometría variable.
2. Linealización por matriz Jacobiana, que es la forma natural de extender Taylor a sistemas multivariables.
3. Comparación lineal vs no lineal, validando que la linealización captura bien la dinámica local.
4. Efecto del feedback sobre la constante de tiempo y la ganancia.
5. Análisis de estabilidad con tiempo muerto vía Routh-Hurwitz y aproximación de Padé, comparando con el comando margin.
6. Error en estado estacionario y aplicación del Principio del Modelo Interno.
7. Robustez ante incrementos del retardo, motivando el uso de estructuras avanzadas como el Predictor de Smith.

Este ejercicio integra conceptos de prácticamente todo un curso de Control de Procesos, y muestra cómo las herramientas teóricas (Taylor, Jacobiano, Routh, Padé) se conectan con MATLAB y Simulink para resolver problemas reales.

Eso es todo por la entrada 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 invitándome a un café super barato, puedes hacerlo en el siguiente link:

👉 Invitar a Sergio a un Café ☕️

Que estén muy bien, nos vemos en la siguiente entrada.

Volver al Curso de Control de Procesos

Sergio A. Castaño Giraldo

¡Hola! Soy Sergio Andrés Castaño Giraldo. Bienvenidos a Control Automático Educación, un espacio dedicado a una de mis grandes pasiones: la Ingeniería de Control y la Automatización.

Mi objetivo aquí es tender un puente entre la teoría académica y la práctica profesional. Como profesor e investigador, he tenido la oportunidad de trabajar en proyectos complejos —desde sistemas de compresión en plataformas offshore hasta modelos de control avanzado— y todo ese aprendizaje lo vuelco en este sitio.

Aquí encontrarás material que crece constantemente: desde fundamentos de instrumentación hasta estrategias de control avanzado y programación profesional. Mi meta es que la información que comparto te sea de utilidad real en tu carrera o formación.

Te invito a que te suscribas, te unas a nuestro canal de YouTube y formemos una comunidad de aprendizaje continuo. ¡Gracias por estar aquí!