cp9_3.m

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

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

clear
disp('Comprehensive Problem 9.3')
disp('Lag control of hydro-turbine system')

dhydro                      % read data
disp('discrete-time transfer function')
Gz = tf(dnum,dden,Ts)       % form tf object
disp('*****>'), pause

Gw = d2c(Gz,'tustin')       % bilinear transform 
%-------------- open-loop frequency response
ww = logspace(-3,2,100)';   % use 100 points for good resolution
disp('*****>'), pause

disp('Apply pre-filter')
Gw_pc = tf(0.05,[1 0.05]);
disp('*****>'), pause

Gw = Gw*Gw_pc;
[mag,ph] = bode(Gw,ww);     % computes magnitude & phase as 1x1x100 arrays
mag = reshape(mag,1,100);   % reshape arrays to be vectors
ph = reshape(ph,1,100);

%--- interpolate to get magnitude for 45 deg phase margin
%--- allows for 5 extra degrees in the design
ph_des = -180 + (45+5);
mag_des = interp1(ph,mag,ph_des) 
Klag = 1/mag_des            % set gain for selected phase 
lfg = Klag*dcgain(Gw)       % lfg for plant + gain
Alag = 19/lfg               % lag alpha supplies rest of req'd low-freq gain
                            % 19 translate into dcgain of 20 and 
                            % steady-state error of 5%
wwc = interp1(ph,ww,ph_des)

figure
bode(Gw,ww),grid
title('Open-loop frequency response ')   
disp('Plotting open-loop frequency response')
disp('*****>'), pause

% Select lag zero & calc phase margin
Zlag = -wwc/10              % set lag zero  
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(dB)   phase crossover frequency(rad/s)')
disp([20*log10(km)  wkm])   % gain margin in dB & omega_gm
disp('phase margin(deg)   gain crossover frequency(rad/s)')
disp([pm wpm])              % phase margin in degrees & omega_pm
ess = 1/(1 + dcgain(KGw))   % verify that steady-state error is OK
Klagz = c2d(Klagw*Gw_pc,Ts,'tustin')      % controller in z domain
disp('*****>'), pause

%========== frequency responses =====
w = ww;
ph180 = -180*ones(length(w));

[mag_K,ph_K] = bode(Klag*Gw,w);
mag_K = reshape(mag_K,1,100);
ph_K  = reshape(ph_K,1,100);
mag_db_K = 20*log10(mag_K);
%-------- then with final lag design
[mag_lag,ph_lag] = bode(KGw,w);
mag_lag = reshape(mag_lag,1,100);
ph_lag  = reshape(ph_lag,1,100);
mag_db_lag = 20*log10(mag_lag);
%--------- now plot both sets
figure
subplot(211)
semilogx(w,mag_db_lag,'-',w,mag_db_K,'--');grid;
xlabel('Frequency (rad/s)')
ylabel('Magnitude (dB)')
title('Compensated frequency response')
subplot(212)
semilogx(w,ph_lag,'-',w,ph_K,'--',w,ph180,':');grid
xlabel('Frequency (rad/s)')
ylabel('Phase angle (deg)')
disp('Plotting compensated frequency response')
disp('*****>'), pause

%==========  CL step responses  =====
Tz = feedback(Klagz*Gz,1)     % closed-loop transfer function  
[y,tt] = step(Tz,50);         % step response to 20 s
figure
plot(tt,y,'o'), hold on
yss = dcgain(Tz);
plot([0 50],[yss yss],'--'), hold off
xlabel('Time (s)')
ylabel('Amplitude')
title('Step response')
grid

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