% Definición de los parámetros de simulación
clc
clear
close all

fr=50;					% Frecuencia fundamental deseada de la tensión de salida en Hz
w=round(2*pi*fr);		% Frecuencia fundamental deseada de la tensión de salida en rad/seg
mf=100;          		% Número de pulsos PWM por periodo de la componente fundamental
fs=mf*fr;               % Frecuencia de conmutación (la "s" es de switching)
Ts=1/fs;				% Periodo de conmutación
nm=1;                   % Veces con la que se muestrea la variable controlada respecto del periodo de conmutación
fm=nm*fs;               % Frecuencia de muestreo
Tm=1/fm;                % Periodo de muestreo y de actualización de la acción de control.
fmi=fm;                 % Frecuencia de muestreo y de actualización del lazo con controlador por Modelo Interno
Tmi=1/fmi;              % Periodo de muestreo y actualización del lazo con controlador por Modelo Interno
ms=round(fm/fr);        % Número de muestras por periodo de la tensión de salida
mr=round(fmi/fr);       % Número de muestras del modelo interno por periodo de la tensión de salida
nc=200;				    % Número de ciclos o periodos de simulación
na=nc*ms;               % Número total de muestras por simulación
p=1024;                 % Número de puntos por muestra (define el paso de simulación)
np=na*p;                % Número total de puntos de la simulación
nr=mr*nc;               % Número total de muestras del modelo interno para el total de periodos de simulación
Vbase=311;              % Valor base de tensión para normalización
Ibase=50;               % Valor base de corriente para normalización
Vcc=400;                % Valor de la tensión CC de entrada al convertidor
N=mr;                   % Define el orden del controlador por Modelo Interno

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parametros de la Planta. Filtro LC y Carga
Lf=500e-6;              % Inductancia del filtro LC en uH
Cf=60e-6;          	    % Capacitancia del filtro LC en uF
Rc=Vbase/Ibase;         % Resistencia de carga 
wo=1/sqrt(Lf*Cf);       % Frecuencia angular de resonancia del filtro LC en rad/seg
fo=wo/2/pi;             % Frecuencia de resonancia del filtro LC en Hz
xita=1/(2*wo*Rc*Cf);      % Factor de amortiguamiento relativo

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parametros de la carga no lineal. Rectificador de puente completo no
% controlado con filtro capacitivo
Cdc=4700e-6;
Rdc=30;
Rs=0.1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=tf('s');
L=Lf; C=Cf; R=Rc;

% Función de transferencia de la planta en tiempo continuo
Gpc=(1/(L*C))/(s^2 + s/(R*C) + 1/(L*C)); 

z=tf('z',Tm);
Gpd=c2d(Gpc,Tm);

F=Gpd.num; D=Gpd.den;
n=F{:,:}; d=D{:,:};
b1=n(2); b2=n(3);
a1=d(2); a2=d(3);

% Parametros del controlador OSAP modificado
p1=a1; p2=a2; q1=b1; q2=b2;
p1m=(p2-a1*p1)/q1;
p2m=-p1*a2/q1;
q1m=1/q1;
q2m=(p1*b1-q2)/q1;
q3m=p1*b2/q1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Definición de las matrices del modelo para simulação sin la carga
% Las variables de estado son: x = [iL;vo]
Carga=0;
 
A=[ 0    -1/Lf
   1/Cf  Carga  ];
     
B=[1/Lf;0];

F=[0;-1/Cf];

C=[0 1];

D=0;

[Gc,Hc]=c2d(A,B,Tm/p);
[Gcc,Fc]=c2d(A,F,Tm/p);

% Matriz de Normalização  
Tn=[1/Ibase    0
       0    1/Vbase];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vectores de conmutacion y Matrices de decomposicion
V0=[0 0]';
V1=[1 0]';
V2=[0 1]';
V3=[1 1]';

M1= 1;
M2=-1;

M=[M1 M2];

SVector0=[0 0
          0 0];

SVectorf=[1 0
   	      0 1];

SVectors=[1 1
  	      1 1];

% Inicialización de las variables utilizadas en la simulación
S=1;
T1=p/2;
T2=p/2;
ACTR0=SVector0(1,:);
ACTRf=SVectorf(1,:);
ACTRs=SVectors(1,:);

kc=1:na*p;          % Tamaño del vector de las variables en tiempo continuo con periodo de muestreo Tm/p
ka=1:na;            % Tamaño de las variables en tiempo discreto con periodo de muestroe Tm

