%% Sistema térmico: modelado y análisis de propiedades
% © 2021, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/termem.html
% Este código funcionó correctamente con Matlab R2020b
% Objetivo:  introducir en Matlab el modelo en representación interna de este sistema térmico. Calcular un punto de funcionamiento. Analizar las propiedades en tiempo y frecuencia.
%
% Tabla de Contenidos
% Modelado
% Análisis de propiedades
%     Ganancia, punto de funcionamiento
% Respuesta temporal
%         Modos de la respuesta temporal
% Respuesta en frecuencia
%%
%% Modelado

%parámetros constantes del sistema
k_hab_ext1=.24; k_hab_ext2=.27;
k_hab_c=.33; k_techo1=0.4; k_techo2=0.35; k_techo_c=0.08; k_tejado=0.55;
c_h1=210; c_h2=180; c_centr=45; c_buh=420;

MassMatr=diag([c_h1  c_h2 c_centr c_buh]) ;
%usamos la sintaxis Matlab   M\N <==> inv(M)*N
A=MassMatr \ [-(k_hab_ext1+k_hab_c+k_techo1)  0   k_hab_c     k_techo1; ...
                   0   -(k_hab_ext2+k_hab_c+k_techo2)   k_hab_c   k_techo2; ...
                k_hab_c   k_hab_c    -(2*k_hab_c+k_techo_c)    k_techo_c; ...
               k_techo1   k_techo2      k_techo_c   -k_techo1-k_techo2-k_techo_c-k_tejado]
B=MassMatr \ [1 0 k_hab_ext1; 0 1 k_hab_ext2; 0 0 0; 0 0 k_tejado]
C=eye(4); %sacamos todas las temperaturas como salida, por el momento.
casa=ss(A,B,C,zeros(4,3));
casa.InputName={'Q1','Q2','Tamb'};
casa.StateName={'T1','T2','Tc','Tb'};
casa.OutputName=casa.StateName; %al ser C=I salidas y estados coinciden

%% Análisis de propiedades
%% Ganancia, punto de funcionamiento

Gan=dcgain(casa)
gancasaT1T2_Q1Q2=Gan(1:2,1:2)
[U,S,V]=svd(gancasaT1T2_Q1Q2)
cond(S)
gancasaT1T2_amb=Gan(1:2,3)

% Por ejemplo, podemos usar la matriz de ganancia para calcular puntos de funcionamiento.
% y=Gu+G_{pert}\cdot T_a, busco 0=Gu+G_{pert}\cdot T_a, esto es u=-G^{-1}G_{pert}\cdot T_a
% La matriz (ganancia, rég. permanente) que indica los incrementos de actuadores Q_1, Q_2 que cancelarían un grado de incremento de T_a es:

Prealim=-inv(gancasaT1T2_Q1Q2)*gancasaT1T2_amb

% Los valores de actuadores necesarios para compensar 15 grados de frío exterior son:

act_15=Prealim*(-15)

%% Respuesta temporal

step(casa), grid on
step(casa(1,:),casa(2,:),casa(3,:),casa(4,:),4000), grid on, legend('hab.1','hab.2','central','buhardilla')

%
%% Modos de la respuesta temporal

[V,D]=eig(A);
[S,idx]=esort(diag(D));
S'
V(:,idx) %modos
-4./S' %tiempos de establecimiento de cada modo

%
%% Respuesta en frecuencia
% 12 diagramas de Bode no son demasiado informativos...

bodemag(casa), grid on %sólo amplitud

% Valores singulares de manipuladas o perturbación a controladas:

sigma(casa(1:2,1:2), casa(1:2,3)), grid on
legend("Q->T ºC/Kw", "Ta->T ºC/ºC")

% Nota: saber si se pueden seguir ciertas referencias o cancelar ciertos incrementos de temperatura ambiente de una cierta amplitud necesita de "escalado", no en los objetivos de este material.
%%