clear
close all

s=tf('s');
% Planta
Gp=1.06/(s*(s+1)*(s+2));
[Z,P,Ks] = zpkdata(Gp,'v');
figure; rlocus(Gp)
Gcl=feedback(Gp,1);

pole_sc=pole(Gcl);

tf=20;
t=0:0.001:tf-0.001;
figure; step(Gcl);

r=t;
[y,t]=lsim(Gcl,r,t);
essv=100*(r-y');
figure; lsim(Gcl,r,t);
figure; plot(t,essv); grid
legend('essv sin compensar en %')
Kvsc=1/(essv(length(t))/100);

%% Primera seleccion
T=10; beta=10;

zc=(1/T); pc=(1/beta/T);
% Aporte de fase del compensador de atraso
s1d=-0.3307+1i*0.5864;

v1=s1d + pc;
v2=s1d + zc;

% Angulos del cero y del polo a s1d
fi1=angle(v1); fi1g=fi1*180/pi;
fi2=angle(v2); fi2g=fi2*180/pi;

Fase_comp=fi1g - fi2g;

% Condicion de Magnitud para calcular Kc'
Kca = abs(v1)/abs(v2);

%% Segunda seleccion
T=20; beta=10;

zc=(1/T); pc=(1/beta/T);
% Aporte de fase del compensador de atraso
s1d=-0.3307+1i*0.5864;

v1=s1d + pc;
v2=s1d + zc;

% Angulos del cero y del polo a s1d
fi1=angle(v1); fi1g=fi1*180/pi;
fi2=angle(v2); fi2g=fi2*180/pi;

Fase_comp=fi1g - fi2g;

v1p=s1d;
v2p=s1d + P(3);
v3p=s1d + P(2);

% Condicion de Magnitud para calcular Kc'
Kca = abs(v1)/(abs(v2));

Gc=Kca*(s + (1/T))/(s + (1/beta/T));

% Se calcula la ganancia estatica Kv'
Kvcc=Kvsc*Kca*beta;
essva=100*(1/Kvcc);

Gla=series(Gc,Gp);
figure; rlocus(Gp,Gla); axis([-2.5 1 -2 2]); axis('square')
legend('sin compensacion','con compensacion de atraso','Location','NorthEast')

Glcc=feedback(Gla,1);
polos_cc=pole(Glcc);

figure; pzmap(Glcc)
legend('Polos y Ceros de LC')

tf=100;
t=0:0.001:tf-0.001;

figure; step(Gcl,Glcc,t);
legend('sin compensacion','con compensacion de atraso','Location','NorthEast')

figure; bode(Gp,Gla)
legend('sin compensacion','con compensacion de atraso')

figure; rlocus(Gp,Gla)
legend('sin compensacion','con compensacion de atraso')
close all
%% 
tf=200;
t=0:0.001:tf-0.001;
r=t;
[y2,t]=lsim(Glcc,r,t);
essvc=100*(r-y2');
figure; lsim(Gcl,Glcc,r,t); axis([65 70 65 70]); grid
legend('sin compensacion','con compensacion de atraso','Location','NorthEast')

figure; plot(t,essvc); grid
legend('essv compensado en %')
Kvc=1/(essvc(length(t))/100);

%% Acción de control
tf=100;
t=0:0.001:tf-0.001;
rstep=ones(length(t),1);
[y3]=step(Glcc,t);
essp=rstep-y3;
figure; lsim(Gc,essp,t);

%% Diseño de los dispositivos electrónicos para la implementación
% Selecciono los capacitores C1 y C2
C1=47e-6;  C2=47e-6;
% Los resistores de la etapa de salida son iguales dado que Kc aproximadamente igual a 1
R3=10e3; R4=10e3;
% Resistores asociados a las constantes de tiempo tau1 y tau2
R1=T/C1;  R2=beta*T/C2;