%% Caso ejemplo de un inversor monofásico PWM
% DEFINICIÓN DE PARAMETROS DE SIMULACIÓN
clear
close all
s=tf('s');
format long

fr=50;					% FRECUENCIA DESEADA DE LA TENSION DE SALIDA EN Hz
wo=2*pi*fr;             % FRECUENCIA ANGULAR DESEADA DE LA TENSION DE SALIDA EN rad/seg
Vbase=311;              % VALOR BASE DE TENSION PARA NORMALIZACIÓN
Ibase=50;               % VALOR BASE DE CORRIENTE PARA NORMALIZACIÓN
Vcc=450;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETROS DE LA PLANTA. Filtro y Carga
Lf=1000e-6;              % Indutancia de entrada en uH
Cf=20e-6;          	   % Capacitancia del filtro de salida en uF
wn=1/sqrt(Lf*Cf);        % Frecuencia angular de resonancia del filtro LC en rad/seg
fn=wn/2/pi;            % Frecuencia de resonancia del filtro LC en Hz
R=5;                   % Resistencia de carga 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Vector de estado: x = [vc  iL]'. Variable de salida: corriente en el inductor del filtro LC
% Matriz de estados A: dinámicas de la planta 
A=[-1/(R*Cf)  1/Cf
     -1/Lf    0  ];
     
B=[0;1/Lf];  C=[0 1];

F=0; D=0;

sistema_ssi=ss(A,B,C,D);
figure; step(Vbase*sistema_ssi)
legend('Respuesta del Modelo en SS con iL')

% Función de transferencia Y(s)/U(s) = C*inv[(sI - A)]*B
Gpi=zpk((R*Cf*s + 1)/(R*Lf*Cf*s^2 + Lf*s + R));

% Funcion de transferencia de la planta en tiempo continuo
[num_Gp,den_Gp] = ss2tf(A,B,C,D);
GpiL = zpk(tf(num_Gp,den_Gp));

figure; step(Vbase*Gpi)
legend('Respuesta de la FT Gpi')
title('Variable de salida: Corriente iL')

%% Eligiendo otra variable como salida
% Vector de estado: x = [vc  iL]'. Variable de salida: tensión sobre el capacitor del filtro LC
C=[1 0];

sistema_ssv=ss(A,B,C,D);
figure; step(Vbase*sistema_ssv)
legend('Respuesta del Modelo en SS con vc')

% Función de transferencia Y(s)/U(s) = C*inv[(sI - A)]*B
Gpv=zpk(R/(R*Lf*Cf*s^2 + Lf*s + R));

% Funcion de transferencia de la planta en tiempo continuo
[num_Gp,den_Gp] = ss2tf(A,B,C,D);
Gpvc = zpk(tf(num_Gp,den_Gp));

figure; step(Vbase*Gpv)
legend('Respuesta de la FT Gpv')
title('Variable de salida: Tensión vc')
close all

%% Obtención de las matrices en tiempo discreto
% z=tf('z',Tm);
% Seleccionamos una frecuencia de muestreo de 10 kHz
fm=20000;
Tm=1/fm;

[G,H]=c2d(A,B,Tm); C=[1 0]; D=0;

sisd_ssv=ss(G,H,C,D,Tm);  % Con ZOH
figure; step(Vbase*sisd_ssv,Vbase*sistema_ssv)
legend('Modelo Muestreado Sin Atraso. Salida vc','Modelo TC. Salida vc')

%% Se considera un tiempo de atraso Td igual a la mitad del periodo de muestreo T
Td=Tm/2;    
[eATD,H1] =c2d(A,B,Td);                      
[eATTD,H2]=c2d(A,B,(Tm-Td));
 
%% Matriz H0. La matriz identidad es 4 x 4 porque el orden de la planta es igual a 4, o sea, 4 estados.
Ho=inv(A)*eATTD*(eATD - eye(2))*B;          
 
%% Matriz H1
H1=inv(A)*(eATTD - eye(2))*B;
 
%% Matrices Gp y Hp de la planta en tiempo discreto, que tienen el cuenta el atraso de la implementación digital
Gp=[G Ho;zeros(1,3)];                                 
Hp=[H1;eye(1)];
Cca=[1 0 0];
sisd_ca_mp=ss(Gp,Hp,Cca,D,Tm);  % Con ZOH

figure; step(Vbase*sisd_ca_mp,Vbase*sistema_ssv)
legend('Modelo TD CA Tm/2. Salida vc','Modelo TC. Salida vc')

sisd_lc_mp=feedback(sisd_ca_mp,1);

%% Se considera un tiempo de atraso Td igual a un periodo de muestreo T
Td=Tm;    
[eATD,H1] =c2d(A,B,Td);                      
[eATTD,H2]=c2d(A,B,(Tm-Td));
 
%% Matriz H0. La matriz identidad es 4 x 4 porque el orden de la planta es igual a 4, o sea, 4 estados.
Ho=inv(A)*eATTD*(eATD - eye(2))*B;          
 
%% Matriz H1
H1=inv(A)*(eATTD - eye(2))*B;
 
%% Matrices Gp y Hp de la planta en tiempo discreto, que tienen el cuenta el atraso de la implementación digital
Gp=[G Ho;zeros(1,3)];                                 
Hp=[H1;eye(1)];
Cca=[1 0 0];
sisd_ca_up=ss(Gp,Hp,Cca,D,Tm);  % Con ZOH

figure; step(Vbase*sisd_ca_up,Vbase*sistema_ssv)
legend('Modelo TD CA Tm. Salida vc','Modelo TC. Salida vc')

sisd_lc_up=feedback(sisd_ca_up,1); sisc_lc=feedback(sistema_ssv,1);

figure; step(Vbase*sisc_lc,Vbase*sisd_lc_mp,Vbase*sisd_lc_up)
legend('Respuesta LC TC','Respuesta LC CA Tm/2','Respuesta LC CA Tm')

figure; bode(sisd_ca_mp,sisd_ca_up)
legend('Respuesta LA CA Tm/2','Respuesta LA CA Tm')

%% Lugar de raices
figure;
pzmap(sisd_ca_mp)
legend('LR CA Tm/2')

figure;
pzmap(sisd_ca_up)
legend('LR CA Tm')
close all

%% Solución en forma recursiva de la ecuación de espacio de estado en tiempo discreto
np=50;
t=0:Tm:np*Tm-Tm;
y=step(sisd_ssv,t)';

x=zeros(2,np);
xm=zeros(2,np);
iL=zeros(1,np);
vc=zeros(1,np);
td=zeros(1,np);
u=zeros(1,np);
r=zeros(1,np);

for k=1:np
    td(k)=(k-1)*Tm;
    
    r(k)=1;
    xm(:,k)=[x(1,k);x(2,k)];
    
    vc(k)=xm(1,k)/Vbase;
    iL(k)=xm(2,k)/Ibase;
    
    u(k)=Vbase*r(k);    
    x(:,k+1)= G*x(:,k) + H*u(k);
end

figure; stairs(t,y)
hold on
stairs(td,vc)

figure; stairs(t,y)
legend('Rta vc con step')

figure; stairs(td,vc)
legend('Rta vc con cálculo numérico')