%% Modelado, linealización y simulación (lineal/no lineal) de un 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/linmixqs.html
% Motivación y objetivos: El objetivo de este vídeo es modelar un tanque de mezclado cuyo caudal de salida es forzado a un cierto valor por la demanda aguas abajo. Se linealizará y se comparará la simulación del modelo linealizado con el modelo original.
%%
%% Ec. físicas

syms v_a v_b  dv_adt dv_bdt  %2 estados y sus derivadas
syms x_a x_b h  %3 otras variables
syms q_s q_ain q_bin %3 variables de entrada
S=0.85; %parámetros constantes (sección transversal)

% El modelo tiene 5 ecs.

Modelo=[ h==(v_a+v_b)/S; 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]

% Bien planteado porque hay 8 señales, con 3 entradas (q_a, q_b, q_s) y cinco variables a calcular (v_a, v_b, x_a, x_b, h).
% Para repr. interna, consideramos estados (v_a,v_b), con lo que podemos despejar 5 letras (dv_adt,dv_bdt,x_a,x_b,h) en función de otras 5: los 2 estados + las 3 entradas.

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

% Sólo concentración de A, x_A y nivel h serán de interés como salidas:

EcSalida=simplify([ReprInterna.x_a;ReprInterna.h])

% La ecuación NO lineal de estado es:

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

% y la ecuación de salida (numérica):

EcSal_num=matlabFunction(EcSalida,'Vars',{[v_a;v_b]});

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

EcsEquilibrio=[Modelo;dv_adt==0;dv_bdt==0;q_ain==3;q_bin==1;q_s==4];
solPF=solve(EcsEquilibrio) %probar con q_s=3???
eval(solPF.h)

% Pero la solución no es correcta... hay infinitas... si sospechamos de múltiples equilibrios, haríamos:

solPF=solve(EcsEquilibrio,  'ReturnConditions', true)

% Donde devuelve infinitas soluciones de altura dependientes de un paŕametro z\neq 0:

solPF.h
solPF.parameters
solPF.conditions

% En conclusión, las 10 ecs. en equilibrio NO son independientes y "una" podría eliminarse... o añadirse otra:

EcsEquilibrio=[h==1.25;dv_adt==0;dv_bdt==0;Modelo;q_ain==3;q_bin==1];
solPF=solve(EcsEquilibrio,  'ReturnConditions', true)
solPF.parameters
solPF.conditions

% Esto ya es, por tanto, una solución "única", sin "condiciones" ni parámetros.
% Recalculamos, pues:

solPF=solve(EcsEquilibrio); 
q_apf=eval(solPF.q_ain)
q_bpf=eval(solPF.q_bin)
v_apf=eval(solPF.v_a)
v_bpf=eval(solPF.v_b)
qs_pf=eval(solPF.q_s)
xa_pf=eval(solPF.x_a)
xb_pf=eval(solPF.x_b)
h_pf=eval(solPF.h)

%% linealizamos en el punto elegido

%1: haciendo derivadas parciales
Asym=simplify(jacobian(derivEstado,[v_a v_b]))
Bsym=simplify(jacobian(derivEstado,[q_ain q_bin q_s]))
Csym=simplify(jacobian(EcSalida,[v_a v_b]))
Dsym=simplify(jacobian(EcSalida,[q_ain q_bin q_s]))
%2: evaluadas en el punto de funcionamiento
Alin=eval(subs(Asym,solPF))
Blin=eval(subs(Bsym,solPF))
Clin=eval(subs(Csym,solPF))
Dlin=eval(subs(Dsym,solPF))

% El linealizado es:

dxdt_lin=@(x,u) Alin*x+Blin*u;
EcSal_lin=@(x,u) Clin*x+Dlin*u;

%% Simulación
% Una vez disponibles ecuaciones de estado y de salida, vamos a simular el llenado del depósito... para evitar dividir por cero, comenzaremos con una "pequeña" cantidad de líquido en él (las ecuaciones no son válidas si no está "cubierto" todo el orificio de salida, etc.)
% Las entradas (caudales de A y B) tendrán la siguiente forma:

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

%
%% Simulación de la ecuación de estado
%
% Preparamos ode45 no lineal:

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

% Preparamos ode45  lineal:

entrada_incrementos=vector_entradas_sim-[q_apf;q_bpf;qs_pf];
CondIni_inc=CondIni-[v_apf;v_bpf];
modeloparasimularlin=@(t,estado) dxdt_lin(estado, interp1(Tiempos,entrada_incrementos',t)');
[tiempossimuladosL,inc_estadosimulado]=ode45(modeloparasimularlin,Tiempos,CondIni_inc);

% Representamos gráficamente:

plot(tiempossimulados,estadosimulado), grid on, hold on
plot(tiempossimuladosL,inc_estadosimulado+[v_apf v_bpf]), hold off
legend('volumen de a','volumen de b','v_alin','v_blin','Location','best'), xlabel('tiempo')

%% Evaluación de la ecuación de salida

salidassimuladas=EcSal_num(estadosimulado');
inc_salidas_linealizado=EcSal_lin(inc_estadosimulado',entrada_incrementos);
plot(Tiempos,salidassimuladas), grid on, hold on
plot(Tiempos,inc_salidas_linealizado+[xa_pf;h_pf]), grid on, hold on
legend('x_a (no lineal)','h (nivel)','x_a (lin.)','h (lin.)','Location','best'), xlabel('tiempo')

% Nota: las ecuaciones de estos modelos dejan de tener sentido si algún volumen almacenado de a o b se hace negativo, o si el nivel y concentraciones difieren mucho del punto de linealización (nivel 1.25, concentración de A 0.75 y de B 0.25)
%
% Uso posterior del modelo linealizado:
% Si almacenamos el modelo linealizado como un "objeto sistema lineal" de la Control System Toolbox, podremos analizarlo y controlarlo con las herramientas de dicha toolbox:

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

%