ex8_1b.m

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

%%%%%%%%%%%%%%%%% Example 8.1(b) %%%%%%%%%%%%%%%%
%   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                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  ---- Step responses with proportional gain sweep  ----
%
clear
disp('Example 8.1(b)')

Ts = 0.1;			               % sampling period
Gp_s = tf(4,conv([2 1],[0.5 1]));  % continuous process
H_s = tf(1,[0.05 1]);		       % continuous sensor
Gp_z = c2d(Gp_s,Ts,'zoh');	       % discrete plant (no sensor)
HGp_s = H_s*Gp_s                   % continous plant
HGp_z = c2d(HGp_s,Ts,'zoh')        % discrete plant, using ZOH option

results = [];                      % array to save results for building table
Pgains = [];                       % array for saving gains
dtime = 0:Ts:6;                    % vector of time points separated by Ts
figure,hold on               	   % turn hold on before 1st plot
symbols = ['o' '+' '*' 'd' 'v'];   % different plot symbols
for ii = 1:5
    KP = 0.5*ii;                   % KP = 0.5, 1.0, ..., 2.5
    T = Gp_z*KP/(1+HGp_z*KP)
    T = minreal(T);                % remove pole-zero cancellations
    y = step(T,dtime); 
    ref = dcgain(T);               % steady-state value of y
    [Mo,tpeak,trise,tsettle,ess] = kstats(dtime,y,ref);  %calc performance
    ymax = y(find(dtime == tpeak));
    yss = y(length(dtime));
    %-- save gains and result for making table at end ----
    results = [results; [KP  ymax  trise  tpeak  yss  Mo]];
    Pgains = [Pgains KP];
    disp(['KP = ',num2str(KP)])
    % plot every other point
    plot(dtime(1:2:length(dtime)),y(1:2:length(dtime)),symbols(ii))
    pause(1)                        % wait 1 second before doing next gain
end
hold off                            % turn hold off after last curve
xlabel('Time (s)')
ylabel('Amplitude')
title(['Example 8.1(b) - Step Responses for Different ',...
       'Proportional Gains for HGp\_z']) 
grid

text(1.3,1.30,'KP = 2.5')
text(1.5,0.55,'KP = 0.5')
legend('KP = 0.5','KP = 1.0','KP = 1.5','KP = 2.0','KP = 2.5')

while 0
%------ show table of KP values and plot symbols -----
text('Units','norm','Pos',[0.5,0.35],'Str','KP symbol')
text('Units','norm','Pos',[0.5,0.30],'Str','---------')
text('Units','norm','Pos',[0.5,0.25],'Str','0.5    o')
text('Units','norm','Pos',[0.5,0.20],'Str','1.0    +')
text('Units','norm','Pos',[0.5,0.15],'Str','1.5    *')
text('Units','norm','Pos',[0.5,0.10],'Str','2.0  diamond')
text('Units','norm','Pos',[0.5,0.05],'Str','2.5  triangle')
end

%----- print results table
disp(' ')
disp('      KP      y_max     t_rise   t_peak     y_ss    M_o (pct)')
disp('    -----------------------------------------------------------')
disp(results)
disp('    -----------------------------------------------------------')
disp('Hit any key to continue..'); pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by trial-and-error, find that 20 pct overshoot occurs 
%	for KP = 1.032 --> 19.9808 pct
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Gain and phase margins
%----- return values but no plot
KP = 1.032
[Gm,Pm,Wcg,Wcp] = margin(KP*HGp_z)
%----- do plot only
figure
margin(KP*HGp_z)
%%%%%%%%%%