%%%%%%%%%%% Comprehensive Problem 8.4 %%%%%%%%%%%
%   Discrete-Time Control Problems using        %
%       MATLAB and the Control System Toolbox   %
%   by J.H. Chow, D.K. Frederick, & N.W. Chbat  %
%         Brooks/Cole Publishing Company        %
%                September 2002                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hydro-turbine System

clear
disp('Comprehensive Problem 8.4')
disp('PID control of hydro-turbine system')

dhydro  % read data
disp('discrete-time transfer function')
Gz = tf(dnum,dden,Ts)  % form tf object
disp('*****>'), pause

disp(' ')
disp('Hydro-turbine system with PID controller')
disp('                     Ts * z          z-1    ')
disp('Gc(z) = K_P (1 + K_I ------ + K_D * ------ )')
disp('                       z-1          Ts * z  ')

disp('In several previous Comprehensive Problems, we have shown that ')
disp('good control design requires K_P = 0.1 - 0.6, K_I = 0.5 - 2,') 
disp('and K_D = 0.1 - 0.6. Thus for this design, we pick K_P = 0.4')
disp('and K_I =1 and vary K_D to satisfy the desired performance specifications.')

disp('*****>'), pause
disp('Derivative gain') 
K_D = [0.1:0.1:0.6]
len_KD = length(K_D);
disp('*****>'), pause

km_mat = zeros(len_KD,1);
pm_mat = km_mat;
Mo_mat = km_mat;
tr_mat = km_mat;
ts_mat = km_mat;
p_mat  = [];
d_mat  = [];
y_mat  = [];
t = [0:Ts:20]';
K_P = 0.4; K_I = 1; 

 for jj = 1:len_KD
  G_c = K_P*(1 + tf(K_I*Ts*[1 0],[1 -1],Ts) ...
               + tf(K_D(jj)*[1 -1],Ts*[1 0],Ts) );
  [km,pm,gw,pw] = margin(Gz*G_c);
  G_cl = feedback(Gz*G_c,1);
  y = step(G_cl,t);
  [wn,d] = damp(G_cl);
  p_mat = [p_mat pole(G_cl)];
  d_mat = [d_mat d];
  y_mat = [y_mat y];
  [Mo,tp,tr,ts,ess] = tstats(t,y,1);
  if isempty(Mo) 
      Mo = 0;
  end
  if isempty(ts) 
      disp('Settling time set to max time of simulation')
      ts = max(t);
  end
  if isempty(tr) 
      tr = max(t);
  end  
  km_mat(jj) = km;
  pm_mat(jj) = pm;
  Mo_mat(jj) = Mo;
  tr_mat(jj) = tr;
  ts_mat(jj) = ts;
 end

disp('Performance values')
disp('    K_D      Mo         tr       ts')
disp([K_D' Mo_mat tr_mat ts_mat])
disp('*****>'), pause

disp('    K_D          damping ratios of the 4 modes')
disp([K_D; d_mat]')
disp('*****>'), pause
 
figure, plot(t,y_mat), grid
title('Step response for K_D = 0.1,0.2,...,0.6')
legend('0.1','0.2','0.3','0.4','0.5','0.6')
xlabel('Time (s)')
ylabel('Amplitude')
disp('Plotting step responses')

disp('*****>'), pause
figure, plot(p_mat,'x'), zgrid(0,[]), grid, axis('equal')
title('Root locus plot of closed-loop system vs K_D')
xlabel('Real part')
ylabel('Imag part')
disp('Plotting root locus') 
disp('*****>'), pause

disp('A PID control with K_P = 0.4, K_I = 1, & K_D = 0.3')
disp('will satisfy the design specifications.')

disp('end of Comprehensive Problem 8.4')
%%%%%%%%%%