Sistemas de Control
ESCUELA SUPERIOR POLITÉCNICA DEL LITORAL
Facultad de Ingeniería en Electricidad y Computación
Sistemas de Control
Tema:
Avance 1 del proyecto
Presentado por:
Aguirre Vacas José David
Vásquez Isaac
Paralelo: 101
Profesor:
Ing. Adriana Aguirre
Término Académico:
II PAO – 2024
Sistemas de Control
Ecuaciones del sistema
Al identificar el caso de estudio como un problema que involucra transferencia de calor
y cambios en las concentraciones, podemos concluir que será pertinente el uso de
ecuaciones fundamentales de conversación de masa y energía.
Conservación de masa:
Conservación de energía:
Ley de proporción de Arrhenius:
Sistemas de Control
Parte no lineal
Subsistema interiormente
Sistemas de Control
Código:
%% PARTE 1A: Modelo no lineal
%Parámetros del Proyecto (Puntos de operación)
Q = 1; % Caudal de la corriente caliente (m3/s)
V = 1; % Caudal de la corriente fría (m3/s)
deltaH = -3000; % Cambio de entalpía en la corriente caliente (J/kg)
UA = 300; % Coeficiente de transferencia de calor multiplicado por el
area del reactor (W/K)
CAiop = 15; % Concentracion del reactivo A ingresando al reactor
(mol/L)
Caini = 15;
Tinop = 310; % Temperatura del reactivo A ingresando al reactor (K)
Tini = 310;
Tcop = 285; % Temperatura del refrigerante circulando por la camisa
(K)
delta_CAi = 0.02; % Cambio en la concentración del reactivo A
ingresando al reactor (10%)
delta_Tin = -0.5; % Cambio en la temperatura del reactivo A
ingresando al reactor (K)
delta_Tc = 0.01; % Cambio porcentual en la temperatura del
refrigerante circulando por la camisa (K)
% Definir la variable 'ruido' para la concentración de A (valor
constante)
ruido = 1.00E-04; % Perturbación constante en la concentración de A
% Definición de constantes del sistema
R = 1.985875; % Constante de los gases ideales (kcal/(kmol·K))
E = 11843; % Energía de activación (kcal/kmol)
k0 = 34930800; % Factor pre-exponencial (1/h)
rhoCp = 500; % Producto densidad-capacidad calorífica
(kcal/(K·m^3))
Tc = 292; % Temperatura del refrigerante (K)
dH = -5960; % Calor de reacción (kcal/kmol)
T0 = 311.267; % Temperatura inicial de A (K)
CA0 = 8.5695; % Concentración de A inicial (kmol/m^3)
% Declaración de variables simbólicas
syms CA T
CAi = CAiop + ruido; % Agregar el ruido a la concentración de A
Tin = Tinop; % Temperatura de entrada del reactivo A (K)
Tc = Tcop; % Temperatura del refrigerante (K)
% Ecuaciones diferenciales para la concentración (CA) y temperatura
(T) en el reactor
dCAdt = (Q/V)*(CAi - CA) - k0*exp(-E/(R*T)) * CA;
dTdt = (Q/V)*(Tin - T) + (UA/(rhoCp*V))*(Tc - T) +
(deltaH/rhoCp)*k0*exp(-E/(R*T)) * CA;
% Definir las ecuaciones diferenciales como funciones anónimas
f1 = @(t, y) (Q/V)*(CAi - y(1)) - k0*exp(-E/(R*y(2))) * y(1);
f2 = @(t, y) (Q/V)*(Tin - y(2)) + (UA/(rhoCp*V))*(Tc - y(2)) +
(deltaH/rhoCp)*k0*exp(-E/(R*y(2))) * y(1);
% Resolver las ecuaciones con ode45
y0 = [CA0, T0]; % Condiciones iniciales con el ruido agregado
tspan = [0 100]; % Tiempo de simulación
Sistemas de Control
[t, y] = ode45(@(t, y) [f1(t, y); f2(t, y)], tspan, y0);
% Salidas
CA = y(:, 1); % Concentración de A en el reactor
T = y(:, 2); % Temperatura del reactor
% Generación de las señales dependientes del tiempo
% Ahora que 't' está definido, calculamos las señales adicionales
CA_in = CAiop + ruido + delta_CAi * t; % Concentración de A
ingresando al reactor
T_in = Tinop + delta_Tin * t; % Temperatura del reactivo A
ingresando al reactor
T_c = Tcop + delta_Tc * Tcop * t; % Temperatura del refrigerante
% Simulación en Simulink
Tarr = 17000; % Tiempo de arranque (segundos)
Tsim = 35000; % Tiempo de simulación (segundos)
Tarr1 = 12000; % Tiempo de arranque (segundos)
Tsim1 = 25000; % Tiempo de simulación (segundos)
sim('DiagramasAB', Tsim);
t1 = tout; % Tiempo exportado desde Simulink
t2 = tout;
% Variables exportadas
CA_in1 = entradas(:,1); % Concentración de A ingresando al
reactor
T_in1 = entradas(:,2); % Temperatura del reactivo A ingresando
al reactor
T_c1 = entradas(:,3); % Temperatura del refrigerante
CA1 = salidas(:,1); % Concentración de A en el reactor
T1 = salidas(:,2); % Temperatura absoluta del producto al
interior del reactor
% Graficas de las salidas principales del sistema
figure;
subplot(2,1,1);
plot(t, CA, 'm', 'LineWidth', 1.5);
title('Concentración de A en el Reactor');
xlabel('Tiempo (h)');
ylabel('Concentración CA (kmol/m^3)');
grid on;
subplot(2,1,2);
plot(t, T, 'g', 'LineWidth', 1.5);
title('Temperatura del Reactor');
xlabel('Tiempo (h)');
ylabel('Temperatura T (K)');
grid on;
% Graficas Adicionales
figure;
% Gráfica 1: Concentración de A ingresando al reactor vs tiempo
subplot(3,1,1);
plot(t1, CA_in1, 'b', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('CA_{in} (kmol/m^3)');
title('Concentración de A ingresando al reactor vs Tiempo');
grid on;
Sistemas de Control
% Gráfica 2: Temperatura del reactivo A ingresando al reactor vs
tiempo
subplot(3,1,2);
plot(t1, T_in1, 'b', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('T_{in} (K)');
title('Temperatura del reactivo A ingresando al reactor vs Tiempo');
grid on;
% Gráfica 3: Temperatura del refrigerante vs tiempo
subplot(3,1,3);
plot(t1, T_c1, 'y', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('T_c (K)');
title('Temperatura del refrigerante vs Tiempo');
grid on;
% Salidas
figure;
% Gráfica 1: Concentración del reactivo A en el reactor vs tiempo
subplot(2,1,1);
plot(t1, CA1, 'y', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('CA(t) (kmol/m^3)');
title('Concentración de A en el reactor vs Tiempo');
grid on;
% Gráfica 2: Temperatura absoluta del producto al interior del reactor
vs tiempo
subplot(2,1,2);
plot(t1, T1, 'g', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('T(t) (K)');
title('Temperatura absoluta del producto al interior del reactor vs
Tiempo');
grid on;
% Análisis de índices de desempeño de la respuesta escalón
% Utilizamos stepinfo para analizar las señales obtenidas
info_CA = stepinfo(CA, t); % Información de la respuesta de CA
info_T = stepinfo(T, t); % Información de la respuesta de T
info_CA_in = stepinfo(CA_in, t); % Información de la respuesta de
CA_in
info_T_in = stepinfo(T_in, t); % Información de la respuesta de
T_in
info_T_c = stepinfo(T_c, t); % Información de la respuesta de T_c
% Mostrar índices de desempeño en la consola
disp('Índices de Desempeño para la Concentración de A en el Reactor
(CA):');
disp(info_CA);
disp('Índices de Desempeño para la Temperatura del Reactor (T):');
disp(info_T);
disp('Índices de Desempeño para la Concentración de A ingresando al
reactor (CA_in):');
disp(info_CA_in);
Sistemas de Control
disp('Índices de Desempeño para la Temperatura del reactivo A
ingresando al reactor (T_in):');
disp(info_T_in);
disp('Índices de Desempeño para la Temperatura del refrigerante
(T_c):');
disp(info_T_c);
%% PARTE 1B
%1. Hallar los puntos de operación del sistema usando vpasolve
% Parámetros
syms CA T s
Q = 2; % Caudal de la corriente caliente (m3/s)
V = 1.5; % Caudal de la corriente fría (m3/s)
deltaH = -5960; % Cambio de entalpía en la corriente caliente
(J/kg)
UA = 300; % Coeficiente de transferencia de calor
multiplicado por el area del reactor (W/K)
R = 1.985875; % Constante de los gases ideales (kcal/(kmol·K))
E = 11843; % Energía de activación (kcal/kmol)
k0 = 34930800; % Factor pre-exponencial (1/h)
ko = 34930800;
rhoCp = 500; % Producto densidad-capacidad calorífica
(kcal/(K·m^3))
CAi = CAiop + ruido; % Agregar el ruido a la concentración de A
Tin = Tinop; % Temperatura de entrada del reactivo A (K)
Tc = Tcop; % Temperatura del refrigerante (K)
% Ecuaciones no lineales
dCAdt = (Q/V)*(CAi - CA) - k0*exp(-E/(R*T)) * CA;
dTdt = (Q/V)*(Tin - T) + (UA/(rhoCp*V))*(Tc - T) +
(deltaH/rhoCp)*k0*exp(-E/(R*T)) * CA;
% Supón que ya has definido las ecuaciones no lineales
eq1 = (Q/V)*(CAi - CA) - k0*exp(-E/(R*T)) * CA == 0;
eq2 = (Q/V)*(Tin - T) + (UA/(rhoCp*V))*(Tc - T) +
(deltaH/rhoCp)*k0*exp(-E/(R*T)) * CA == 0;
% Resolver las ecuaciones no lineales usando vpasolve
sol = vpasolve([eq1, eq2], [CA, T], [CAi, Tin]); % [CA, T] son las
incógnitas
% Convertir las soluciones simbólicas a números
CA_op = double([Link]);
T_op = double(sol.T);
% Verifica la solución
disp(sol);
Sistemas de Control
Gráficas:
Sistemas de Control
Parte Lineal
Subsistema interiormente
Sistemas de Control
Código:
%% PARTE 1B
%1. Hallar los puntos de operación del sistema usando vpasolve
% Parámetros
syms CA T s
Q = 2; % Caudal de la corriente caliente (m3/s)
V = 1.5; % Caudal de la corriente fría (m3/s)
deltaH = -5960; % Cambio de entalpía en la corriente caliente
(J/kg)
UA = 300; % Coeficiente de transferencia de calor
multiplicado por el area del reactor (W/K)
R = 1.985875; % Constante de los gases ideales (kcal/(kmol·K))
E = 11843; % Energía de activación (kcal/kmol)
k0 = 34930800; % Factor pre-exponencial (1/h)
ko = 34930800;
rhoCp = 500; % Producto densidad-capacidad calorífica
(kcal/(K·m^3))
CAi = CAiop + ruido; % Agregar el ruido a la concentración de A
Tin = Tinop; % Temperatura de entrada del reactivo A (K)
Tc = Tcop; % Temperatura del refrigerante (K)
% Ecuaciones no lineales
dCAdt = (Q/V)*(CAi - CA) - k0*exp(-E/(R*T)) * CA;
dTdt = (Q/V)*(Tin - T) + (UA/(rhoCp*V))*(Tc - T) +
(deltaH/rhoCp)*k0*exp(-E/(R*T)) * CA;
% Supón que ya has definido las ecuaciones no lineales
eq1 = (Q/V)*(CAi - CA) - k0*exp(-E/(R*T)) * CA == 0;
eq2 = (Q/V)*(Tin - T) + (UA/(rhoCp*V))*(Tc - T) +
(deltaH/rhoCp)*k0*exp(-E/(R*T)) * CA == 0;
% Resolver las ecuaciones no lineales usando vpasolve
sol = vpasolve([eq1, eq2], [CA, T], [CAi, Tin]); % [CA, T] son las
incógnitas
% Convertir las soluciones simbólicas a números
CA_op = double([Link]);
T_op = double(sol.T);
% Verifica la solución
disp(sol);
%% Linealizar las ecuaciones no lineales del sistema utilizando los
puntos de operacion
% Derivadas parciales de la ecuación de concentración (dCAdt)
dCAdt_CA = diff(dCAdt, CA); % Derivada con respecto a CA
dCAdt_T = diff(dCAdt, T); % Derivada con respecto a T
% Derivadas parciales de la ecuación de temperatura (dTdt)
dTdt_CA = diff(dTdt, CA); % Derivada con respecto a CA
dTdt_T = diff(dTdt, T); % Derivada con respecto a T
% Evaluación en el punto de operación (CA_op, T_op)
dCAdt_CA_op = double(subs(dCAdt_CA, [CA, T], [CA_op, T_op]));
dCAdt_T_op = double(subs(dCAdt_T, [CA, T], [CA_op, T_op]));
dTdt_CA_op = double(subs(dTdt_CA, [CA, T], [CA_op, T_op]));
dTdt_T_op = double(subs(dTdt_T, [CA, T], [CA_op, T_op]));
% Ecuaciones linealizadas (en términos de desviaciones)
Sistemas de Control
% Desviaciones de CA y T
CA_dev = CA - CA_op;
T_dev = T - T_op;
% Ecuación linealizada para dCAdt
dCAdt_lin = dCAdt_CA_op * CA_dev + dCAdt_T_op * T_dev;
% Ecuación linealizada para dTdt
dTdt_lin = dTdt_CA_op * CA_dev + dTdt_T_op * T_dev;
% Funciones de transferencia del sistema (en el dominio de Laplace)
% Función de transferencia para CA(s)
G_CA = (dCAdt_CA_op * s + dCAdt_T_op) / (s + (Q/V + k0*exp(-
E/(R*T_op))));
% Función de transferencia para T(s)
G_T = (dTdt_CA_op * s + dTdt_T_op) / (s + (Q/V + UA/(rhoCp*V) -
deltaH/(rhoCp) * k0 * exp(-E/(R*T_op)) * (E/(R*T_op^2))));
% Definir las ecuaciones linealizadas
dCAdt =
(61365869812029708840186307560829/2535301200456458802993406410752) -
(5320510841784427 * T)/144115188075855872 - (3128339998353961 *
CA)/2251799813685248;
dTdt =
(204847116040504692153090945040053/316912650057057350374175801344) -
(4894064821336779 * T)/2251799813685248 - (6004830966166693 *
CA)/9007199254740992;
% Simplificar las ecuaciones
dCAdt_simplified = simplify(dCAdt);
dTdt_simplified = simplify(dTdt);
% Funciones de transferencia
CA_s = (-((3128339998353961 * s)/2251799813685248 +
5320510841784427/144115188075855872))/(s +
3128339998353961/2251799813685248);
T_s = (-((6004830966166693 * s)/9007199254740992 +
4894064821336779/2251799813685248))/(s +
8012739659711277/4503599627370496);
% Simplificar las funciones de transferencia
CA_s_simplified = simplify(CA_s);
T_s_simplified = simplify(T_s);
% Graficas con respecto al tiempo - Sistema Lineal
% Variables exportadas
CA_in2 = entradas_lineal(:,1); % Concentración de A ingresando al
reactor
T_in2 = entradas_lineal(:,2); % Temperatura del reactivo A
ingresando al reactor
T_c2 = entradas_lineal(:,3); % Temperatura del refrigerante
CA2 = salidas_lineal(:,1); % Concentración de A en el reactor
T2 = salidas_lineal(:,2); % Temperatura absoluta del producto
al interior del reactor
% Entradas
figure;
% Gráfica 1: Concentración de A ingresando al reactor vs tiempo
subplot(3,1,1);
Sistemas de Control
plot(t2, CA_in2, 'm', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('CA_{in} (kmol/m^3)');
title('Concentración de A ingresando al reactor vs Tiempo');
grid on;
% Gráfica 2: Temperatura del reactivo A ingresando al reactor vs
tiempo
subplot(3,1,2);
plot(t2, T_in2, 'g', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('T_{in} (K)');
title('Temperatura del reactivo A ingresando al reactor vs Tiempo');
grid on;
% Gráfica 3: Temperatura del refrigerante vs tiempo
subplot(3,1,3);
plot(t2, T_c2, 'c', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('T_c (K)');
title('Temperatura del refrigerante vs Tiempo');
grid on;
% Salidas
figure;
% Gráfica 1: Concentración del reactivo A en el reactor vs tiempo
subplot(2,1,1);
plot(t2, CA2, 'm', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('CA(t) (kmol/m^3)');
title('Concentración de A en el reactor vs Tiempo');
grid on;
% Gráfica 2: Temperatura absoluta del producto al interior del reactor
vs tiempo
subplot(2,1,2);
plot(t2, T2, 'g', 'LineWidth', 1.5);
xlabel('Tiempo (h)');
ylabel('T(t) (K)');
title('Temperatura absoluta del producto al interior del reactor vs
Tiempo');
grid on;
% Análisis de Índices de desempeño de la respuesta escalón - Sistema
Lineal
% Funciones de transferencia previamente calculadas
% 1. Función de transferencia para la concentración de A (CA)
CA_s_num = [-3128339998353961, -5320510841784427];
CA_s_den = [1, 3128339998353961/2251799813685248];
CA_s_tf = tf(CA_s_num, CA_s_den);
% 2. Función de transferencia para la temperatura del reactor (T)
T_s_num = [-6004830966166693, -4894064821336779];
T_s_den = [1, 8012739659711277/4503599627370496];
T_s_tf = tf(T_s_num, T_s_den);
% 3. Función de transferencia para la concentración de A ingresando al
reactor (CA_in)
CA_in_s_num = [1]; % Entrada directa para perturbaciones
CA_in_s_den = [1]; % Sin dinámica interna
CA_in_s_tf = tf(CA_in_s_num, CA_in_s_den);
Sistemas de Control
% 4. Función de transferencia para la temperatura del reactivo
ingresando al reactor (T_in)
T_in_s_num = [UA / (rhoCp * V)];
T_in_s_den = [1, (Q/V + UA / (rhoCp * V))];
T_in_s_tf = tf(T_in_s_num, T_in_s_den);
% 5. Función de transferencia para la temperatura del refrigerante
(T_c)
T_c_s_num = [UA / (rhoCp * V)];
T_c_s_den = [1, (Q/V + UA / (rhoCp * V))];
T_c_s_tf = tf(T_c_s_num, T_c_s_den);
% Generar las respuestas escalón
[CA_step_response, t_CA] = step(CA_s_tf);
[T_step_response, t_T] = step(T_s_tf);
[CA_in_step_response, t_CA_in] = step(CA_in_s_tf);
[T_in_step_response, t_T_in] = step(T_in_s_tf);
[T_c_step_response, t_T_c] = step(T_c_s_tf);
% Cálculo de Índices de desempeño
info_CA = stepinfo(CA_s_tf);
info_T = stepinfo(T_s_tf);
info_CA_in = stepinfo(CA_in_s_tf);
info_T_in = stepinfo(T_in_s_tf);
info_T_c = stepinfo(T_c_s_tf);
% Mostrar Índices de desempeño
disp('Índices de Desempeño para la Concentración de A en el Reactor
(CA):');
disp(info_CA);
disp('Índices de Desempeño para la Temperatura del Reactor (T):');
disp(info_T);
disp('Índices de Desempeño para la Concentración de A ingresando al
reactor (CA_in):');
disp(info_CA_in);
disp('Índices de Desempeño para la Temperatura del reactivo A
ingresando al reactor (T_in):');
disp(info_T_in);
disp('Índices de Desempeño para la Temperatura del refrigerante
(T_c):');
disp(info_T_c);
% Mostrar resultados
% Imprimir los puntos de operación de CA y T
disp(['Punto de operación de CA: ', num2str(CA_op)]);
disp(['Punto de operación de T: ', num2str(T_op)]);
disp('Ecuaciones linealizadas:');
%disp(['dCAdt linealizada: ', char(dCAdt_lin)]);
%disp(['dTdt linealizada: ', char(dTdt_lin)]);
disp('dCAdt linealizada simplificada:');
disp(dCAdt_simplified);
disp('dTdt linealizada simplificada:');
disp(dTdt_simplified);
Sistemas de Control
disp('Funciones de transferencia:');
%disp(['Función de transferencia para CA(s): ', char(G_CA)]);
%disp(['Función de transferencia para T(s): ', char(G_T)]);
disp('Función de transferencia para CA(s) simplificada:');
disp(CA_s_simplified);
disp('Función de transferencia para T(s) simplificada:');
disp(T_s_simplified);
Gráficas del Modelo Lineal
Sistemas de Control
Sistemas de Control
Conclusiones:
• El modelo no lineal fue representado por ecuaciones diferenciales que describen
el comportamiento del reactor en términos de concentración de A y temperatura,
con perturbaciones en las entradas como la concentración del reactivo y la
temperatura del refrigerante. Este modelo resultó ser más detallado, pero al mismo
tiempo más complejo, debido a su naturaleza no lineal.
• Los resultados obtenidos del modelo lineal se aproximaron muy bien a los del
modelo no lineal, especialmente para pequeñas desviaciones alrededor de los
puntos de operación. Esto validó que la linealización es una herramienta útil para
sistemas que operan cerca de un punto de operación, donde las variaciones son
pequeñas.
• La linealización fue realizada utilizando las derivadas parciales de las ecuaciones
no lineales respecto a las variables de interés (concentración y temperatura). Esto
resultó en un conjunto de ecuaciones lineales que describen las desviaciones de
las variables respecto a los puntos de operación.
Recomendaciones:
➢ Aunque la linealización es útil, es importante realizar una validación constante del
modelo con datos experimentales o simulaciones detalladas para asegurarse de
que las aproximaciones lineales sean precisas. En caso de grandes desviaciones
del punto de operación, el modelo no lineal debe ser utilizado
➢ El modelo lineal es especialmente útil cuando se desea analizar el sistema en
términos de control, ya que permite el uso de técnicas estándar como las funciones
de transferencia y los diagramas de Bode (que se estudiarán a futuro en el apartado
teórico de la materia) para analizar la estabilidad y el comportamiento transitorio
del sistema.
