% DISEÑO DE CONTROLADORES PI EN EL DOMINIO DE LA FRECUENCIA - Tema 10 %
clear all
close all
clc

%% Caso a: la K del controlador no está fijada por las especificaciones

num=[1];
den=conv([1 1 0],[1 10])

G=tf(num,den)
[Gm,Pm,Wgm,Wpm] = margin(G)

figure(1)
bode(G)
grid on

%Especificaciones: ev<1% y SO<20
SO=0.2; ev=0.01;

d = sqrt(log(SO)^2/(pi^2+log(SO)^2))
alpha=acos(d)

% wn = (pi-alpha)/(ts*sqrt(1-d^2))

raiz_doble=sqrt(sqrt(1+4*d^4)-2*d^2);

Mf_rad=atan((2*d)/raiz_doble)  % approximado Mf=70-SO
Mf=Mf_rad*180/pi

% Kv es limite cuando s->0 de sC(s)G(s)
% ev=0.01 ---> 1/Kv=0.01 ---> Kv=Kc*Kg y Kv=1/0.01

Kg=0.1;
Kv=1/0.01;
Kc=Kv/Kg

Kcdb=real2deb(Kc)

G_P=Kc*G
[Gm,Pm,Wgm,Wpm] = margin(G_P)

figure(1)
hold on
bode(G_P)

Delta_Mf=Mf-Pm+10  % >90 grados: un PD no es suficiente

% Entonces usamos un PI

% Primero elegimos Kc para garantizar el margen de fase deseado

Phase_wc=180-Mf

% en el bode buscamos la frecuencia
figure(2)
bode(G)
grid on

wc=0.77

% calculamos Kc para que esa se la frecuencia de corte

syms s

Gw= 1/(s^3 + 11 * s^2 + 10 * s)

s=i*wc;
Gw_wc=eval(Gw)

Kc=1/abs(Gw_wc)
Kcdb=real2deb(Kc)

G_P=Kc*G
[Gm,Pm,Wgm,Wpm] = margin(G_P)

% Ti se elige de manera que el margen de fase no decrezca: 1/Ti < wc/10 -->
% Ti > 10/wc

Ti=(10/wc)*1.5

C=tf(Kc*[Ti 1],[Ti 0])

G_PI=C*G
[Gm,Pm,Wgm,Wpm] = margin(G_PI)

figure(2)
hold on
bode(G_PI)

Glc=feedback(G_PI,1)
figure(3)
step(Glc)
datos = stepinfo(Glc,'RiseTimeLimits',[0 1]);

fprintf('ts = %.4f \nSO = %.4f \n', datos.RiseTime, datos.Overshoot)

% para ajustar tomamos un Mf un poco más alto

Mf=55

Phase_wc=180-Mf

% en el bode buscamos la frecuencia
figure(4)
bode(G)
grid on

wc=0.6 % desde el bode

% calculamos Kc para que esa se la frecuencia de corte

syms s

Gw= 1/(s^3 + 11 * s^2 + 10 * s)

s=i*wc;
Gw_wc=eval(Gw)

Kc=1/abs(Gw_wc)
Kcdb=real2deb(Kc)

G_P2=Kc*G
[Gm,Pm,Wgm,Wpm] = margin(G_P2)

% Ti se elige de manera que el margen de fase no decrezca: 1/Ti < wc/10 -->
% Ti > 10/wc

Ti=(10/wc)*1.5

C=tf(Kc*[Ti 1],[Ti 0])

G_PI2=C*G
[Gm,Pm,Wgm,Wpm] = margin(G_PI)

figure(4)
hold on
bode(G_PI2)

Glc2=feedback(G_PI2,1)
figure(3)
hold on
step(Glc2)
datos = stepinfo(Glc2,'RiseTimeLimits',[0 1]);

fprintf('ts = %.4f \nSO = %.4f \n', datos.RiseTime, datos.Overshoot)



%% Caso b: la K del controlador está fijada por las especificaciones
clear all
close all
clc
% Paso 1: diseñar K=Kc/Ti para satisfacer las especificaciones
% Paso 2: diseñar un controlador C(s)=(1+Ti*s) para el sistema
% G'(s)=K/s*G(s)

num=[1];
den=conv([1 1 0],[1 10])

G=tf(num,den)
[Gm,Pm,Wgm,Wpm] = margin(G)

figure(1)
bode(G)
grid on

I=tf(1,[1 0])

%Especificaciones: ea<1% y SO<20
SO=0.2; ea=0.01;

% Ka es limite cuando s->0 de s^2*C(s)G(s)
% ea=0.01 ---> 1/(Ka)=0.01 ---> Ka=K*Kg y Kg es la ganancia estatica del
% sistema compensado con 1/s

Ka=1/0.01

Kg=0.1 

K=Ka/Kg

Gp=K*I*G
[Gm,Pm,Wgm,Wpm] = margin(Gp)

figure(1)
hold on
bode(Gp)

% Pm es demasioado pequeño, un cero aporta 90 grados con lo cual no vamos a
% poder satisfacer las especificaciones

clear all
close all
clc
% Paso 1: diseñar K=Kc/Ti para satisfacer las especificaciones
% Paso 2: diseñar un controlador C(s)=(1+Ti*s) para el sistema
% G'(s)=K/s*G(s)

num=[10];
den=conv([1 1],[1 10])

G=tf(num,den)
[Gm,Pm,Wgm,Wpm] = margin(G)

figure(1)
bode(G)
grid on

I=tf(1,[1 0])

%Especificaciones: ev<1% y SO<60
SO=0.6; ev=0.01;

% Kv es limite cuando s->0 de s*C(s)G(s)
% ev=0.01 ---> 1/(Kv)=0.01 ---> Kv=K*Kg y Kg es la ganancia estatica del
% sistema compensado con 1/s

Kv=1/0.01

Kg=1 

K=Kv/Kg

Gp=K*I*G
[Gm,Pm,Wgm,Wpm] = margin(Gp)

figure(1)
hold on
bode(Gp)

d = sqrt(log(SO)^2/(pi^2+log(SO)^2))
alpha=acos(d)

% wn = (pi-alpha)/(ts*sqrt(1-d^2))

raiz_doble=sqrt(sqrt(1+4*d^4)-2*d^2);

Mf_rad=atan((2*d)/raiz_doble)  % approximado Mf=70-SO
Mf=Mf_rad*180/pi

% El cero se agrega como en un PD
Delta_Mf=Mf-Pm+10
    
Delta_Mf_rad=(Delta_Mf*pi)/180
Ti=tan(Delta_Mf_rad)/Wpm

% Entonces K=Kc/Ti --> Kc=K*Ti

Kc=K*Ti

C=tf(Kc*[Ti 1],[Ti 0])

G_PI=C*G
[Gm,Pm,Wgm,Wpm] = margin(G_PI)

figure(1)
hold on
bode(G_PI)

Glc=feedback(G_PI,1)
figure(3)
step(Glc)
datos = stepinfo(Glc,'RiseTimeLimits',[0 1]);

fprintf('ts = %.4f \nSO = %.4f \n', datos.RiseTime, datos.Overshoot)

% para ajustar tomamos un Mf un poco más alto

%% funciones que necesito
function b=real2deb(a)

b=20*log10(a);
end

function b=deb2real(a)

b= 10^(a/20);

end