Section 36-04.html

Chapter VII: Complex Variables

Unit 2: Applications

Section 36.4: the root locus

Copyright

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> with(linalg):
with(plots):
with(plottools):
with(inttrans):
with(student):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

Warning, the name hilbert has been redefined

> alias(X=laplace(x(t),t,s), Y=laplace(y(t),t,s), F=laplace(f(t),t,s)):

>

Introduction

In the study and design of feedback control systems, engineers make use of the root-locus diagram which shows the variation of eigenvalues as a function of a system parameter. The eigenvalues are actually found as poles of an appropriate transfer function, a function of the complex variable arising from the Laplace transform of the differential equations of the system. For each value of the system parameter, the poles are plotted in the complex plane. The resulting locus of points in the complex plane is called the root locus .

This section will give some mathematical insight into the construction of a root-locus diagram. We begin with a study of the "locus of roots," and reserve the phrase "root locus" specifically for the case pertinent to controls engineering. (For a more complete engineering perspective, see [1] or [2].)

>

Example 36.20

In Chapter 12, the mixing-tank problem led to x ' = A x , the first-order linear system of ordinary differential equations with constant coefficients. The solution is given by x , where x = x (0), and is the fundamental matrix. The eigenvalues of determine stability of the equilibrium point at the origin.

Suppose the matrix A contains the parameter , making the system performance dependent on the value of . For example, let A be the matrix

> A := matrix(2,2,[1,c,3,4]);

>

The eigenvalues of A are then functions of , and are given by

> q := eigenvals(A):
Q := simplify(subs(c=0,[q])):
if Q[1]>Q[2] then
lambda[1] := q[1]:
lambda[2] := q[2]:
else lambda[1]:=q[2]:lambda[2]:=q[1];fi;

>

(The additional Maple commands are to guarantee the eigenvalues are consistently labeled.)

These eigenvalues are graphed as functions of in Figure 36.7.

> p1 := plot([lambda[1],lambda[2]], c=-1..2, color=[black,red], linestyle=[1,2], labels=[c,``], labelfont=[TIMES,ITALIC,12], xtickmarks=4, ytickmarks=6):
p2 := textplot({[.2,5.2,`l`],[1.8,5,`l+`],[.7,.7,`l-`]}, font=[SYMBOL,12]):
display([p||(1..2)], scaling=constrained);

>

The solid (black) curve is , whereas the dotted (red) curve is . For , the eigenvalues are complex with positive real part, so the equilibrium point at the origin is unstable. For , there are two positive real eigenvalues, the origin is a node (out) and is again unstable. For , one eigenvalue is poistive and the other is negative, so the origin is a saddle, again unstable.

We have just seen that when causes both eigenvalues to fall into the left-half of the complex plane, the origin is an asymptotically stable equilibrium point, but when causes at least one eigenvalue to fall into the right-half of the complex plane, the origin is an unstable equilibrium point. A plot of the eigenvalues as points in the complex plane is an important tool in the design of feedback control systems which also can be modeled with systems of linear differential equations. In the complex -plane, the curve traced by the eigenvalues as they vary with an appropriate parameter of the system, is called a root-locus . Figure 36.8 shows, for this example, the locus of roots as varies in the interval . (Our parameter is not one which would be found in a typical problem in controls engineering.)

> p1 := complexplot(lambda[1], c=-2..4, color=black, thickness=3):
p2 := complexplot(lambda[2], c=-2..4, color=cyan, thickness=3):
p3 := textplot({[2.5,2,`c = -2`], [2.5,-2.1,c = -2], [-1,.3,`c = 2`], [6,.3,`c = 2`]}):
p4 := textplot({[2.7,1.4,`l+`],[2.7,-1.4,`l-`]},font=[SYMBOL,12], align=RIGHT):
display([p||(1..4)],scaling=constrained,xtickmarks=5,ytickmarks=4, labels=[` Re`,`Im `]);

>

