ex8_2a.m

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

%%%%%%%%%%%%%%%%% Example 8.2(a) %%%%%%%%%%%%%%%%
%   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                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   ----  Proportional-integral control  ----
%
clear
disp('Example 8.2(a)')

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

% --- compute CL step responses for range of KP values
time  = 0:Ts:8;
KI = input('Enter integral gain KI ==> [0.5] ');  % user specifies KI, default is  0.5
if isempty(KI), KI = 0.5; end
figure, hold on
symbols = ['o' '+' '*' 'd' 'v'];         % different plot symbols
results = [];
Pgains = [];
for ii = 1:5
    KP = 0.3+0.2*(ii-1);                  % KP = 0.3, 0.5, ..., 1.1
    Gc  = KP*tf([1+KI*Ts  -1],[1 -1],Ts); % PI control law
    T = Gp_z*Gc/(1+HGp_z*Gc);             % closed-loop system
    T = minreal(T);                       % remove pole-zero cancellations
    y = step(T,time);                     % closed-loop step response
    ref = dcgain(T);
    [Mo,tp,tr,ts,ess] = kstats(time,y,ref);  % response stats
    ymax = y(find(time == tp));
    yss = y(length(time));
    results = [results; [KP  ts  ymax  tr  tp  yss  Mo]];
    Pgains = [Pgains KP];
    disp(['KP = ',num2str(KP)])
    % plot every other point
    plot(time(1:2:length(time)),y(1:2:length(time)),symbols(ii))
end
xlabel('Time Sample (s)')
ylabel('Amplitude')
title(['Example 8.2(a) - Step Responses for Different ',...
       'Values for KP, where KI = ',num2str(KI)]) 
grid
legend('KP = 0.3','KP = 0.5','KP = 0.7','KP = 0.9','KP = 1.1',4)
hold off                                
%------ show table of KP values and plot symbols -----
while 0
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.3    o')
text('Units','norm','Pos',[0.5,0.20],'Str','0.5    +')
text('Units','norm','Pos',[0.5,0.15],'Str','0.7    *')
text('Units','norm','Pos',[0.5,0.10],'Str','0.9  diamond')
text('Units','norm','Pos',[0.5,0.05],'Str','1.1  triangle')
end

text('Units','norm','Pos',[0.75,0.45],'Str',['KI = ',num2str(KI)])

%----- print results table
disp(' ')
disp('      KP      tsettle   y_max     t_rise   t_peak     y_ss    M_o (pct)')
disp('    -------------------------------------------------------------------')
disp(results)
disp('    -------------------------------------------------------------------')
%%%%%%%%%%