Example H6.14

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

Rootlocus and determination of instability gain for a feedback system with delay

Date : 14/07/2018
Revision : 1.0

Table of Contents

Rootlocus and determination of instability gain for a feedback system with delay The consequence of time delay on stability: the assignment Solution MATLAB code for this example Create models: Analysis without delay Analysis with delay Conclusion Rootlocus of the system with Pade approximation Show the step response of the found gain value where is oscillates Generate MATLAB figure(s) for usage in the book Init create Enhanced Figures Enhance the figures: Figure 1

The consequence of time delay on stability: the assignment

Figure 6.19a shows a heat exchanger. This process has the following properties:

the time constant of the heat transfer (first order process) is 50 seconds
the time constant of the valve control (first order process) is 10 seconds
the time delay of the liquid is 6 seconds
feedback is applied from sensor to valve with gain K

This creates a control system according to the block diagram of figure 6.19b.

The assignment is to construct the root locus for

and to determine the value of K for which the system becomes unstable.

Solution

The following transfer function applies to the process:

As the root locus in figure 6.20 shows, without time delay (dashed line) that system could never become unstable.

When the pure time delay of 6 seconds is replaced by a second order Padé approximation, the transfer function becomes:

The Padé appoximation yields two complex adjugated poles:

and two complex adjugated zeros

. In figure 6.20 these poles and zeros are added to those of the heat exchanger and the root locus has been constructed using the computer (solid line). If we compare the solid line with the dashed line, we see for instance that for an often desired

(this means

) the influence of the time delay is quite small. But stability is strongly influenced by the time delay. For the intersection with the imaginary axis we see:

rad/s. If we determine the gain

for this point using the graphical method according to formula 6.46 we find:

The exact values can be determined by simulation and turn out to be:

rad/s and

. Especially the deviation in the value of

seems quite large, but for a practical setting with

we find a value for K of approximately 4.

MATLAB code for this example

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

Create models:

First we create the models for the proces and the valve.

% Create proces model
Tau_heat_transfer=50;   % Time constant heat transfer
Tau_valve=10;           % Time constant of valve control
Td=6;                   % Transporttime fluid
% create model 
H_valve=1/(Tau_valve*s+1);
H_1=exp(-s*Td)*(1/(Tau_heat_transfer*s+1));
% Create total open loop model
H_tot=H_valve*H_1;

Then we create a model which uses Pade approximation of the delay.

% Create proces model with pade approximation of delay
Order=2;
H_1_pade=pade(exp(-s*Td),Order)*(1/(Tau_heat_transfer*s+1));
H_tot_pade=H_valve*H_1_pade;
 

Analysis without delay

We are assuming that the temperature sensor works as a unit feedback and that there is gain controller that actuates the valve. See figure 6.19b of the book.

First we look at the root locus of the systeem if the delay would not be present:

H_1_nodelay=H_valve*(1/(Tau_heat_transfer*s+1));

figure(101);

rlocusplot(H_1_nodelay);

The rootlocus shows that this system without delay would always be stable (for every positive value of K') as the rootlocus does not cross the vertical axis.

Analysis with delay

Next we look at the system with delay present. Without an approximation of the delay the rootlocus can not be shown. But it is offcourse possible to show the bode diagram of the open system to find the gain margin at a phase shift of -180 degrees.

First we will show the bode plot:

figure(102);

bodeplot(H_tot,H_tot_pade);

grid on;

title('Bode diagram of H_{tot}');

legend('No Pade','With Pade','Location', 'SouthWest')

Zooming in around the area of -180 deg phase shift gives:

figure3=figure(103);
ZoomArea={0.1, 0.2};
bodeplot(H_tot,H_tot_pade,ZoomArea);
grid on;
title('Bode diagram of H_{tot} zoomed in around -180 deg');
legend('No Pade','With Pade','Location', 'SouthWest')
% Place textarrow 
annotation(figure3,'textarrow',[0.432 0.514],[0.693 0.693],...
    'String',{'Gain at 180 degrees'});

You can notice that the pade has a phase difference of 360 degrees. For instablity this can be disregarded. You can now read the gain margins from the plot. We could also use a MATLAB command for this:

%%Determine the gain margin of H_tot from the bode plot;

[Gm,~,Wgm,~] = margin(H_tot); % [Gm,Pm,Wgm,Wpm] but Pm and Wpm not used thus ~

MsgStr=['Gain margin H_{tot} is : ', num2str(Gm), ' at ', num2str(Wgm), ' rad/s'];

% Add the answer as text to plot

annotation(figure3,'textbox',...

[0.17 0.53 0.39 0.064],...

'String',{MsgStr},...

'LineStyle','none',...

'FitBoxToText','on');

%%Determine the gain margin of H_tot_pade from the bode plot;

[Gm,~,Wgm,~] = margin(H_tot_pade); % [Gm,Pm,Wgm,Wpm] but Pm and Wpm not used thus ~

MsgStr=['Gain margin H_{tot\_pade} is : ', num2str(Gm), ' at ', num2str(Wgm), ' rad/s'];

% Add the answer as text to plot

annotation(figure3,'textbox',...

[0.17 0.48 0.39 0.064],...

'String',{MsgStr},...

'LineStyle','none',...

'FitBoxToText','on');

Conclusion

It shows that the difference between found value of the gain for oscillation without and with the 2nd order Pade approximation is very small. It can be concluded that the second order Pade approximation of the delay works very well in this case.

Rootlocus of the system with Pade approximation

Next we look again at the system when we approximate the delay with Pade. This allows us to also show the rootlocus.

figure(104);

h=rlocusplot(H_tot_pade);

grid on;

% get the handle to the plot options

p=getoptions(h);

% Zoom in

p.XLim=[-1 1];

p.YLim=[-1 1];

% Adjust Title

p.Title.String='Root Locus of H_{tot\_Pade}';

% write the options to the figure

setoptions(h,p);

Notice that the poles and zeros of the pade approximation are symmetrical around the imaginary axis with the zeros on the right and the poles on the left.

Show the step response of the found gain value where is oscillates

As a last step we will verify the oscillation using the step response. According to the bode diagrams, ocillation occurs around K_r=11.14. Here we take a look at the stepresponse of the feedback system for that gain.

K_r=11.14;

Tsim=1000;

H_cl=feedback(K_r*H_tot,1);

figure(105);

step(H_cl,Tsim);

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: Figure 1

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(H_tot_pade,'k',H_1_nodelay,'k--');

% Determine Kpn

[Z,P,Kpn]=zpkdata(H_tot_pade);

% 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_{tot\_pade}';

% 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');