The initial and terminal points on each locus are indicated on the graph. The black curve is the locus of , whereas the cyan curve is the locus of . When , both eigenvalues are complex. In fact, they are complex conjugates. As increases towards , the two eigenvalues approach the real axis. For , the eigenvalues remain along the positive real axis. For , remains on the positive real axis, but is now found on the negative real axis.

The following animation shows the dynamic development of these loci as increases.

> f1 := t -> plot([Re(lambda[1]), Im(lambda[1]), c=-2..t], color=black, thickness=3):
f2 := t -> plot([Re(lambda[2]), Im(lambda[2]), c=-2..t], color=cyan, thickness=3):
p4 := display([seq(f1(k/5),k=-10..10)],insequence=true):
p5 := display([seq(f2(k/5),k=-10..10)],insequence=true):
p6 := disk([2.5,2],.1,color=black):
p7 := disk([2.5,-2],.1,color=cyan):
p8 := textplot({[2.9,2,`c = -2`],[2.9,-2,`c = -2`]}):
display([p||(4..8)],scaling=constrained, labels=[`Re`,`Im`], xtickmarks=6, ytickmarks=5);

>

Example 36.21

The first-order linear system x ' = A x for which the matrix A is given by

> A := matrix(2,2,[c,-5,3,4]);

>

has eigenvalues

> q := eigenvals(A):
Q:=simplify(subs(c=20,[q])):
if Q[1]>Q[2] then
lambda[1] := q[1]:
lambda[2] := q[2]:
else lambda[1]:=q[2]:lambda[2]:=q[1];fi;

>

Figure 36.9 contains a graph of these eigenvalues as functions of the parameter with shown in black, and shown in cyan.

> pp1 := plot([lambda[1],lambda[2]], c=-10..20, color=[black,cyan], linestyle=[1,2], numpoints=500):
pp2 := textplot({[1,17,`l`],[18,14.5,`l+`],[18,4,`l-`], [-8,4,`l+`],[-8,-9,`l-`]},font=[SYMBOL,12]):
display([pp||(1..2)], scaling=constrained,labels=[c,``], labelfont=[TIMES,ITALIC,12], xtickmarks=7, ytickmarks=5);

>

There are definitely values of for which the eigenvalues are of opposite sign. There are also values of for which the eigenvalues are complex. It is not yet clear if there are values of for which both eigenvalues remain positive. Thus, we compute

that is,

> Limit(lambda[2],c=infinity) = limit(lambda[2],c=infinity);

>

to determine if, for large , the eigenvalue remains positive. Hence, for large enough, both eigenvalues are positive.

To determine the specific ranges of for which each outcome occurs, obtain the characteristic polynomial for A . Thus, obtain

> char_poly := charpoly(A,s);

>

for the characteristic polynomial, and

> discriminant := discrim(char_poly,s);

>

for its discriminant. The transition between real and complex takes place where the discriminant vanishes, so solve the equation

> q1 := discriminant = 0;

>

for + as the values for transition between real and complex eigenvalues.

Using Maple, these values are found in exact and floating-point form by

> solve(q1,c);
fsolve(q1,c);

>

Figure 36.10 shows the locus of roots for .

> p1 := complexplot(lambda[1], c=-5..15, color=black, thickness=3):
p2 := complexplot(lambda[2], c=-5..15, color=cyan, thickness=3):
p3:=textplot({[-3,.3,`c = -5`], [2,.3,c = -5], [6,.3,`c = 15`], [13.5,.3,`c = 15`]}):
p4:=textplot({[7,3.1,`l+`],[7,-3.1,`l-`]},font=[SYMBOL,12],align=RIGHT):
display([p||(1..4)], scaling=constrained, xtickmarks=5, ytickmarks=4, labels=[` Re`,`Im `]);

>

