%Este programa simula un contolador  PWM  para  el
%convertidor  Boost  (estabilización). Emplea  los
%programas   auxiliares   Raboost0.m,  Raboost1.m, 
%Rauboost.m y Ramboost.m

clear;
global R L E C rf cf Yave Alpha Phi1 Phi2 x i Ym ru eta m1 m2 T
R=50; L=20e-3; E=126; C=0.1e-6; Yave=3; T=1/8e3;
Tfin=0.002; Alpha=0.3; eta=0.75;rf=10;cf=20e-6;
x0=[0.78 -31 -31];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculo de los puntos de equilibrio
Ueq=0.2;
Ta1=1/L;
Ta2=1/C;
Ta3=1/(R*C);
Ta4=E/L;
X1eq=(Ta3*Ta4*Ueq)/(Ta2*Ta1*((1-Ueq)^2));
X2eq=-(Ta4*Ueq)/(Ta1*(1-Ueq));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Calculo de las raices
%
Coef=[1;(1/(R*C));(1/(L*C))];
Ra=roots(Coef);
m1=Ra(1);
m2=Ra(2);

%Calculo de los parámetros
%
Phi1=exp(m1*T);
Phi2=exp(m2*T);

%Calculo de U en estado estacionario
%
options(5)=1;
Ue=fsolve('ssubb',0.5,options);

%Calculo de y(kt) en estado estacionario
%
va1=Phi1^(1-Ue);
va2=Phi2^(1-Ue);
Den=m2*(1-va1)-m1*(1-va2);
Nun1=m2*va1-m1*va2;
Nun2=va2-va1;
Ym=(Nun1/Den)*((E*Ue*T)/L)+(Nun2/Den)*(E/L);

%Inicialización del reloj y de los vectores de tiempo, 
%estado y control (discretos y continuos)
%
nn=Tfin/T;
t=zeros(2*nn+1,1);
x=zeros(2*nn+1,3);
u=ones(2*nn,1);
ru=zeros(2*nn+1,1);
X=[];
U=[];
T1=[];
MU=[];
reloj=0;
i=1;
x(i,:)=x0;
t(i)=reloj;
ru(i)=0;

%Solucion de las ecuaciones del PWM
%
for j=1:nn,

%Cálculo de la razón de trabajo
%
options(5)=1; 
ui=fsolve('ssmbb',0.5,options);
uu=ui;

%Función limitadora de la razón de trabajo
if uu>1, Mu=1;
elseif uu>0, Mu=uu;
else Mu=0;
end

%Primer ciclo del PWM, U=1
%
MuT=Mu*T;
reloj=reloj+MuT;
dt=MuT/10;
if MuT==0
x(i+1,:)=x(i,:);
else
disp(reloj)
z0=x(i,:);
[tpri,z]=ode23('ssbbf1',reloj-MuT,reloj,z0,0.001);
n=length(tpri);
x(i+1,:)=z(n,:);
uax=ones(n,1);
Muax=Mu*ones(n,1);
X=[X;z];
T1=[T1;tpri];
U=[U;uax];
MU=[MU;Muax];
end

u(i)=1;
i=i+1;
t(i)=reloj;
ru(i)=ru(i-1);

%Segundo ciclo del PWM, U=0
%
resT=T-Mu*T;
reloj=reloj+resT;
if resT==0
x(i+1,:)=x(i,:);
else
z0=x(i,:);
%dt=resT/10;
[tpri,z]=ode23('ssbbf0',reloj-resT,reloj,z0,0.001);
n=length(tpri);
x(i+1,:)=z(n,:);
uax=zeros(n,1);
Muax=Mu*ones(n,1);
X=[X;z];
T1=[T1;tpri];
U=[U;uax];
MU=[MU;Muax];
end

u(i)=0;
i=i+1;
t(i)=reloj;
u(i)=u(i-1);
ru(i)=ru(i-1);
end

%Graficas
%
clg

YAVE=Yave*ones(length(T1),1);

subplot(221),
plot(T1,X(:,1),'y-',T1,YAVE,'r-.')
%title('input inductor current')
title('corriente en el inductor')
xlabel('tiempo[s]')
ylabel('Amperios')
axis([0 0.002 0 4])

subplot(222),
plot(T1,X(:,2),'y',T1,X(:,3),'r')
%title('capacitors voltage output')
title('salida de voltajes en capacitores')
xlabel('tiempo[s]')
ylabel('Voltios')
axis([0 0.002 -200 50])


subplot(223),
plot(T1,MU)
%title('duty ratio')
title('razón de trabajo')
xlabel('tiempo[s]')
axis([0 0.002 0.2 1.2])

%subplot(224),
%plot(T1,U)
%title('input (switch) sequence')
%axis([0 0.002 -0.2 1.2]),

%pause;
%clg

%subplot(221),
subplot(224),plot(t,ru)
%title('perturbation noise')
title('perturbación estocástica externa')
axis([0 0.002 -50 50]),

end


