ex9_3.m

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

%%%%%%%%%%%%%%%%%% Example 9.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                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   ---- Lead controller design ----
%
clear
disp('Example 9.3')
%---- build model with cmds from Ex. 9.2(a)
tauP1 = 2; tauP2 = 0.5; tau_sen = 0.05; Kproc = 4; % parameters
Gp = tf(Kproc,conv([tauP1 1],[tauP2 1]));          % plant Gp(s)
H  = tf(1,[tau_sen 1]);                            % sensor H(s)
%--- OL transfer function is plant * sensor in s-domain
GH = H*Gp;           
Ts = 0.1;                                      % sampling period
GHz = c2d(GH,Ts,'zoh')         
GHw = d2c(GHz,'tustin')                        % bilinear transform 
disp('******>'), pause 

%----- Part (a) determine lead gain -----
Kp_reqd = (1/0.1) - 1                          % required gain to give 4% ess
plant_lfg = dcgain(GHw)
Klead = Kp_reqd/plant_lfg

%----- Part (b) determine lead alpha -----
[gm,pm,wwgm,wwpm] = margin(Klead*GHw);
pm_target = 65                                 % phase-margin target in degrees
del_pm = pm_target - pm                        % max phase angle in degrees 
Alead = (1+sin(del_pm*pi/180))/(1-sin(del_pm*pi/180)) % Eq (9.21)
disp('******>'), pause 

%----- Part (c) determine lead zero -----
ww = logspace(-1,2,100);                       % set up frequency scale
[mag_db, ph] = bodedb(Klead*GHw,ww);           % compute frequency response
mag_ratio = 10.^(mag_db/20);                   % convert to ratio
wwc = interp1(mag_ratio,ww,1/Alead)            % find wwc
Zlead = -wwc/sqrt(Alead)                       % compute zero
disp('******>'), pause 

Gc1 = tf(Alead*[1 -Zlead],[1 -Alead*Zlead])    % lead section
Gcw = Klead*Gc1                                % lead controller
GcGHw = Gcw*GHw;                               % OL tf of lead, plant, & sensor
[kmc,pmc,wwkm,wwpm] = margin(GcGHw)            % find new phase margin
disp('******>'), pause 

figure
margin(GcGHw)
disp('******>'), pause 

%--------------  OL frequency responses for Fig. 9.7 ---------
w = logspace(-2,3);
ph180 = -180*ones(length(w));
%---- first, with gain only ----------------
[mag_db_K,ph_K] = bodedb(Klead*GHw,w);         % returns vectors, w/ mag in dB
%-------- then with final lag design --------
[mag_db_lead,ph_lead] = bodedb(GcGHw,w);       % returns vectors, w/ mag in dB
%--------- now plot both sets ----------------
figure
subplot(211)
semilogx(w,mag_db_lead,'-',w,mag_db_K,'--');grid;
xlabel('Frequency (rad/s)')
ylabel('Magnitude (dB)')
title('Bode plots of initial (proportional) and final (lead) designs for Example 9.3')
text(2,-38,'Proportional')
text(14,-10,'lead')
subplot(212)
semilogx(w,ph_lead,'-',w,ph_K,'--',w,ph180,':');grid
text(0.7,-220,'Proportional')
text(20,-200,'lead')
xlabel('Frequency (rad/s)')
ylabel('Phase angle (deg)')
disp('******>'), pause 

%----- Part (d) CL system for ref input with sensor in fdbk path -----
Gcz = c2d(Gcw,Ts,'tustin')            % convert to DTS controller
Gz = c2d(Gp,Ts,'zoh')
Tz = Gcz*Gz/(1+Gcz*GHz);              % closed-loop transfer function
[y,tt] = step(Tz,3);                  % step response to 3 s

figure
plot(tt,y,'*');grid
hold on
yss = dcgain(Tz);
plot([0 3],[yss yss],'--')
hold off
xlabel('Time (s)')
ylabel('Amplitude')
title('Unit-step response for the final lead design of Example 9.3')
%%%%%%%%%%