% ***** MATLAB Code Starts Here *****
%
olzp = [-1; -1]; % plant poles and zeros
olpp = [0; -0.01; -0.02+j*0.9998; -0.02-j*0.9998; -10];
%
nump = 0.01 * real(poly(olzp)); % plant numerator
denp = real(poly(olpp)); % and denominator
%
w = logspace(-4,2,601); % frequency vector
%
[magp,php] = bode(nump,denp,w); % plant frequency response
figure(1),clf,semilogx(w,20*log10(magp),w,php),grid,...
xlabel('Frequency (r/s)'),ylabel('Magnitude (db) & Phase (deg)'),...
title('Bode Plots of Uncompensated System'),...
hold on,semilogx([min(w) max(w)],[-180 -180],'r--'),hold off
%
ess_spec = 0.1; % specifications
pm_spec = 50;
wx_spec = 0.05;
%
Kvp = 0.01/(0.01*10); % plant velocity error constant
Kv_spec = 1/ess_spec;
%
Kc1 = Kv_spec / Kvp; % calculate compensator gain
%
% Original and gain compensated Bode plots
%
figure(2),clf,semilogx(w,20*log10(Kc1*magp),w,20*log10(magp),'k--',w,php),grid,...
xlabel('Frequency (r/s)'),ylabel('Magnitude (db) & Phase (deg)'),...
title('Bode Plots of Uncompensated System with K_c = 100'),...
hold on,semilogx([min(w) max(w)],[-180 -180],'r--'),hold off
%
% Find uncompensated phase shift at desired gain crossover frequency
%
[i0,i1] = min(abs(w-wx_spec));
wxc1 = w(i1);
phpwxc = php(i1);
%
% Design a single stage phase lead compensator to provide the
% necessary positive phase shift at the desired gain crossover frequency
%
phimax = (pm_spec + 10) - (180 + phpwxc);
alfad = (1 - sin(phimax*pi/180)) / (1 + sin(phimax*pi/180));
%
zcd = wxc1 * sqrt(alfad); % lead compensator zero & pole
pcd = zcd / alfad;
%
numcd = [1/zcd 1]; % lead compensator numerator
dencd = [1/pcd 1]; % and denominator
%
numf1 = Kc1 * conv(numcd,nump);
denf1 = conv(dencd,denp);
%
[magf1,phf1] = bode(numf1,denf1,w); % frequency response
%
figure(3),clf,semilogx(w,20*log10(magf1),w,20*log10(Kc1*magp),'k--',w,phf1,...
w,php,'k--'),grid,...
xlabel('Frequency (r/s)'),ylabel('Magnitude (db) & Phase (deg)'),...
title('Bode Plots of Lead Compensated System'),...
hold on,semilogx([min(w) max(w)],[-180 -180],'r--'),hold off
%
% Design a two stage phase lag compensator to drop the magnitude
% curve of the lead-compensated system down to 0 db at the
% desired gain crossover frequency
%
alfag = sqrt(magf1(i1));
%
zcg = wxc1 / 20; % lag compensator zero & pole
pcg = zcg / alfag;
%
numcg = conv([1/zcg 1],[1/zcg 1]); % lag compensator numerator
dencg = conv([1/pcg 1],[1/pcg 1]); % and denominator
%
numf = conv(numf1,numcg);
denf = conv(denf1,dencg);
%
[magf,phf] = bode(numf,denf,w); % frequency response
%
figure(4),clf,semilogx(w,20*log10(magf),w,20*log10(magf1),'k--',w,phf,...
w,phf1,'k--'),grid,...
xlabel('Frequency (r/s)'),ylabel('Magnitude (db) & Phase (deg)'),...
title('Bode Plots of Final Lag-Lead Compensated System'),...
hold on,semilogx([min(w) max(w)],[-180 -180],'r--'),hold off
%
[i0,i1] = min(abs(magf - 1));
wxc = w(i1); % actual compensated gain crossover
pmc = 180 + phf(i1); % frequency & phase margin
%
t = linspace(0,1000,2001); % time vector
%
dclp = denp + [0 0 0 nump]; % closed-loop denominators
dclf = denf + [0 0 0 numf];
%
ysp = step(nump,dclp,t); % step responses
ysf = step(numf,dclf,t);
%
figure(5),clf,plot(t,ysp,t,ysf),grid,xlabel('Time (s)'),...
ylabel('Step Responses'),title('Uncompensated and Final Compensated Systems'),...
text(140,1.45,'Uncompensated')
%
yrp = step(nump,[dclp 0],t); % ramp responses
yrf = step(numf,[dclf 0],t);
yrp1 = step(nump,[dclp 0],3*t); % ramp responses
yrf1 = step(numf,[dclf 0],3*t);
%
figure(6),clf,subplot(121),plot(t,yrf,t,yrp,[0 1000],[0 1000],'r--'),grid,...
xlabel('Time (s)'),ylabel('Ramp Responses'),...
title('Uncompensated and Final Compensated Systems'),...
axis([0 1000 0 1000]),...
%
subplot(122),plot(3*t,yrf1,3*t,yrp1,[0 3000],[0 3000],'--'),grid,...
xlabel('Time (s)'),ylabel('Ramp Response'),...
title('Expanded View'),...
axis([2980 3000 2980 3000])
%
% ***** MATLAB Code Stops Here *****
Click the
icon to return to the Dr. Beale's home page
Lastest revision on
Friday, May 26, 2006 9:07 PM