% distred.m -illustration of disturbance attenuation
% KJA, 2012

aminit;					% initialize FBS macros

%   Process model
b = 2; a = 1;

%   Open Loop
P = tf(b, [1 a]);
tmax = 8; t = 0:0.02:tmax;
yol = step(P, t);

figure(1); subplot(321);
pl = plot(t, yol, 'b-.');
set(pl,  'LineWidth',  AM_data_linewidth);
hold on;

%   Proportional control
Str = {'r:', 'r--', 'r-'}; jjj = 1;
for kp = [0.5 1 2]
  Gyv = tf(b, [1 a+b*kp]);
  Guv = tf(b*kp, [1 a+b*kp]);
  y = step(Gyv, t);
  u = -step(Guv, t);
  str = Str{jjj}; jjj = jjj+1;

  subplot(321); 
  pl = plot(t, y, str); hold on;
  set(pl,  'LineWidth',  AM_data_linewidth);
  amaxis([0 tmax -0.1 2.1]);
  ylabel('y');
% title('P')

  subplot(323); 
  pl = plot(t, u, str); hold on;
  set(pl,  'LineWidth',  AM_data_linewidth);
  amaxis([0 tmax -1.2 0]);
  xlabel('t'); ylabel('u');
end

amprint('distred-P.eps')

%   PI Control
figure(2); subplot(321);
pl = plot(t, yol, 'b-.');
set(pl,  'LineWidth',  AM_data_linewidth);

Gains = [kp; 0];			% Design parameters
w0 = 1; z = 0.707; wc = a/2/z;		% Controller parameters
jjj = 1;

for w0 = [wc 1 2]
  str = Str{jjj}; jjj = jjj+1;
  kp = (2*z*w0-a)/b;
  ki = w0^2/b;
  Gains = [Gains [kp; ki]];
  den = [1 a+b*kp b*ki];
  Gyv = tf([b 0], den);
  Guv = -tf([b*kp b*ki], den);
  Gupv = -tf([b*kp b*ki 0], den);
  y = step(Gyv, t);
  u = step(Guv, t);
  up = step(Gupv, t);

  subplot(321); 
  pl = plot(t, y, str); hold on;
  set(pl,  'LineWidth',  AM_data_linewidth);
  amaxis([0 tmax -0.1 2.1]);
  ylabel('y');
% title('PI');

  subplot(323); 
  pl = plot(t, u, str); hold on;
  set(pl,  'LineWidth',  AM_data_linewidth);
  amaxis([0 tmax -1.2 0]);
  ylabel('u'); xlabel('t');
end

amprint('distred-PI.eps')




