cp7_3.m

离散控制系统设计的MATLAB 代码

%%%%%%%%%%% Comprehensive Problem 7.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                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hydro-turbine System

clear
disp('Comprehensive Problem 7.3')
disp('PI 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 PI controller')
disp('                     Ts * z  ')
disp('Gc(z) = K_P (1 + K_I ------ )')
disp('                       z-1   ')

disp('Varying K_P and K_I')
disp('Integral gain K_I')
K_I = [0.5:0.1:1.2]
len_KI = length(K_I);
disp('*****>'), pause

disp('Proportional gain') 
K_P = [0.1:0.05:0.4]
len_KP = length(K_P);
disp('*****>'), pause

disp('Sweeping of parameters ...')
km_mat = zeros(len_KP,len_KI);
pm_mat = km_mat;
Mo_mat = km_mat;
tr_mat = km_mat;
ts_mat = km_mat;
t = [0:Ts:50]';
for ii = 1:len_KI
 for jj = 1:len_KP
  G_c = K_P(jj)*(1 + tf(K_I(ii)*Ts*[1 0],[1 -1],Ts));
  [km,pm,gw,pw] = margin(Gz*G_c);
  G_cl = feedback(Gz*G_c,1);
  y = step(G_cl,t);
  [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,ii) = km;
  pm_mat(jj,ii) = pm;
  Mo_mat(jj,ii) = Mo;
  tr_mat(jj,ii) = tr;
  ts_mat(jj,ii) = ts;
 end
end

figure, mesh(K_I,K_P,km_mat)
xlabel('K_I'), ylabel('K_P')
zlabel('Gain margin (dB)') 
title('Gain margin k_m vs K_I and K_P')
disp('Plotting gain margin k_m vs K_I and K_P')
disp('*****>'), pause

figure, mesh(K_I,K_P,pm_mat)
xlabel('K_I'), ylabel('K_P')
zlabel('Phase margin (deg)') 
title('Phase margin p_m vs K_I and K_P')
disp('Plotting phase margin p_m vs K_I and K_P')
disp('*****>'), pause

figure, mesh(K_I,K_P,Mo_mat)
xlabel('K_I'), ylabel('K_P')
zlabel('Overshoot (%)') 
title('Overshoot M_o vs K_I and K_P')
disp('Plotting overshoot M_o vs K_I and K_P')
disp('*****>'), pause


figure, mesh(K_I,K_P,tr_mat)
xlabel('K_I'), ylabel('K_P')
zlabel('Rise time (s)') 
title('Rise time t_r vs k_I and K_P')
disp('Plotting rise time t_r vs K_I and K_P')
disp('*****>'), pause

figure, mesh(K_I,K_P,ts_mat)
xlabel('K_I'), ylabel('K_P')
zlabel('Settling time (s)') 
title('Settling time t_s vs K_I and K_P')
disp('Plotting settling time t_s vs K_I and K_P')
disp('*****>'), pause

% selection of gain
figure
[cpm,hpm] = contour(K_I,K_P,pm_mat,'--'); % plot level contour lines
clabel(cpm,hpm);                          % label contour lines
grid
xlabel('K_I')
ylabel('K_P')
hold on
[cts,hts] = contour(K_I,K_P,ts_mat);     % plot level contour lines
clabel(cts,hts);                         % label contour lines
title('t_s in s (solid), pm in deg (dashed)')
hold off
disp('Contour plot of settling time and phase margin vs K_I and K_P')
disp('*****>'), pause

disp('From the contour plot, minimizing t_s while maintaining pm > 60 deg,')
disp('we need to select K_P = 0.2 and K_I = 1')
disp(' ')

jj = 3; ii =6;
G_c = K_P(jj)*(1 + tf(K_I(ii)*Ts*[1 0],[1 -1],Ts));
figure, margin(Gz*G_c); 
disp('Bode plot to show gain and phase margins')
disp('*****>'), pause
G_cl = feedback(Gz*G_c,1);
figure
[y,dt] = step(G_cl);
plot(dt,y,'o'),grid
xlabel('Time (s)')
ylabel('Amplitude')
title('CP7.3: Closed-loop Response of Power P to Unit Step of Reference Input')

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