%Este programa simula un contolador   PWM   para   el
%convertidor Buck (seguimiento). Emplea los programas
%auxiliares buck0.m, buck1.m, bucku.m y ubuck.m

clear;
%
global R L E ru i eta j z Phi1 Phi2 Alpha xd uu Z0 T Mu Mua
R=2.8e-2;L=1e-5;E=126;T=1/12000; eta=0.7;
Tfin=.006;x0=0;Alpha=0.3;

%Calculo de los parámetros
%
Theta1=R/L;
Theta2=E/L;
Phi1=exp(-Theta1*T);
Phi2=Theta2/Theta1;

%Calculo de xd(kt), funcion deseada
%
nn=Tfin/T;
xd=zeros(nn+1,1);
na=nn/6;
for i=0:na,
	xd(i+1)=1237*i/(na);
end
for i=na:2*na,
	xd(i+1)=1237;
end
for i=2*na:3*na,
	xd(i+1)=1237-1237/(na)*(i-2*na);
end
for i=3*na:4*na,
	xd(i+1)=-1237/(na)*(i-3*na);
end
for i=4*na:5*na,
	xd(i+1)=-1237;
end
for i=5*na:6*na,
	xd(i+1)=-1237+1237/(na)*(i-5*na);
end
xd=-xd;
%Inicialización del reloj y de los vectores de tiempo, 
%estado y control (discretos y continuos)
%
t=zeros(2*nn+1,1);
tiempo=zeros(nn+1,1);
x=zeros(2*nn+1,1);
z=zeros(nn+1,1);
uu=zeros(nn+1,1);
Mu=zeros(2*nn+1,1);
ru=zeros(2*nn+1,1);
X=[];
U=[];
T1=[];
MU=[];
T2=[];
reloj=0;
runT=0;
x(1)=x0;
t(1)=reloj;
%ue0=1/2*x0*Theta1/Theta2;
Mu(1)=-0.2748; Mua=Mu(1);
uu(1)=Mu(1);
i=0;

%Solucion de las ecuaciones del PWM
%
for j=1:nn-1,
	tiempo(j)=tiempo(j)+j*T;

%Primer ciclo del PWM, U=1 o U=-1
%
	i=i+1;
	MuT=abs(Mu(i))*T;
	reloj=reloj+MuT;
	t(i)=reloj;
	if MuT==0
		x(i+1)=x(i);
	else
		zc0=x(i);
		[tpri,zc]=ode23('buck1',reloj-MuT,reloj,zc0,0.001);
		n=length(tpri);
		x(i+1)=zc(n);
		uax=ones(n,1)*sign(Mu(i));
		X=[X;zc];
		T1=[T1;tpri];
		U=[U;uax];
		Mu(i+1)=Mu(i);
	end

%Segundo ciclo del PWM, U=0
%
	i=i+1;
	resT=T-abs(Mu(i))*T;
	reloj=reloj+resT;
	t(i)=reloj;
	if resT==0
		x(i+1)=x(i);
	else
		zc0=x(i);
		dt=resT/10;
		[tpri,zc]=ode23('buck0',reloj-resT,reloj,zc0,0.001);
		n=length(tpri);
		x(i+1)=zc(n);
		uax=zeros(n,1);
		X=[X;zc];
		T1=[T1;tpri];
		U=[U;uax];
	end

% Obtener la función promedio
%
	z(j)=1/2*(x(i)+x(i-1));

%Cálculo de la razón de trabajo un paso adelante
%
	M1=fsolve('ubuck',0);
	uu(j+1)=M1;
	
%Función limitadora de la razón de trabajo
%
  if xd(j+1)>0,
 	if uu(j+1)>1, Mu(i+1)=1;
	elseif uu(j+1)>0, Mu(i+1)=uu(j+1);
	else Mu(i+1)=0;
	end
	Mua=Mu(i+1)
  else
 	if uu(j+1)>0, Mu(i+1)=0;
	elseif uu(j+1)>-1, Mu(i+1)=uu(j+1);
	else Mu(i+1)=-1;
	end
	Mua=Mu(i+1)
  end
end  % Fin del ciclo FOR

for h=1:length(Mu)
	if Mu(h)==0, Mu(h)=Mu(h-1);
	end
end
%Graficas
%
clg
subplot(221),plot(T1,X)
      title(' inductor current trajectory')
%      axis([0 (reloj-T) 0 max(z)+300])
subplot(222),plot(t(1:2*nn-2,1),Mu(1:2*nn-2,1))
      title('duty ratio input')
%     axis([0 Tfin -0.01 0.40])
subplot(223), plot(T1,U)
      title('input(switch) sequence')
	axis([0 Tfin -1.10 1.10]),
%      axis([0 Tfin -0.10 1.10]),
subplot(224),plot(tiempo(1:nn-1,1),xd(1:nn-1,1))
      title('estimated trajectory')
%      axis([0 Tfin 0 max(xd)+300]),

pause;
figure;
clg
subplot(221),plot(tiempo(1:nn-1,1),xd(1:nn-1,1))
      title('estimated trajectory')
%      axis([0 Tfin 0 max(xd)+300]),
subplot(222),plot(tiempo(1:nn-1,1),z(1:nn-1,1))
      title('average value of the inductor current trajectory')
%      axis([0 Tfin 0 max(xd)+300]),
subplot(223),plot(t(1:2*nn-2,1),ru(1:2*nn-2,1))
      title('external perturbation signal')
%     axis([0 Tfin -50 50])
end


