%   DISEÑO DE UN COMPENSADOR EN ADELANTO POR EL
%     METODO DE RESPUESTA EN FRECUENCIA
                   

%   DADA LA SIGUIENTE PLANTA ( Ogata pag. 595)
%
% 	              4
%           G(s)= --------
%	           s(s+2)

%   Parametros de diseño:
%   Margen de Fase, PM=50º
%   Margen de Ganancia, GM=10 dB 
%   Ganancia Estatica Kv= 20 seg-1 

% 1.  Ajustamos la ganancia K para cumplir con el coeficiente estatico
%     de control requerido.
   
%    Kv=lim S Gc(s)G(s)=lim   S*K*4
%       s-o             s-o -------- = 20
%                            S(S+2) 
%   Se obtiene un K=10 con este K el sistema satisface el requisito de Kv
%   La planta queda:
%                               40
%                       G(s)= ---------
%                              S(S+2)
% Implementamos las transformaciones de radianes a grados

rad=180./pi;
grad=1./rad;
echo on 


PMD=50
GMD=10
num=4
den=[1 2 0]

% Ahora se dibuja bode del sistema sin compensar G(s).
% Con los valores de margen de fase, margen de ganancia, f.corte

figure(1)
bode(num,den)
[Gm,Pm,Wcg,Wcp]=margin(num,den);
margin(num,den)
pause

disp('El margen de ganacia sin compensar es:')
Gm
disp('El margen de Fase sin compensar es:')
Pm

% Segun el diagrama de bode el sistema sin compensar ni con ganancia K
% da un margen de Fase Pm=  51.8273
% Sacamos Bode con K=10
K=10

num2=40
den2=[1 2 0]
%
figure(2)
bode(num2,den2)
[Gm,Pm,Wcg,Wcp]=margin(num2,den2);
margin(num2,den2)
pause
%
% El nuevo margen de fase del sistema no compensado es
% Pm= 17.9642 (segun bode figura 2)

Pm2= 17.9642
Phides=50

% Con estos parametros diseñamos un compensador en adelanto
% Se halla el angulo necesario para compensar (+5º de seguridad)
% Phi=Phides-Pm2+5=   17.9642
 
Phi=(Phides-Pm2+5).*grad
pause
% Los 5 deg es el aumento que se hace para compensar el desplazamiento
% en la frecuencia de cruce de ganancia

% Calculamos el factor de atenuacion de alpha
%
%		1-sin(Phi)
%	alpha= ------------
%		1+sin(Phi)
%
alpha=(1-sin(Phi))./(1+sin(Phi))

% Obtenemos un Alpha de = 0.2482

pause

% Determinamos las frecuencias de corte
% Halllado el factor de atenuacion,
% por diagrama de bode "figura 2" hallamos wm.
%
% sea q=(1/sqrt(alpha)),
%
q=(1/sqrt(alpha));
%
% de donde se obtiene q= 2.0073(GANANCIA)
% con un valor de q= 2.0073 y utilizando bode Fig.2
%
dB=-20*log10(2.0073)
%
% de donde obtenemos dB= -6.0522
%
% por lo tanto:( De la gráfica de Bode, haciendo un zoom)

% OTRA FORMA DE HALLAR LA FRECUENCIA DE CRUCE:
% La frecuencia a la cual la magnitud del sistema
% no compensado es igual a -20*log10(1/sqrt(alpha))
[mag,phase,w]=bode(num2,den2);
i=find(mag<q+.1&mag>q) %El comando bode no da la mag en dB, ojo!!!
% ans=18
% w(18) = 4.1596
%
% luego  z= 4.1596
 
z=4.1596
p=z/alpha 
%
% sea:
%

%
%                                      s + 4.1596
% entoces el compensador será Gc(s) = -------------
%                                      s + 16.7594
%
% para compensar la atenuacion del compensador se
% se incrementa la ganancia en un factor (1/alpha)
%
Kc=K/alpha
%
% por esto Kc = 40.2909
%                                       s + 4.1596
% el compensador será Gc(s) = 40.2909 ---------------
%                                       s + 16.7594

%
% Con este podemos decir que la función de transferencia
% del sistema compensado será:
%
%                       s +4.1596          4
%  Gc(s)*G(s) = 40.2909 -------------- --------------
%                       s + 16.7594       s(s+2)   
% La funcion de trasferencia del compensador será:
%
num3=4.02909.*[1 4.1596]
den3=[1 16.7594]
%
% Hallamos bode "sistema compensado" 
%
num4=conv(num3,num2)
den4=conv(den3,den2)
%
figure(3)
bode(num4,den4)
[Gm,Pm,Wcg,Wcp]=margin(num4,den4);
margin(num4,den4)
pause
%