% Vectores de las variables de la planta en tiempo discreto
Xm=zeros(2,length(ka));
xe=zeros(3,length(ka));
vc=zeros(1,length(ka));
uk=zeros(1,length(ka));
upi=zeros(1,length(ka));
uc=zeros(1,length(ka));
e=zeros(1,length(ka));
y=zeros(1,length(ka));
r=zeros(1,length(ka));
td=zeros(1,length(ka));
iL=zeros(1,length(ka));
urp=zeros(1,length(ka));
Se=zeros(1,length(ka));
Id=zeros(1,length(ka));
Id1=zeros(1,length(ka));
uffc=zeros(1,length(ka));
dt=zeros(2,length(ka));
uosap_m=zeros(1,length(ka));
Vrms=zeros(1,length(ka));
vs=zeros(1,length(ka));

% Vectores de las variables de la planta en tiempo continuo
x=zeros(2,length(kc));
Io=zeros(1,length(kc));
Ior=zeros(1,length(kc));
upwm=zeros(length(kc),1);
PWM1=zeros(length(kc),1);
PWM2=zeros(length(kc),1);
io=zeros(length(kc),1);
t=zeros(length(kc),1);
portadora=zeros(length(kc),1);
Vdc=zeros(length(kc),1);
Vdc(1:100,1)=278*ones(100,1);
Vo=zeros(1,length(kc));
PWM=zeros(length(kc),2);

% Variables enteras o booleanas
Timer1 = 0;
flag_ADC = 1;
direccion = 1;
ComparReg1 = p/2;
ComparReg2 = p/2;

km=3;
i=0;
Tf=1/fr;
Npd=Tf/Tm;
Npc=Tf/(Tm/p);

MA=0;
CR=0;
flag_perturbacion=0;
kr=0.15;
Na=2;
Qr=1;

% Referencia de la tensión de salida generada para el control OSAP
tf=nc/fr;
tsin=0:Tm:tf;
ref=sin(2*pi*fr*tsin);
% figure; plot(tsin,ref);

% Inicia Barra de Progressão
perc=0.01;
passo=0.01;
wb = waitbar(0,'Aguarde, Simulando...');

% Se inicia la simulación del inversor: Loop externo simula al inversor con
% la planta en tiempo continuo, con paso Tm/p y el loop interno, simula el
% cálculo de la estrategia de control dentro de la computadora cada Tm
% segundos.
%% Loop externo
for k=2:np
    t(k)=k*Tm/p;
    Timer1=Timer1 + direccion;
    portadora(k)=Timer1;
    
    %% Loop interno: Muestreo y cálculo de la acción de control
    if flag_ADC == 1    % La CPU lee la bandera de interrupción del ADC
        flag_ADC = 0;   % pone a cero la bandera de IRQ del ADC
        km=km+1;
        i=i+1;
        td(km)=km*Tm;   % Vector de tiempo discreto
        
        Xm(:,km)=[x(1,k)/Ibase;x(2,k)/Vbase];
        iL(km)=Xm(1,km);
        vc(km)=Xm(2,km);
        vs(i)=vc(km);
        
        % Calcula el valor eficaz de la tensión muestrada
        Vrms(km)=sqrt(1/Npd*sum(vc(1:Npd).^2));
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Vector de estado
        xe(:,km)=[iL(km);vc(km);uk(km-1)];

        % Referencia de Tensión
%         r(km)=(311/Vbase)*sin(2*pi*fr*km*Tm);

%         if k<0.5*np
%             r(km)=(75/Vbase)*sin(2*pi*fr*km*Tm);
%         end

        % Salida de la planta
        y(km)=vc(km);

        % Erro de Tensao
        r(km)=ref(i);
        r(km+1) = ref(i+1);
        e(km)=r(km) - y(km);

        % Calculo del controlador repetitivo con acción integral (Modelo Interno)
