%%%%%%%%%%% Comprehensive Problem 8.3 %%%%%%%%%%%
%   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                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Electric Power Generation System

clear all, close all
disp('Comprehensive Problem CP8_3')

dpower2    % read data
disp('*****>'), pause
disp('Power system model state matrices, with Vterminal, ')
disp('  w(speed), and P(power) as output')
disp('Continuous-time model')
Gs
disp('*****>'), pause

disp('discretized model with a sampling period of 0.01 s')
Ts = 0.01
Gz = c2d(Gs,Ts) 
disp('Note that most of the zero entries of the A & B matrices')
disp('  of the continuous time model have become nonzero in the')
disp('  discrete time models. The C & D matrcies are unchanged.')
disp('*****>'), pause

disp('System model with Vterminal as the only output')
G1 = Gz(1,1)
disp('*****>'), pause

%%%%%%%
disp('Part a')
disp('Applying P controller to Vterminal control')
disp('Gain sweep on proportional gain K_P')
K_P = [24:4:40]
len_KP = length(K_P);
disp('*****>'), pause

disp('Sweeping of gains...')

t = [0:0.01:15]';
yall = [];
disp('Prop gain    DC gain   Settling time (s)')
for jj = 1:len_KP
  G1cl = feedback(G1*K_P(jj),1);
  y = step(G1cl,t);
  y = 0.05*y;
  [Mo,tp,tr,ts,ess] = kstats(t,y,0.05*dcgain(G1cl));
  disp([K_P(jj) dcgain(G1cl) ts])
  yall = [yall y];
end
figure, plot(t,yall,'-'),grid
xlabel('Time (s)'), ylabel('Amplitude')
title('Sweep of Proportional Gain')
text(3, 0.04, 'KP = 24')
text(3, 0.055, 'KP = 40')
disp('*****>'), pause
hold off

disp('A proportional gain of K_P = 28 will result in')
disp('a dcgain of lower than 0.95 and a settling time')
disp('of less than 5 s.')

%%%%%%%
disp('Part b')
disp('*****>'), pause
disp('Applying PI controller to Vterminal control')
disp('PI controller is of form')
disp('             (1+K_I T_s) z - 1  ')
disp('K_V(z) = K_P(-----------------) ')
disp('                    z-1         ')

disp('Varying K_I and K_P')
disp('Integral gain K_I')
K_I = [0.05:0.05:0.5];
len_KI = length(K_I);
disp('*****>'), pause
disp('Proportional gain') 
K_P = [10:5:40]
len_KP = length(K_P);
disp('*****>'), pause
disp('Sweeping of gains...')

km_mat = zeros(len_KP,len_KI);
pm_mat = km_mat;
Mo_mat = km_mat;
tr_mat = km_mat;
ts_mat = km_mat;
feas_mat = km_mat;
t = [0:0.01:30]';
for ii = 1:len_KI
 for jj = 1:len_KP
  K_V = K_P(jj)*tf([(1+K_I(ii)*Ts) -1],[1 -1],Ts);
  [km,pm,gw,pw] = margin(G1*K_V);
  G1cl = feedback(G1*K_V,1);
  y = step(G1cl,t);
  y = 0.05*y;
  [Mo,tp,tr,ts,ess] = kstats(t,y,0.05);
  if isempty(ts) == 1
      ts = max(t);
  end
  km_mat(jj,ii) = km;
  pm_mat(jj,ii) = pm;
  Mo_mat(jj,ii) = Mo;
  tr_mat(jj,ii) = tr;
  ts_mat(jj,ii) = ts;
      if tr <= 1
          if ts <= 5
              feas_mat(jj,ii) = 1;
          end
      end
  end
end
disp('Values of KP as rows of feasibility matrix')
disp(K_P)
disp('Values of KI as columns of feasibility matrix')
disp(K_I)
disp('Feasibility matrix')
disp(feas_mat)

% selection of gain
disp('A proportional gain of K_P = 25 and K_I = 0.1 s will result in')
disp('a rise time of less than 0.8 s,') 
disp('and a settling time of less than 5 s.')

disp('*****>'), pause
K_V = 25*tf([(1+0.1*Ts) -1],[1 -1],Ts)
[km,pm,gw,pw] = margin(G1*K_V)
G1cl = feedback(G1*K_V,1);
[y,t] = step(G1cl);
y = 0.05*y;
figure
plot(t,y,'o'), 
xlabel('time (s)')
ylabel('voltage amplitude')
title('Step response for Comprehensive Problem 8.3')
hold off, grid

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