ep9_1.m

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

%%%%%%%%%%%% Exploratory problem 9.1 %%%%%%%%%%%%
%   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                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   ---- Lag-compensator design for type-0 plant ----
%
clear
disp('Exploratory Problem 9.1')

Gp = tf(160*[1 10],conv(conv([1 2],[1 4]),conv([1 5],[1 20])))
Ts = 0.02                         % sampling period
Gz = c2d(Gp,Ts,'zoh')         
[zGz,pGz,kGz] = zpkdata(Gz,'v')   % zeros, poles, & gain of GH(z)
Gw = d2c(Gz,'tustin')             % bilinear transform 
[zGw,pGw,kGw] = zpkdata(Gw,'v')   % zeros, poles, & gain of GH(w)

%-------------- open-loop frequency response
ww = logspace(-1,1,100)';         % use 100 points for better resolution
[mag_db,ph] = bodedb(Gw,ww);      % magnitude & phase as 1:1:100 arrays
mag_ratio = 10.^(mag_db/20);      % convert to ratio
% desired phase margin is 60 deg
%-------- display table of index, phase, mag, & freq values
%           by selecting -120 <= phase <= -110 deg
for ii = find((ph <= -110)& (ph >= -120)),   
  disp([ ii  ph(ii)  mag_ratio(ii)  ww(ii)]) 
end 
%--- interpolate to get magnitude for 60+5 deg phase margin
mag115 = interp1(ph,mag_ratio,-115) 
Klag = 1/mag115                  % set gain for selected phase 
lfg = Klag*dcgain(Gw)            % lfg for plant + gain
Alag = 49/lfg   % lag alpha supplies rest of req'd low-freq gain
wwc = interp1(ph,ww,-115)
figure, bode(Gw,ww)

w_Zlag = input('Enter corner freq for lag zero 0.1*2.676 ..... ') 
Zlag = -w_Zlag                   % lag zero must be negative
Klagw = tf(Klag*[1 -Zlag],[1 -Zlag/Alag])   % lag controller
KGw = Klagw*Gw     % connect lag in series with plant and sensor
lfg = dcgain(KGw)                % low-freq gain of current design
[km,pm,wkm,wpm] = margin(KGw);   % open-loop frequency response
disp('Gain margin in dB and omega_gm')
disp([20*log10(km)  wkm])        % gain margin in dB & omega_gm
disp('Phase margin in degrees and omega_pm')
disp([pm wpm])                   % phase margin in degrees & omega_pm
ess = 1/(1 + dcgain(KGw))        % verify that steady-state error is OK
Klagz = c2d(Klagw,Ts,'tustin')   % controller in z domain

w = logspace(-3,1);
ph180 = -180*ones(length(w));
%---- first, with gain only -----
[mag_db_K,ph_K] = bodedb(Klag*Gw,w); % bodedb returns mag in dB
%-------- then with final lag design
[mag_db_lag,ph_lag] = bodedb(KGw,w);
%--------- now plot both sets
figure
subplot(211)
semilogx(w,mag_db_lag,'-',w,mag_db_K,'--');grid;
text(0.004,7,'Proportional')
text(0.05,30,'lag')
xlabel('Frequency (rad/s)')
ylabel('Magnitude (dB)')
subplot(212)
semilogx(w,ph_lag,'-',w,ph_K,'--',w,ph180,':');grid
text(0.3,-25,'Proportional')
text(0.02,-70,'lag')
xlabel('Frequency (rad/s)')
ylabel('Phase angle (deg)')

% reference response
Gz = c2d(Gp,Ts,'zoh')         % discretize plant
Klagz = ss(Klagz)             % convert to ss form for better numerics
Tz = Klagz*Gz/(1+Klagz*Gz)    % closed-loop transfer function
[y,tt] = step(Tz,8);          % step response to 8 s

yss = dcgain(Tz);

figure
plot(tt,y,'o');grid
hold on
plot([0 8],[yss yss],'--')
hold off
xlabel('Time (s)')
ylabel('Amplitude')
%%%%%%%%%%