For in the range , (in black) is along the positive real axis, but (in cyan) is along the negative real axis. For in the range , the eigenvalues are complex conjugates. The locus of roots consists of the upper and lower halves of a circle with radius 4 and center at . For , both eigenvalues are on the positive real axis. We have already seen that tends towards the point . The locus of roots indicates that as , the eigenvalue also tends to infinity along the positive real axis.

The following animation shows the dynamic development of this locus of roots.

> f1 := t -> plot([Re(lambda[1]),Im(lambda[1]),c=-5..t], color=black, thickness=3):
f2 := t -> plot([Re(lambda[2]),Im(lambda[2]),c=-5..t], color=cyan, thickness=3):
p4 := display([seq(f1(k/2),k=-10..30)],insequence=true):
p5 := display([seq(f2(k/2),k=-10..30)],insequence=true):
P1 := subs(c=-5,lambda[1]):
P2 := subs(c=-5,lambda[2]):
p6 := disk([P1,0],.2,color=black):
p7 := disk([P2,0],.2,color=cyan):
display([p||(4..7)],scaling=constrained, xtickmarks=8, ytickmarks=5);

>

Transfer Functions and Block Diagrams

Recall that for a differential equation of the form

that is,

> q := diff(y(t),t,t)+2*diff(y(t),t) + 10*y(t) = f(t);

>

with the inert initial conditions = ' , the Laplace transform leads to

that is,

> q1 := subs(y(0)=0,D(y)(0)=0,laplace(q,t,s));

>

and

that is

> q2 := isolate(q1,Y);

>

The function is an "input" to the system, providing the system with a driving force or excitation function. The solution is the system "output," so the function

that is,

> T = coeff(rhs(q2),F);

>

is called the transfer function for the differential equation. Since , the transfer function is characterized as the ratio of Laplace transforms of the system output over the system input. (See Section 6.11 for a determination of the transfer function by means of the Dirac delta function.) Written as , the transfer function determines how the action of the input, , will be transferred into the output, .

The engineering approach to feedback control systems is couched in the language of transfer functions and block diagrams. To illustrate this approach, consider the first-order linear system on the left in Table 36.2.

______________________________________________________________________

Table 36.2

The system parameter is often called a gain in the controls literature. The Laplace transform of each equation, in concert with the inert initial conditions = 0 leads to the equations in the center of Table 36.2. Solving the first of these equations for , and the second for gives the equations on the right in Table 36.2, where

are transfer functions.

This form of the differential equations, namely,

> q3 := diff(x(t),t) = x(t)+2*y(t);
q4 := diff(y(t),t) = -K*x(t)+4*y(t)+3*f(t);

>

is modeled by the block diagram in Figure 36.11.

> p1 := arrow([0,0],[1,0],.1,.3,.2,color=black):
p2 := rectangle([1,.5],[2,-.5]):
p3 := arrow([2,0],[3.25,0],.1,.3,.2,color=black):
p4 := circle([3.5,0],.25):
p5 := arrow([3.75,0],[5,0],.1,.3,.2,color=black):
p6 := rectangle([5,.5],[6,-.5]):
p7 := arrow([6,0],[8,0],.1,.3,.2,color=black):
p8 := rectangle([5,-1.5],[6,-2.5]):
p9 := arrow([7,0],[7,-2],.1,.3,.2,color=black):
p10 := arrow([7,-2],[6,-2],.1,.3,.2,color=black):
p11 := arrow([5,-2],[3.5,-2],.1,.3,.2,color=black):
p12 := arrow([3.5,-2],[3.5,-.25],.1,.3,.2,color=black):
p13 := textplot({[1.5,0,`G3`],[5.5,0,`G1`],[5.5,-2,`K G2`], [.2,.2,`F`], [4.25,.2,`Y`], [7.3,.2,X], [3.25,-.45,`-`], [3.15,.3,`+`]}, font=[TIMES,BOLD,10]):
display([p||(1..13)], scaling=constrained, axes=none);

>

