⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 cp9_2.m

📁 离散控制系统设计的MATLAB 代码
💻 M
字号:
%%%%%%%%%%% Comprehensive Problem 9.2 %%%%%%%%%%%
%   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                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Electric Power Generation System

clear all, close all
disp('Comprehensive Problem CP9.2')

dpower2    % read data
disp('*****>'), pause

disp('Power system model state matrices, with Vterminal, ')
disp('  w(speed), and P(power) as output')
disp('Continuous-time model')
Gs
disp('*****>'), pause

disp('discretized model with a sampling period of 0.01 s')
Ts = 0.01
Gz = c2d(Gs,Ts) 
disp('Note that most of the zero entries of the A & B matrices')
disp('  of the continuous time model have become nonzero in the')
disp('  discrete time models. The C & D matrcies are unchanged.')
disp('*****>'), pause

disp('System model with Vterminal as the only output')
Gv = Gz(1,1)
disp('*****>'), pause

Gvw = d2c(Gv,'tustin');           % bilinear transform 
%-------------- open-loop frequency response
ww = logspace(-2,2,100)';         % use 100 points for good resolution
[mag,ph] = bode(Gvw,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 90 deg phase margin
%--- allows for 5 extra degrees in the design
ph_des = -180 + (90+5);
mag_des = interp1(ph,mag,ph_des) 
Klag = 1/mag_des                  % set gain for selected phase 
lfg = Klag*dcgain(Gvw)            % lfg for plant + gain
Alag = 49/lfg                     % lag alpha supplies rest of req'd low-freq gain
                                  % 49 translate into dcgain of 50 and 
                                  % steady-state error of 2%
wwc = interp1(ph,ww,ph_des)

figure
bode(Gvw,ww),grid
title('Open-loop frequency response ')   
disp('Plotting open-loop frequecy 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
GKvw = Gvw*Klagw;                 % connect lag in series with plant and sensor
lfg = dcgain(GKvw)                % low-freq gain of current design
[km,pm,wkm,wpm] = margin(GKvw);   % open-loop frequency response
disp('gain margin     phase crossover frequency')
disp([20*log10(km)  wkm])         % gain margin in dB & omega_gm
disp('phase margin    gain crossover frequency')
disp([pm wpm])                    % phase margin in degrees & omega_pm
ess = 1/(1 + dcgain(GKvw))        % verify that steady-state error is OK
Klagz = c2d(Klagw,Ts,'tustin');   % controller in z domain

%-----------  OL frequency responses ---------
w = ww;
ph180 = -180*ones(length(w));
%---- first, with gain only
[mag_K,ph_K] = bode(Gvw*Klag,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(GKvw,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;
legend('compensated','uncompensated')
xlabel('Frequency (rad/s)')
ylabel('Magnitude (dB)')
title('Compensated frequency responses')
subplot(212)
semilogx(w,ph_lag,'-',w,ph_K,'--',w,ph180,':');grid
legend('compensated','uncompensated')
xlabel('Frequency (rad/s)')
ylabel('Phase angle (deg)')
disp('Plotting compensated frequecy response')
disp('*****>'), pause

%-----------  CL step responses  -------
Tz = feedback(Gv*Klagz,1);      % closed-loop transfer function
[y,tt] = step(Tz,10);       % step response to 10 s
figure
plot(tt,0.05*y,'o'), hold on
yss = dcgain(Tz);
plot([0 10],[0.05*yss 0.05*yss],'--'), hold off
xlabel('Time (s)')
ylabel('Amplitude')
title('Step response')
grid
disp('Plotting closed-loop system step response')
disp('end of Comprehensive Problem 9.2')
%%%%%%%%%%

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -