%% Simulación intermuestreo de un bucle cerrado con control por computador (sdlsim), efecto filtro antialias
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentaciones en vídeo:
% http://personales.upv.es/asala/YT/V/sdlsim1.html , http://personales.upv.es/asala/YT/V/sdlsim2.html .
% Este código funcionó correctamente con Matlab R2021b
% Objetivo: Comprender cómo realizar simulaciones intermuestreo de un sistema discreto de control usando sdlsim. Período de muestreo intencionamente "no demasiado bien elegido": pasarán cosas que se verán en la simulación intermuestreo.
%
% Tabla de Contenidos
% Modelo y controlador a probar
% Prestaciones en tiempo discreto
% Prestaciones intermuestreo
% Conclusiones
%
%% Modelo y controlador a probar

s=tf('s');
sysc=3.1/(s^2+1.1*s+3)+1.1*s/(s^2+s+60);
pole(sysc)
sysc.InputName='ureal';sysc.OutputName='y';
Ts=0.75; %0.1 o incluso 0.06 sería recomendable 
w_Nyquist=pi/Ts %Máximo ancho de banda de las señales sin "alias".

% Discreticemos el proceso:

FiltroAntiAlias=1/(s/(1.0*w_Nyquist)+1)^2 %paso bajo
%bodemag(FiltroAntiAlias)
sysdNoAA=c2d(sysc,Ts,'zoh');
sysd=c2d(FiltroAntiAlias*sysc,Ts,'zoh');
sysd.InputName='ureal';sysd.OutputName='y';
step(sysd,sysdNoAA,sysc), grid on, legend("Discreto+antialias","Discreto", "Continuo")
bodemag(sysd,sysdNoAA,sysc), grid on, legend("Discreto+antialias","Discreto no AA", "Continuo",Location="best")
z=tf('z',Ts);
Kd=0.7*(z^2-0.44*z+0.43)/(z-1)/(z+0.32);

%% Prestaciones en tiempo discreto
% Cerramos bucle con "connect", y podemos hacer step, bode, lsim, initial...
%

Kd.InputName='e';
Kd.OutputName='u';
s1=sumblk('e=r-y');s2=sumblk('ureal=u+du');
BC=connect(s1,s2,Kd,sysd,{'r','du'},{'y','u'});
BCnoAA=connect(s1,s2,Kd,sysdNoAA,{'r','du'},{'y','u'});
step(BC,BCnoAA), grid on
bodemag(BC,BCnoAA), grid on

%% Prestaciones intermuestreo
% El proceso controlado es realmente continuo, sdlsim puede simular respuesta intermuestreo.
% Debemos separar todo lo que es contínuo y todo lo que es discreto del diagrama original: esto sería en discreto
%
% pero realmente es:
%
% de modo que es una interconexión "lower Linear Fractional Transformation (LFT)" entre un sistema contínuo y uno discreto:
%
% Siendo w_1=(r,d_u) y z_1=(y,u). El sistema "2" tiene un único elemento (el regulador discreto), el sistema 1 tiene un diagrama de bloques, que podemos simplificar de varias formas:
% Opción 1: con el comando "connect":

PGc=connect(s1,s2,sysc,{'r','du','u'},{'y','u','e'});
size(PGc)

% Opción 2: operando "simbólicamente". En efecto, podemos deducir con facilidad que:
% \left(\begin{array}{c}y \\ u \\e\end{array}\right)=\left(\begin{array}{cc}0 &G(s) &G(s) \\ 0&0&1\\1 &-G(s) &-G(s)\end{array}\right)\cdot \left(\begin{array}{c}r \\ du \\ u\end{array}\right)
% con lo que podríamos programar esa matriz:

PGopcion2= [0 0 0;0 0 1;1 0 0]+[1;0;-1]*sysc*[0 1 1]; %escribo sysc sólo 1 vez para que no duplique estados
size(PGopcion2)
norm(PGc-PGopcion2) %son lo mismo

% Se recomienda introducir un "Antialias analógico" antes del muestreo para control digital... luego veremos por qué:

%PGc=blkdiag(1,1,FiltroAntiAlias)*PGc; %ANTIALIAS analógico

Tfin=60;
Tsimc=0:0.04:Tfin; %tiempos de "interpolación" continuos deseados

% Esto es un escalón de referencia

%UUc=ones(length(Tsimc),1)*[1 0]; %entradas generalizadas continuas (interpoladas): no necesariamente ha de ser "escalón", obviamente.

% Esto es escalón en "du":

%UUc=ones(length(Tsimc),1)*[0 1]; 

% Esto son otras cosas

%UUc=sin(0.15*Tsimc')*[0 1]; 
UUc=sin(8.6*Tsimc')*[0 1]; %alias, más allá de frec. Nyquist. Dará problemas!

% Simulamos

Y0=lsim(lft(PGc,0),UUc,Tsimc); %bucle abierto.

Y1=sdlsim(PGc,Kd,UUc,Tsimc);   %simulación HÍBRIDA (contínuo+discreto: sampled-data)

Tiempos_sdlsim=Y1{1};
Salidas_sdlsim=Y1{2};
Ts_sdlsim=Tiempos_sdlsim(2)-Tiempos_sdlsim(1)
N=round(Ts/Ts_sdlsim)
Tmuestras=Tiempos_sdlsim(1:N:end);
plot(Tiempos_sdlsim,Salidas_sdlsim(:,1),LineWidth=2), hold on, grid on
plot(Tsimc,Y0(:,1)), plot(Tmuestras,Salidas_sdlsim(1:N:end,1),'.',MarkerSize=18)
hold off
legend("intermuestreo","B.abto","Muestras (sin AA, con AA retrasaría un poco)",Location="best")
title("Salida")
plot(Y1{1},Y1{2}(:,2),'r',LineWidth=2), grid on
title("Acc. control")

%% Conclusiones
% Período de muestreo demasiado grande:
% -- Frecuencia Nyquist / AntiAliasing cercana a ancho banda bucle cerrado.
% -- Resonancia del proceso a alta frecuencia no capturada por el modelo. *Esto a veces es desconocido hasta que se experimenta y se ve que el control no funciona (error de modelado). La vida es dura.
% -- Existencia "experimental" de componentes de alta frecuencia en la señal contínua (si se dispone de equipo de más ancho de banda para poder verlos, si no, nunca se sabrá -- o se creerá que es "alias"--).
%
% Pero, es lo que hay... Al menos sdlsim permite simular todo eso. El filtro antialias elimina (en parte, no es ideal) los componentes de alta frecuencia para evitar que el control introduzca alias a frecuencias bajas que NO cancelan las perturbaciones a frecuencias altas y, por tanto, sea peor el remedio que la enfermedad...
% En resumen, bajar el período a 0.1, ancho banda Nyquist = 30 rad/s capturaría todo lo simulado aquí  sería posiblemente lo recomendable (requeriría rediseñar regulador discreto, pero el efecto "antialias" sería despreciable o quizás innecesario dado que el proceso controlado ahí ya filtra todo)...
% Al menos, en este caso, si no se puede bajar T_s, se debería perseguir unas especificaciones más lentas, entonces el antialias no molesta tanto a las prestaciones nominales (step) y si hay perturb. de alta frecuencia pues no se cancelan pero no se meten como "alias" de baja frecuencia en el bucle...
% NOTA: Esto es un "ejercicio académico"; en un bucle de control por computador bien diseñado, lo que se vea "en muestras" debe ser "todo lo importante" y la simulación "intermuestreo" no debería aportar nada de interés.