%=================================================== 
%		DISEÑO EN EL PLANO W (BODE DISCRETO)
%===================================================
%Ing. Jorge Duque
%CUTB May/2001 Cartagena
%Curso de Control Digital


%EJEMPLO 4.12 (OGATA, SISTEMAS DE CONTROL EN TIEMPO DISCRETO)

%Considere la siguiente planta:

%			           k
%         G(s) = -------
%                 s(s+1)
% Diseñe un controlador digital, en el plano w, de tal forma que el margen de
% fase sea de 50°, elmargen de ganacia sea por lo menos 10dB, y la constante
% de error de velocidad estática Kv sea 2 (seg)-1
% El tiempo de muestreo es T= de 0.2 seg


%SOLUCION:


%1. SELECCION DEL PERIODO DE MUESTREO T:

% En este caso el período de muestreo está dado:

T=0.2

%Una regla práctica es muestrear con una frecuecnia de 10 veces el ancho de banda del 
%sistema en lazo cerrado
%
num=1
den=conv([1 0],[1 1]);
[numc,denc]=cloop(num,den)
bode(numc,denc)
% De la gráfica se determina que el ancho de banda es 1,5 rad/seg
% La frecuencia de muestreo debe ser como mínimo Ws=1,5*10= 15 rad/seg
% Con t=0.2 seg, Ws=2*pi/T = 31,7  rad/seg 
% se comprueba que el período de muestreo es satisfactorio


%2. OBTENCION DE LA FUNCION DE TRANSFERENCIA DISCRETA:
%Con un período de muestreo de 0.2 seg, la función de transferencia discreta, precedida
%de un retenedor de orden cero,se obtiene como sigue:

num=1;
den=conv([1 0],[1 1]);
printsys(num,den,'s')

[numd,dend]=c2dm(num,den,T,'zoh');
printsys(numd,dend,'z')

% numd/dend = 
 
%   0.018731 z + 0.017523
%  ------------------------
%  z^2 - 1.8187 z + 0.81873

%3. TRANSFORMACION DE LA FUNCION DE TRANSFERENCIA DISCRETA AL PLANO W :
% Ahora se transforma la función G(z) en G(w) Mediante la 
% transformación bilineal  o Tustin:

[nw,dw]=d2cm(numd,dend,T,'tustin');
printsys(nw,dw,'w')

%num/den (w):

 %  -0.00033201 w^2 - 0.096348 w + 0.99668
 % --------------------------------------
 %       w^2 + 0.99668 w 
 
 
 %4. DISEÑO DEL COMPENSADOR:
 %Se ensayará con un compensador en adelanto de la forma:
 
 %            1  +  tao*w
 %   Gd(w)= -----------------  donde  z=tao*w, p= alfa*tao*w  y 0 < alfa < 1
 %            1  +  alfa*tao*w
 
 % Ahora se compensa primero por la constante de error de velocidad, Kv=2
 %  Kv= lím w*Gd(w)*G(w)
 %      w->0
 %  
 %   		  	     w(1 + tao*w)(-0.00033201 w^2 - 0.096348 w + 0.996689)*K
 %  Kv=  lím     -----------------------------------------------------
 %       w->0   	w(1 + alfa*tao*w)(w + 0.99668) 	
 % luego Kv=~K
 
 K=2
 
 % El diagrama de bode sin compensar, pero compensado en ganancia:
 margin(K*nw,dw);
 
 % De la gráfica se observa que:
 % Pm=31,56° y Gm= 14.27 dB
 % Debido a que las especificaciones exigen un margen de fase de 50°, el ángulo 
 % adicional del adelanto de fase necesario para satisfacer este requisito
 % es de 20°.
 % Debido al desplazamiento en la frecuencia de cruce de ganancia se toma un
 % margen de seguridad de 8°, luego el ángulo de adelanto de fase requerido es de
 % 28°:
 %                   1-alfa
 %   Phi deseado= -----------
 %                   1+alfa
 
 phides=28*pi/180
 alfa=(1-sin(phides))/(1+sin(phides))
 
 %alfa = 0.3610
 
 % Ahora se encuentra la frecuencia donde la magnitud del sistema no compensado
 % es -20*log(1/sqrt(alfa)
 
 -20*log10(1/sqrt(alfa))
 
 %ans = -4.4245 dB

  
 % En el diagrama de Bode, por tanteo, se encuentra la frecuencia correspondiente
 % a -4.4.dB ( aprox.),la cual es de: aprox. 1.8 rad/seg
 
 % Luego la frecuencia de cruce de ganancia
 
 vc= 1.8
 
 %Otra forma de obtener la frecuencia de cruce:
 
 % -20*log10(1/sqrt(alfa)) en magnitud es sqrt(alfa)
 
 sqrt(alfa);%con este valor hallamos en el vector la frecuencia correspondiente
 % a  0.6009 (=-4.4245dB)
 [magw,phasew,ww]=bode(K*nw,dw,logspace(-1,1,1000));
 i=find(0.6<magw & magw<0.61)
 %i =

  %615
  % 616
  magw(615:616)
  ww(615:616)
  %ans =

  %  0.6065
  %  0.6017


%ans =

%   1.6952
%   1.7031
% luego la frecuencia de cruce es: 1.7 rad/seg
% se obtienen tao:
 
 tao= 1/(1.7*sqrt(alfa))
 
 % alfa*tao
 
 alfa*tao
 
 % ans =0.3338
 

 % Por tanto el compensador será:
 
 numcd=[tao 1]
 dencd=[alfa*tao 1]
 printsys(numcd,dencd,'w')
 
% num/den = 
 
%   0.97899 w + 1
%   -------------
%   0.35345 w + 1
 
% Gráfica de Bode del sistema compensado:

[nums,dens]=series(K*numcd,dencd,nw,dw)
margin(nums,dens)

%Hacemos la transformación en z del compensador

[ncz,dcz]=c2dm(K*numcd,dencd,T,'tustin');
printsys(ncz,dcz,'z')
[numt,dent]=series(ncz,dcz,numd,dend);
printsys(numt,dent,'z')
pause
disp('Para el lazo cerrado')
[ncl,dcl]=cloop(numt,dent,-1);
printsys(ncl,dcl)
disp('Observamos la respuesta del sistema compensado a una entrada paso')
figure(5)
t=[0:0.2:10];
dstep(ncl,dcl,t);

zoom
pause