Example H8.3

MATLAB code for example 8.3 from the book "Regeltechniek voor het HBO"

Controller gain for an given overshoot of a third order system with unit feedback

Copyright (c) 2021, Studieboeken Specialist Permission is granted to copy, modify and redistribute this file, provided that this header message is retained.
Table of Contents
Controller gain for an given overshoot of a third order system with unit feedback Assignment Solution MATLAB code for this example The given situation and assignment Define model parameters Create proces and plot stepresponse Plot root locus Create controller and analyze closed loop Generate MATLAB figure(s) for usage in the book Init create Enhanced Figures Enhance the figures Function pool

Assignment

In figure 8.8 a third order process is controlled by a P-controller with gain .
  1. What should be the K-factor of this controller if the overshoot in a step response is allowed to be about 20%?
  2. For this also determine the settling time .

Solution

The root locus for is shown in figure 8.9 (see figure 6.11-l, or simulate the root locus in MATLAB as shown in example 6.12). The root locus equation is:
.
From the overshoot can be determined:
→ (note: here you can see what we do very often: using a formula that applies to a pure second order proces for a higher order process that is "dominant second order" in the points of figure 8.9, because the third pole is already (far?) away to the left on the real axis).
From we get , but we want to continue with the , λ, ω parameters.
Because , we can describe in figure 8.9 as where . Substituted in the root locus equation we get:
→ .
These are two equations in one:
and
.
From the latter we get the (mathematical) solutions: , and .
But we can't accept physically as we can see directly in figure 8.9. The same holds for . That would mean in figure 8.9 that ! So is the only one we can acccept:
and .
With the sum rule (see formula 6.35 and ) we easily find the third pole: → .
Substituting this value in the root locus equation gives us in the points of figure 8.9:
→ , and that means an area of proportionality of the controller of (in practice the area of proportionality is an important aspect of the controller).
For and thus we can determine the settling time in a step response:
→ → seconds.
At a step response we also are interested in the static behaviour of the controlled process:
for .
So the accuracy is and the static error .
A larger give a smaller error for the static state, but a larger overshoot and a higher frequency in the transient. And the risk of oversteering the actuator is present. Like in real live, there is also give and take in control engineering!
In this example we had to do a lot of calculations. In fact, we need to let the computer do that as much as possible so that we can focus on understanding. That's all that really matters.

MATLAB code for this example

% clear all variables from Workspace and close all figures.
clear variables;
close all;
 
% Define 's' variable
s=tf('s');

The given situation and assignment

Assume that a third order system with unit feedback is controlled by a P-controller (figure 8.8 in the book). The third order system is given as:
Find the gain of the proportional controller, such that the steprespons has a overshoot of 20%. Determine also the setlling time for this gain.

Define model parameters

p1=0.5;
Kpn=0.1;

Create proces and plot stepresponse

Hp=Kpn / ( (s+p1)^3 );
figure(101);
step(Hp);
grid on;
title('Stepressponse openloop process');

Plot root locus

Determine gain of the P-controller
figure(102);
h=rlocusplot(Hp);
grid on;
% get the handle to the plot options
p=getoptions(h);
% Zoom in
p.XLim=[-1 0.5];
p.YLim=[-1 1];
% Adjust Title
p.Title.String='Root Locus of H_{p}';
% write the options to the figure
setoptions(h,p);
Read the gain from the plot where the rootlocus intersects with the line from the origin of the plot under 45 degrees. That is where Beta = 0.45 and the overshoot is max 20%. To do this in MATLAB you might have to go the command Window and paste the code above again to get a picture/plot outside the LiveScript environment where you can right click on it to read the rootlocus gain values.
Kreg=1.57;

Create controller and analyze closed loop

REG=Kreg;
 
Hcl=feedback(REG*Hp,1);
figure(2);
step(Hcl,50);
grid on;
title('Stepresponse closed loop');
displayCharacteristics(Hcl,0.02)

Generate MATLAB figure(s) for usage in the book

Init create Enhanced Figures

Close all the earlier enhanced figures with a certain tag
EnhancedFig = findobj(0, 'Tag', 'EnhancedImage');
close(EnhancedFig);

Enhance the figures

