clear; close all;
format long

% Nuevos ejercicios año 2017
% Control del nivel de líquido en una cisterna
% Análisis transitorio y estacionario: modelo linealizado con variables
% siendo q1, h y q2

s=tf('s');

% Modelo de la válvula
Tv=5;                % Constante de tiempo de la válvula en segundos
Kv=2;                % Ganancia estática de la válvula
Gv=Kv/(s*Tv + 1);    % Función de transferencia de la planta

ro=1000;            % Densidad del agua en Kg/m3
g=9.8;              % Aceleración de la gravedad
A1=pi/4;            % Area entrada tanque en m2
A2=A1/100;          % Area orificio salida tanque en m2
He=1;               % Altura líquido en régimen estacionario
Qe=34.77;           % Caudal de líquido en régimen estacionario
qo=1;               % 1 Kg/s
t=0:0.01:400-0.01;

% Modelo de la cisterna o tanque: DH(s)/Dq1(s). 
% K2 está multiplicada x 100 para dar en cm
K2=100/(ro*A1); W=(A2^2/A1)*(ro*g/Qe);
Gt=K2/(s + W);
dq1=qo*ones(length(t),1);
figure; lsim(Gt,dq1,t)
legend('Variación de Altura en cm respecto a q1=1Kg/s')

figure; step(qo*Gt);
legend('Variación de Altura en cm respecto a q1=1Kg/s')

% Respuesta en frecuencia
figure; bode(qo*Gt)
legend('Función de Transferencia dh(s)/dq1(s)')

% Dinamica entre el fluido entrante y el saliente: Dq2(s)/Dq1(s)
Gq=W/(s + W);
figure; step(qo*Gq)
legend('Variación de caudal dq2(t) respecto a dq1(t)')

Qss=ro*A2*sqrt(2*g)*sqrt(He);
Hss=sqrt(Qe/(ro*A2*sqrt(2*g)));

% Función de transferencia de lazo abierto sin compensador
% Gla_act=series(Gt,Gv);
% figure; step(Gla_act)
% legend('Respuesta a Lazo Abierto con actuador')

% Respuesta en frecuencia
figure; bode(qo*Gq)
legend('Función de Transferencia dq2(s)/dq1(s)')

% Compensador PID
% Gc=Kp + (Ki/s) + (s*Kd);

%% Compensador proporcional
% Gc=Kp;
Gla=Gt; 
Gc1=0.5;
% Función de transferencia de lazo abierto con compensador
Gla1=series(Gc1,Gla);
% Función de transferencia de lazo cerrado con compensador
Glc1=feedback(Gla1,1);

Gc2=1;
Gla2=series(Gc2,Gla);
Glc2=feedback(Gla2,1);

Gc3=1.2;
Gla3=series(Gc3,Gla);
Glc3=feedback(Gla3,1);

Gc4=1.5;
Gla4=series(Gc4,Gla);
Glc4=feedback(Gla4,1);

Gc5=2;
Gla5=series(Gc5,Gla);
Glc5=feedback(Gla5,1);

Gc6=5;
Gla6=series(Gc6,Gla);
Glc6=feedback(Gla6,1);

% Respuesta de lazo cerrado
figure
step(Glc1,Glc2,Glc3,Glc4,Glc5,Glc6)
title('Respuesta a lazo cerrado')
legend('Kp=0,5','Kp=1','Kp=1,2','Kp=1,5','Kp=2','Kp=5')
%% --------------------------------------------------------------
% Compensador proporcional-integral
% Probar este conjunto de ganancias para observar los efectos:
Kp=2;             % Ganancia proporcional
Ki=0.05;         % Ganancia integral

% cero_PI=Ki/Kp;
% cero_PI=0.005714;
% Ki=Kp*cero_PI;

Gc=Kp + (Ki/s);
% Gcpi=0.010092*(s+0.005714)/s;

% Función de transferencia de lazo abierto
Glac=series(Gla,Gc);

% Función de transferencia de lazo cerrado
Glc=feedback(Glac,1);

[num,den]=tfdata(Glc,'v');

raices_num=roots(num)
raices_den=roots(den)

zpk(Glc)

polos_la=pole(Glac);
polos_lc=pole(Glc);

t=0:0.01:100-0.01;
figure; pzmap(Glc)
figure; step(3*Glc,t)
legend('Respuesta LC con PI')
[W,Z]=damp(Glc)
figure; bode(Glc)
%% --------------------------------------------------------------
% Compensador proporcional-derivativo
% Probar este conjunto de ganancias para observar los efectos:
Kp=6;             % Ganancia proporcional
Kd=4.75;            % Ganancia derivativa

Gc=Kp + s*Kd;

% Función de transferencia de lazo abierto
Glac=series(Gla,Gc);

% Función de transferencia de lazo cerrado
Glc=feedback(Glac,1);

[num,den]=tfdata(Glc,'v');

raices_num=roots(num)
raices_den=roots(den)

zpk(Glc)

polos_la=pole(Glac);
polos_lc=pole(Glc);

figure; pzmap(Glc)
figure; step(3*Glc)
legend('Respuesta LC con PD')

[W,Z]=damp(Glc)
figure; bode(Glc)
%% --------------------------------------------------------------
% Compensador PID
close all
clc
Gcpi=0.08*(s+0.18)/s;
                  
Gcpd=4.75*(s+1.263);
Gcpid=Gcpi*Gcpd;

% Reajustado
Gcpid=0.39478*(s+1.263)*(s+0.18)/s;

 % Función de transferencia de lazo abierto
Glac=series(Gla,Gcpid);

% Función de transferencia de lazo cerrado
Glc=feedback(Glac,1);

[num,den]=tfdata(Glc,'v');

raices_num=roots(num)
raices_den=roots(den)

zpk(Glc)

polos_la=pole(Glac);
polos_lc=pole(Glc);

figure; pzmap(Glc)
figure; step(3*Glc)
[W,Z]=damp(Glc)

legend('Respuesta LC con PID')
figure; bode(Glc)
figure; bode(Glac)
figure; rlocus(Glac)