The circle is a junction node, the squares are "blocks" representing the transfer functions, and the secret to reading the diagram is to realize that represents a system input , and , a system output . The loop containing the blocks representing and is the feedback loop. The surprise is that inhomogeneous first-order systems arising from mixing-tank problems, and studied in Chapter 12, can be interpreted as feedback systems. Hence, feedback is inherently found in the coupling of the differential equations, and not necessarily in an array of sensors and amplifiers.

The Laplace transform of each equation, in concert with the inert initial conditions

> inits := {x(0)=0,y(0)=0};

>

leads to

> q5 := subs(inits,laplace(q3,t,s));
q6 := subs(inits,laplace(q4,t,s));

>

Solve the first equation for , and the second, for , thereby obtaining

> q7 := expand(isolate(q5,X));
q8 := isolate(q6,Y);

>

Define the transfer functions

then rewrite the transformed differential equations in terms of these transfer functions, obtaining

> q9 := X = G[1]*Y;
q10 := Y = -K*G[2]*X+G[3]*F;

>

Finally, solve for and in terms of the transfer functions, obtaining

> q11 := solve({q9,q10},{X,Y});

>

The transfer function connecting the input to the output is then

The poles of are values of for which the denominator of vanishes. The equation whose roots are the poles of is the characteristic equation for the system matrix

>

> A := genmatrix([rhs(q3),rhs(q4)],[x(t),y(t)]);

>

In fact, the characteristic polynomial for this matrix is

> charpoly(A,s);

>

and the eigenvalues of A are

> eigenvals(A);

>

The denominator of is

> q12 := 1+2*K/(s-1)/(s-4);

>

which, if written with a common denominator, is

> normal(q12);

>

The denominator of vanishes when the characteristic polynomial of A vanishes! Thus, the poles of the transfer function are

> solve(q12 = 0,s);

>

which are precisely the eigenvalues of A . Thus, the stability of the system is determined by the location of the poles of the transfer function because these poles are the eigenvalues of the system matrix.

The prevalence of the form in the denominator of the transfer function leads to the Maple rootlocus command which draws the root locus for solutions of the equation . The inputs to the command are and a range for . The use of this command, and the mathematics it takes to draw the root locus, will be explored in the following examples.

>

Example 36.22

Suppose the denominator of a transfer function for a feedback system is , where is given by

> h1 := (s+1)/s^2;

>

Figure 36.12, created by Maple's rootlocus command, provides a graph of that part of the root locus corresponding to the range provided for the gain . Thus, we have the graph

> rootlocus(h1, s, 0..5, scaling=constrained,thickness=3, xtickmarks=4, ytickmarks=3, labels=[` Re`,`Im `]);

>

In any design process based on this root-locus diagram, it would be important to know where in the complex s-plane the roots (eigenvalues) were for various values of . There are a variety of ways to obtain this information.

First, form the denominator of the transfer function and set it equal to zero, obtaining

> unassign('k'):
q := 1 + k*h1 = 0;

>

Adding fractions leads to

> normal(q);

>

so the equation is equivalent to the vanishing of the numerator. Thus,

> q1 := numer(normal(lhs(q))) = 0;

>

is actually the characteristic equation for the system.

The roots as a function of , that is, , would then be given by

> q2 := solve(q,s);

>

The characteristic equation is quadratic in s so there are two roots. However, each root depends on so each root is a function of . As we did earlier, we could animate the tracing of the root locus, obtaining the following animation.

> f1 := t -> plot([Re(q2[1]), Im(q2[1]), k=-1..t], color=black, thickness=3):
f2 := t -> plot([Re(q2[2]), Im(q2[2]), k=-1..t], color=cyan, thickness=3):
p4 := display([seq(f1(n/2),n=-2..10)],insequence=true):
p5 := display([seq(f2(n/2),n=-2..10)],insequence=true):
p6 := disk([-.62,0],.1,color=cyan):
p7 := disk([1.62,0],.1,color=black):
p8 := textplot({[-.62,.2,`k = -1`],[1.62,.2,`k = -1`]}):
display([p||(4..8)],scaling=constrained, labels=[`Re`,`Im`], xtickmarks=5, ytickmarks=3);