%         if km>N
%             Se(km)=e(km)+Se(km-N);
%             urp(km)=kr*Se(km-N+Na);
%         end

        % Calculo del Controlador Plug-in Repetitivo (Modelo Interno)
        if km>N
           urp(km)=kr*e(km - N + Na) + Qr*urp(km - N);
        end

        % Cálculo del controlador OSAP 
        uosap_m(km)=q1m*r(km+1) + p1m*y(km-1) + p2m*y(km-2) + q2m*uosap_m(km-1) + q3m*uosap_m(km-2);
        uk(km)=uosap_m(km) + urp(km);
        
    if MA==1
       uk(km)=r(km);
    end
    
    % Se implementa la modulación vectorial "Space Vector"
    % Identificación del sector en el cual se encuentra uk
    if uk(km)>0
        S=1;
    else
        S=2;
    end

    % Cálculo de los tiempos (ciclos útiles) de aplicación en cada sector T1 y T2
    T1=p*[M((S-1)+1)]*uk(km)*(Vbase/Vcc);
    T2=p-T1;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % MAQUINA DE ESTADOS:
    ACTR0=SVector0(S,:);
    ACTRf=SVectorf(S,:);
    ACTRs=SVectors(S,:);
    end
    
    % Se implementa el contador o timer del módulo PWM
    if Timer1==p/2
        direccion=-1;    
    end

    if Timer1==0
        direccion=1;
        flag_ADC=1;      % Se activa la interrupción para el muestreo
    end

    if Timer1 < T2/4
        PWM(k,:)=ACTR0;
    else
        if Timer1 < (T2/4 + T1/2)
            PWM(k,:)=ACTRf;
        else
            if Timer1 < (T2/2 + T1/2)
                PWM(k,:)=ACTRs;
            end
        end
    end
                
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Fin de la modulación Space Vector
    % Obtención de la tensión PWM aplicada al filtro LC
    upwm(k)=(PWM(k,1)-PWM(k,2))*Vcc;   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Calcula la perturbación de carga lineal resistiva
    if flag_perturbacion == 1
        Ior(k)=x(2,k)/Rc/10;      % Condición de vacio
        
        if t(k) >= 2.005  %0.085
            Ior(k)=x(2,k)/Rc*2.0; % Produce un escalón de carga
        end
        
        if t(k) >= 2.065  %0.145
            Ior(k)=x(2,k)/Rc/10;  % Vuelve a la condición de vacio
        end
    else 
        Ior(k)=x(2,k)/Rc;         % Condición de carga nominal
    end
    
    %% Simula un rectificador no controlado a diodos en puente completo con filtro capacitivo
    if abs(x(2,k))>=Vdc(k)
        Io(k)=sign(x(2,k))*(abs(x(2,k))-Vdc(k))/Rs;
        Vdc(k+1)= Vdc(k) + Tm/p/Cdc*(( abs(x(2,k)) - Vdc(k) )/Rs - Vdc(k)/Rdc );
    else
        Io(k)=0;
        Vdc(k+1)= Vdc(k) + Tm/p/Cdc*(-Vdc(k)/Rdc );
    end
    
    % Si se cumple la condición, se aplica a la planta una perturbación de
    % carga lineal, en caso contrario, se aplica la carga no lineal
    if CR == 1
        Io(k) = Ior(k);
    end
    
    % Calculo del valor eficaz de la tensión de salida
    Vo(k)=x(2,k);
    
    x(:,k+1)= Gc*x(:,k) + Hc*upwm(k) + Fc*Io(k);
    
    % Atualiza a barra de progressão
    if fix((round(np))*perc)==k
        waitbar(perc)
        perc=perc+passo;
    end
end
close(wb)

% figure
% plot(t,Switp(:,1),t,portadora/p*2,t,Switp(:,2)-1.2,t,(upwm/Vcc)+2.2)
% axis([0 0.02 -1.2 1.2])
 
figure
plot(t,upwm)
axis([3.88 max(t) -Vcc-10 Vcc+10])
 
figure
plot(td(1:km),e(1:km))
legend('error')
axis([0 max(t) -1.0 1.0])

figure
plot(td,uk)
legend('Acción de Control')
axis([3.88 max(td) -1.2 1.2])

figure
plot(td,vc,'r')
hold
plot(td,r(1:km),'m')
legend('vo','ref')
axis([3.88 max(td) -1.1 1.1])

figure
plot(t,x(2,1:k),'r')
hold
plot(t,Io)
grid
legend('vo','io')
axis([3.88 max(t) -320 320])

if flag_perturbacion == 1
    figure
    plot(t,x(2,1:k),'r')
    hold
    plot(t,Io)
    grid
    legend('vo','io')
    axis([1.85 2.25 -320 600])
end


THD_tension_salida
RMSvo=rms(Vo)