%% Ejemplo modelado/análisis de temperatura en una tubería: ruido coloreado de entrada
% © 2019, 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/tuborui.html
% Este código ejecutó sin errores en Matlab R2019a
% Motivación/objetivos: comprender el significado de una perturbación "ruido coloreado" a un proceso físico; simular realización, media, varianza del mismo.
%% Modelado
%
% Suponiendo T_g (temperatura gas) y T_a (temperatura ambiente) como entradas, tendríamos unas constantes de transmisión de calor y capacidades caloríficas que darían lugar a la ecuación de estado:
% \frac{d}{dt}\left(\begin{array}{c}T_s \\ T_t \end{array}\right)=
% \left(\begin{array}{cc}(-k_{st}-k_{sa})/C_s &k_{st}/C_s \\ k_{st}/C_t &(-k_{st}-k_{ta}-k_{tg})/C_t \end{array}\right)
% \left(\begin{array}{c}T_s \\ T_t \end{array}\right)
% +
% \left(\begin{array}{cc}0 &k_{sa} /C_s\\ k_{tg}/C_t &k_{ta}/C_t\end{array}\right)
% \left(\begin{array}{c}T_g \\ T_a \end{array}\right)
% Nota: este es un modelo simplificado de parámetros concentrados; la transmisión de calor en una tubería se rige por ecuaciones en derivadas parciales fuera de los objetivos de este modelado. Supondremos, por tanto, que la tubería es "corta" y que las pérdidas de entalpía del gas (calor transmitido a la tubería) son despreciables.
% Vamos a dar valores numéricos

kst=.3; ktg=.25; kta=.08; ksa=0.01; Cs=.4; Ct=1.5;
A=diag([1/Cs 1/Ct])*[-kst-ksa kst; kst -ktg-kta-kst];
B=diag([1/Cs 1/Ct])*[0 ksa; ktg kta];
C=eye(2); %queremos simular las dos temperaturas (sonda y tubo)
sys=ss(A,B,C,0);
sys.InputName={'Tgas','Tambiente'};
sys.OutputName={'Tsonda','Ttuberia'};
bodemag(sys(1,:),sys(2,:),logspace(-2.5,0.5)), grid on, legend('Tsonda','Ttubería')

%
%% Incorporación de modelo de perturbaciones
% Temperatura del gas va a ser ruido "coloreado":

anchobandagas=0.16; s=tf('s');
GeneradorTgas=1/(s+anchobandagas)^2; %multipl. por ruido blanco...

% Por simplicidad, descartemos variaciones de Tamb...

sys2=ss([sys(:,1);1]*GeneradorTgas);
sys2.OutputName{3}='Tgas';
sys2

% Discretizamos:

Ts=0.2; w_nyq=pi/0.2
sysd=c2d(sys2,Ts,'zoh'); 
%existe una forma teóricamente mejor de discretizar ruidos continuos,
%omitida por brevedad [es la que hace el comando lqed]

%% Análisis de propiedades

Q=300; %varianza ruido de proceso
Zba=gram(sysd,'c')*Q %vzas-covarianzas estado en bucle abierto (estacionarias):

% Si y=Cx+v, entonces E[yy^T]=E[(Cx+v)(x^TC^T+v^T)]=E[Cxx^TC+Cxv^T+vx^TC^T+vv^T]=CE[xx^T]C^T+E[vv^T]:

VzaSalidas=sysd.c*Zba*sysd.c' %suponiendo que no hay ruido medida, salidas "limpias"
dt_salidas=sqrt(diag(VzaSalidas))'
dt_Tgas=dt_salidas(3) %desviación. tip. temperatura gas a priori
norm(sysd(3)*sqrt(Q),2) %también es la norma 2 de un sistema

%% Simulación
% Vamos a simular una realización (trayectoria) cualquiera:

N=2000;Tiempo=(0:(N-1))*Ts;
perturbacionRuidoBlanco=randn(1,N)*sqrt(Q);
Y=lsim(sysd,perturbacionRuidoBlanco);
plot(Tiempo,Y), grid on %luego superpondremos los intervalos de confianza

% Si la perturbación tiene media cero, la simulación de la ecuación de medias es simplemente cero. Vamos a simular la ecuación de varianzas.
%% Ecuación de varianzas

SigmaX=zeros(4,4); %la simulación empezó en estado = 0, conocido perfectamente.
std_Y=zeros(3,N);

% Si x_+=Ax+Bv, entonces E[x_+x^T_+]=E[(Ax+Bv)(x^TA^T+v^TB^T)]=E[Axx^TA+Axv^TB^T+Bvx^TA^T+Bvv^TB^T]=AE[xx^T]A^T+BE[vv^T]B^T:

for i=1:N
    SigmaX=sysd.a*SigmaX*sysd.a'+sysd.b*Q*sysd.b';
    Sigmay=sysd.c*SigmaX*sysd.c';
   std_Y(:,i)=sqrt(diag(Sigmay));%guardamos desv típica
end
%representemos gráficamente la realización y sus intervalos de conf. (sólo Tgas)
hold on
plot(Tiempo,[2*std_Y(3,:);-2*std_Y(3,:)])
yline(2*dt_Tgas);yline(-2*dt_Tgas); %deviación típica estacionaria calculada antes
hold off
 legend('Temp. Sensor','Temp. Tubo','Temp. Gas','95% no estac','95% no estac','95% estac','Location','best')

%%