figure1=figure(1001);
% Give it a tag
set(gcf, 'Tag','EnhancedImage');
% create the plot again but now with rlocusplot to return a figure handle
h=rlocusplot(Hp,'k');
% Determine Kpn
[Z,P,Kpn]=zpkdata(Hp);
% get the handle to the plot options
p=getoptions(h);
% Zoom in
Scale=0.6;
p.XLim=[-Scale Scale];
p.YLim=[-Scale Scale];
% Adjust Title
p.Title.String='Root locus plot of H_p';
% write the options to the figure
setoptions(h,p);
 
%%Add text and labels
% Calculate location axis
Xoffset=0.14; % left part of the figure window not used by the graph
Scale=0.8; % part of the whole figure window used by the graph
Vaxis_loc=Xoffset+0.4;
% Place axis labels
omega = char(969); % Greek letter omega
lambda = char(955); % Greek letter lambda
annotation(figure1,'textarrow',[Vaxis_loc Vaxis_loc],...
[0.79 0.86],'String',{omega});
annotation(figure1,'textarrow',[0.9-0.07 0.9],...
[0.48 0.48],'String',{lambda});
annotation(figure1,'textbox',...
[Vaxis_loc 0.52 0.077 0.051],...
'String',{'s-vlak'},...
'LineStyle','none',...
'FitBoxToText','off');
annotation(figure1,'textbox',...
[Vaxis_loc 0.46 0.0306 0.042],...
'String','0',...
'LineStyle','none',...
'FitBoxToText','off');
annotation(figure1,'textbox',...
[Vaxis_loc+0.01 0.81 0.077 0.0517],...
'String',['Kpn=', num2str(Kpn)],...
'LineStyle','none',...
'FitBoxToText','off');

Function pool

=============
The code below is only meant to supply functions that can be used in the script above
=============
function displaySettlingTime(sys, threshold)
[y,t] = step(sys);
stepResults = stepinfo(y,t,'SettlingTimeThreshold',threshold);
indexSettling = find(t>=stepResults.SettlingTime,1);
hold on;
plot(t(indexSettling),y(indexSettling),'ro');
plot([0 t(indexSettling)],[y(indexSettling) y(indexSettling)],'g-.');
plot([t(indexSettling) t(indexSettling)],[0 y(indexSettling)],'g-.');
hold off;
end
function displayPoles(sys,zeta,wn)
sigma = zeta*wn;
wd = wn*sqrt(1-zeta^2);
pzmap(sys)
hold on;
plot([0 -sigma], [0 wd],'r:');
text(-sigma/2+0.3,wd/2+0.5,['\omega_n']);
text(-sigma/3,0.5,['\beta = cos^{-1}(\zeta)']);
plot([-sigma -sigma], [0 wd],'r:');
text(-sigma-0.3,-0.5,['-\sigma']);
plot([0 -sigma], [wd wd],'k:');
plot([0 -sigma], [0 0],'r:');
text(0.2,wd,['j\omega_d']);
hold off;
end
function displayCharacteristics(sys,threshold)
[y,t] = step(sys);
stepResults = stepinfo(y,t,'SettlingTimeThreshold',threshold);
indexSettling = find(t>=stepResults.SettlingTime,1);
hold on;
plot(t(indexSettling),y(indexSettling),'ro');
plot([0 t(indexSettling)],[y(indexSettling) y(indexSettling)],'k-.');
plot([t(indexSettling) t(indexSettling)],[0 y(indexSettling)],'k-.');
text(stepResults.SettlingTime+0.1,0.9,['t_s \approx ',num2str(stepResults.SettlingTime),' sec']);
plot(stepResults.PeakTime,stepResults.Peak,'ro');
plot([stepResults.PeakTime stepResults.PeakTime], [0 stepResults.Peak],'k-.');
plot([0 stepResults.PeakTime], [stepResults.Peak stepResults.Peak],'k-.');
text(stepResults.PeakTime+0.1,0.1,['t_p = ',num2str(stepResults.PeakTime),' sec']);
text(stepResults.PeakTime+0.1,stepResults.Peak+0.03,['M_p = ',num2str(stepResults.Overshoot),'%']);
hold off;
end