%% Linealización de un modelo de tanque de mezclado
% © 2019, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Este código ejecutó correctamente en Matlab R2019a
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/linmix.html
% Motivación y objetivos: El objetivo de este vídeo es modelar un tanque de mezclado, incorporando no-linealidades en caudal de salida (raíz cuadrada) y concentración (cociente). Se linealizará y se comparará la simulación del modelo linealizado (vbles. incrementales) con el modelo original.
%
%% Modelado
%%
%% Ec. físicas

syms v_a v_b dv_adt dv_bdt %2 estados y 2 deriv. estado
syms x_a x_b q_s h %4 otras variables que aparecen en ecs modelo
syms q_ain q_bin %2 variables de entrada
S=0.85;k_out=4; %parámetros constantes
Modelo=[ h==(v_a+v_b)/S; q_s==k_out*sqrt(h); dv_adt==-q_s*x_a+q_ain;...
         dv_bdt==-q_s*x_b+q_bin; x_a==v_a/(v_a+v_b); x_b==1-x_a];

%% Linealización
%% 1.-Calculamos punto de funcionamiento
% Por diseño, nos dicen que caudales entrada son 3 y 1 unidades. Resolvemos numéricamente ec. estado en equilibrio:

sol=vpasolve([Modelo;dv_adt==0;dv_bdt==0;q_ain==3;q_bin==1])
q_apf=eval(sol.q_ain)
q_bpf=eval(sol.q_bin)
v_apf=eval(sol.v_a)
v_bpf=eval(sol.v_b)
qs_pf=eval(sol.q_s)
xa_pf=eval(sol.x_a)
xb_pf=eval(sol.x_b)
h_pf=eval(sol.h)

%% Modelo linealizado en representación interna:
% La linealización de las seis ecuaciones del modelo sería hacer deriv. parciales por incremento, pero para simulación, control, etc. conviene usar repr. interna normalizada, estados v_a y v_b, entradas q_{a,in}, q_{b,in}.

ReprInterna=solve(Modelo,dv_adt,dv_bdt,q_s,x_a,x_b,h);
derivEstado=[ReprInterna.dv_adt;ReprInterna.dv_bdt]

% Sólo nos interesa simular composición x_a y nivel h (elegidas arbitrariamente):

EcSalida=[ReprInterna.x_a;ReprInterna.h]
%linealizamos
%1: haciendo derivadas parciales
Asym=simplify(jacobian(derivEstado,[v_a v_b]))
Bsym=simplify(jacobian(derivEstado,[q_ain q_bin]))
Csym=simplify(jacobian(EcSalida,[v_a v_b]))
Dsym=simplify(jacobian(EcSalida,[q_ain q_bin]))
%2: las deriv. parciales son evaluadas en punto de funcionamiento
Alin=eval(subs(Asym,sol))
Blin=eval(subs(Bsym,sol))
Clin=eval(subs(Csym,sol))
Dlin=eval(subs(Dsym,sol))

% El modelo linealizado tiene esta ecuación de estado:

dxdt_lin=@(incx,incu) Alin*incx+Blin*incu;

% y esta ecuación de salida:

EcSalidaLinealizado=@(incx,incu) Clin*incx+Dlin*incu;

%% Simulación del modelo original y el linealizado
% Las entradas (caudales de A y B) tendrán la siguiente forma:

Tiempos=0:0.05:10;
q_ain_simula=q_apf+0.25*sign(sin(Tiempos));
q_bin_simula=q_bpf+0.4*sign(sin(Tiempos.^2*0.19));
vector_entradas_sim=[q_ain_simula;q_bin_simula];
plot(Tiempos,vector_entradas_sim), grid on
legend('q_ain','q_bin','Location','best'), xlabel('tiempo'), ylabel('Caudal entrada')

%
%% Simulación de la ecuación de estado
% La ecuación NO lineal de estado es:

dxdt=matlabFunction(derivEstado,'Vars',{[v_a;v_b],[q_ain;q_bin]});

%
% Preparamos ode45 para simular el no lineal:

CondInicial=[0.05;0.2];
modeloparasimular=@(t,estado) dxdt(estado, interp1(Tiempos,vector_entradas_sim',t)');
[tiempossimulados,estadosimulado]=ode45(modeloparasimular,Tiempos,CondInicial);

% Preparamos ode45 para simular el linealizado (entradas, estados y salidas son incrementales):

entrada_incre=vector_entradas_sim-[q_apf;q_bpf];
CondInic_incre=CondInicial-[v_apf;v_bpf];
modelsimulalinealiz=@(t,estado) dxdt_lin(estado, interp1(Tiempos,entrada_incre',t)');
[tiempossimulaL,inc_estadosimul]=ode45(modelsimulalinealiz,Tiempos,CondInic_incre);

% Comparamos ambos resultados (trayectorias vbles. de estado v_a, v_b):

plot(tiempossimulados,estadosimulado), grid on, hold on
plot(tiempossimulaL,inc_estadosimul+[v_apf v_bpf]), hold off %deshacemos incrementos!
yline(v_apf,'-.');yline(v_bpf,'-.','Color','#559977');
legend('volumen de a','volumen de b','v_alin','v_blin','v_apf','v_bpf','Location','best'), xlabel('tiempo')

%% Evaluación de la ecuación de salida
% Una vez simuladas las ecuaciones de estado, hay que sustituirlas en la de salida.
% En el modelo no-lineal:

EcSal_numerica=matlabFunction(EcSalida,'Vars',{[v_a;v_b]});
salidassimuladas=EcSal_numerica(estadosimulado');

% En el modelo linealizado resulta

inc_salidas_linealizado=EcSalidaLinealizado(inc_estadosimul',entrada_incre);

% Lo visualizamos todo en una gráfica para comparar:

plot(Tiempos,salidassimuladas), grid on, hold on
plot(Tiempos,inc_salidas_linealizado+[xa_pf;h_pf]), hold off %deshacemos incrementos
yline(xa_pf,'-.','Color','#556699');yline(h_pf,'-.','Color','#998855');
legend('x_a (no lineal)','h (no lineal)','x_a (lin.)','h (lin.)','x_apf','h_pf','Location','best'), xlabel('tiempo')

% Si almacenamos el modelo linealizado como un "objeto sistema lineal" de la Control System Toolbox, podremos analizarlo y controlarlo (no en objetivos de este material) con las herramientas de dicha toolbox:

TanqueMezcladoLin=ss(Alin,Blin,Clin,Dlin);
TanqueMezcladoLin.InputName={'q_ain','q_bin'};
TanqueMezcladoLin.OutputName={'x_a','h'};
TanqueMezcladoLin.StateName={'vol_a','vol_b'};
size(TanqueMezcladoLin)
TanqueMezcladoLin
save ModeloTanqueLin.mat TanqueMezcladoLin

%