>

Another approach to studying the dependence of the characteristic roots on the gain is to solve the equation

> q;

>

for rather than . This gives

> K := solve(q,k);

>

But in the complex plane, is the complex variable . What about ? For this, we get

> q3 := subs(s = x+I*y, K);

>

What do we know about ? Only that is real. Hence, split into its real and imaginary parts and set the imaginary part equal to zero.

> q4 := simplify(evalc(Re(q3)));
q5 := simplify(evalc(Im(q3)));

>

The equation

= 0

implicitly defines a curve y = y(x). In this example we can obtain this curve explicitly.

> solve(q5,y);

>

However, the vanishing of the numerator leads to the equation

> completesquare(numer(q5)/(-y),x) = 0;

>

which is the implicit form of a circle with radius 1 and center at .

The standard approach to constructing a root locus is numeric. For a fixed value of the gain , the equation

is solved numerically for the two corresponding characteristic roots . For example, if , the corresponding roots are

> fsolve(subs(k=1,q1),s,complex);

>

If this process is repeated for a sequence of values of , a collection of points in the complex -plane will result. For example, for , the collection of solutions would be

> q6 := [seq(fsolve(subs(k=n/10,q1),s,complex),n=0..50)];

>

A plot of these points would then give an outline of the root locus. Here, Maple's complexplot command from the plots package will plot a list of complex numbers as points in the complex plane, giving the following graph.

> complexplot(q6, style=point, scaling=constrained);

>

This is, in fact, the algorithm used by Maple's rootlocus command, and by the equivalent functionality in other commercial packages which provide the root locus. Maple's rootlocus command also contains heuristics for joining the points together into curves.

>

Example 36.23

Let the function be

> h2 := (s+1)/s^3;

>

The root locus (Figure 36.13), that is, the curve defined implicitly by the equation

can be obtained by

> rootlocus(h2, s, 0..5, xtickmarks=2, scaling=constrained, labels=[` Re`,`Im `], thickness=3);

>

Maple's rootlocus command draws the root locus by numerically computing solutions for a range of values of . To see this, we first obtain the characteristic equation. The denominator of the transfer function is

> q7 := 1+k*h2;

>

and the characteristic equation is

> q8 := numer(normal(q7)) = 0;

>

Points on the root locus are obtained by fixing at some value, and computing the roots . For example, if , then the three roots of the characteristic equation would be

> fsolve(subs(k=.1,q8),s,complex);

>

Computing and plotting in the complex -plane a large number of such points leads to the root-locus diagram.

Alternatively, since the characteristic equation is just a cubic in , we can obtain the analytic solution

> q9 := solve(q8,s);

>

There are three branches, each of which can be plotted in the complex plane via Maple's complexplot command. The result is the following version of the root locus, where the three branches are shown in black, red, and green.

> S:=[seq(n*.0011,n=1..100),seq(n*.1,n=1..50)]:
p1 := complexplot(q9[1],k=0..5,color=black, thickness=3, sample=S):
p2 := complexplot(q9[2],k=0..5,color=red, thickness=3, sample=S):
p3 := complexplot(q9[3],k=0..5,color=green, thickness=3, sample=S):
display([p||(1..3)], scaling=constrained, xtickmarks=[-1,0,1], view=[-1..1,-3..3]);

>

For , the roots are all , so the three branches of the root locus start at the origin, and become unbounded as increases.

>

References

[1] Charles L. Phillips and Royce D. Harbor, Basic Feedback Control Systems, Alternate Second Edition, Prentice-Hall, Inc., 1991.

[2] Gene F. Franklin, J. David Powell and Abbas Emami-Naeini, Feedback Control of Dynamic Systems, Addison-Wesley Publishing Company, 1986.

>