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

%   DADA LA SIGUIENTE PLANTA ( Ogata pag.610)
%
% 	                1
%           G(s)= ---------------
%	           s(s+1)(0.5s+1)

%   Parametros de diseño:
%   Margen de Fase, PM=40º
%   Margen de Ganancia, GM=10 dB 
%   Ganancia Estatica Kv= 5 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*1
%       s-o             s-o --------------- = 5
%                            S(S+1)(0.5s+1) 
%   Se obtiene un K=5 con este K el sistema satisface el requisito de Kv
%   La planta queda:
%                                   5
%                       G(s)= ----------------
%                              S(S+1)(0.5s+1)
% Implementamos las transformaciones de radianes a grados

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

K=5

% Sacamos Bode con K=5

num=K
den=conv([1 1 0],[0.5 1])
%
bode(num,den)
[Gm,Pm,Wcg,Wcp]=margin(num,den);
margin(num,den)
pause
%
Pm
Gm
% El nuevo margen de fase del sistema no compensado es
%Pm = -12.9972

Phides=40

% Con estos parametros diseñamos un compensador en adelanto
% Se halla el angulo necesario para compensar (+5º de seguridad)

 
Phi=(-180+Phides+5)
%Phi =  -135

pause
% Los 5 deg es el aumento que se hace para compensar el desplazamiento
% en la frecuencia de cruce de ganancia


% Determinamos la nueva frecuencia de corte

[mag,phase,w]=bode(num,den);
i=find(phase<-132&phase>-135)


% FRERCUENCIA DE CRUCE wc:
% i=9
% w(9)     ===>   ans =  0.5337
% phase(9) ====>  ans = -133.0279

% Se elige el cero por debajo de esta frecuencia de cruce

z=0.1

% Se busca en BODE la atenuacion requerida para obtener 0 dB
% a esta frecuencia

mag(9)
%ans = 7.9863
%luego la atenuacion  es :
betha=8
% El polo:
p=z/betha
Kc=K/betha

                                  s + z
% Sea la red de atraso Gc(s) = Kc---------------
%                                 s + p

%
% Con este podemos decir que la función de transferencia
% del sistema compensado será:
%
%                       s + 0.1000          1
%  Gc(s)*G(s) =  0.6250 -------------- --------------
%                       s + 0.0125       s(s+2)   
% La funcion de trasferencia del compensador será:
%
num2=0.625.*[1 0.1]
den2=[1 0.0125]
%
% Hallamos bode "sistema compensado" 
%
num=1
den=conv([1 1 0],[0.5 1])
num3=conv(num2,num)
den3=conv(den2,den)
%
figure(3)
bode(num3,den3)
[Gm,Pm,Wcg,Wcp]=margin(num3,den3);
margin(num3,den3)
pause

%
%