Download Control System Design: Finding Limits and Steady State Errors and more Study notes Dynamics in PDF only on Docsity!

Professor Fearing EECS C128/ME C134 Problem Set 10 Fall 2011

1 Steady State Error (30 pts)

Given the following continuous time (CT) system

x ˙ = A 1 x + B 1 u =

[

− 12 − 7

]

x +

[

0

]

u(t), y = [10 0] x (1)

a) Given error e(t) = r(t) − y(t) where r(t) is a scalar, find limt→∞ e(t) for input r(t) a unit

step, and control law u(t) = r(t) − y(t).

Recall that

ess =

(

I + C A¯

− 1 ¯

B

)

rss

A^ ¯ = A − BKC

B^ ¯ = BK

when K is the controller gain. In this case, K = 1. Therefore:

A¯ =

[

− 12 − 7

]

−

[

0

]

[

]

=

[

− 12 − 7

]

−

[

0 0

]

=

[

− 22 − 7

]

A^ ¯−^1 = 1

22

[

− 7 − 1

]

B^ ¯ =

[

0

]

1 =

[

0

]

I + C A¯

− 1 ¯

B = 1 +

[

]

1 22

[

− 7 − 1

] [

0

]

= 1 −

5 11

ess =

6 11

rss =

6 11

b) Evaluate the steady state error limt→∞ e(t) for input r(t) a unit step, with state feedback,

that is, u = −K 1 x + r, where K 1 is chosen so that the closed loop poles are at si = − 5 , − 5.

First, find K 1. We want the eigenvalues of A 1 − B 1 K 1 to be (− 5 , −5). Find the characteristic

equation, then make it equal (s + 5)(s + 5):

A 1 − B 1 K 1 =

[

− 12 − 7

]

−

[

0

]

[

k 1 k 2

]

=

[

− 12 − k 1 − 7 − k 2

]

∆(s) =

∣ ∣ ∣ ∣ ∣ ∣

[

−s 1

− 12 − k 1 − 7 − k 2 − s

]

∣ ∣ ∣ ∣ ∣ ∣

= s(s + 7 + k 2 ) + k 1 + 12

= s

2

• (7 + k 2 )s + (k 1 + 12)

= (s + 5)(s + 5) = s

2

• 10s + 25

k 1 + 12 = 25

k 1 = 13

k 2 + 7 = 10

k 2 = 3

Now find the steady state error. Note that the method we used in (a) had u = K(r − y). This

time, u = r − K 1 x. The difference is that r does not get multiplied by K 1. We need to re-derive

the steady state error relationship.

x˙ = (A 1 − B 1 K 1 )x + B 1 r

Note there’s no C 1 , because we’re doing full state feedback. Assume that ˙x = 0 in steady state.

0 = (A 1 − B 1 K 1 )xss + B 1 rss

xss = −(A 1 − B 1 K 1 )

− 1 B 1 rss

ess = rss − C 1 xss

=

(

I + C 1 (A 1 − B 1 K 1 )

− 1 B 1

)

rss

Now we can evaluate ess.

A 1 − B 1 K 1 =

[

− 12 − 7

]

−

[

0

]

[

]

=

[

− 12 − 7

]

−

[

0 0

]

=

[

− 25 − 10

]

(actually, we could have jumped right to this, since A 1 , B 1 are in controllable canonical form, so

A 1 − B 1 K 1 will have the coefficients of the characteristic equation in the bottom row.)

(A 1 − B 1 K 1 )

− 1

1 25

[

− 10 − 1

]

I + C 1 (A 1 − B 1 K 1 )

− 1 B 1 = 1 +

[

]

1 25

[

− 10 − 1

] [

0

]

= 1 −

2 5

ess =

3 5

rss =

3 5

For parts c) and d) below, we want to add an integral controller to the system ( 1) to get zero

steady state error, but still have state feedback to allow arbitrary pole placement.

c) Let u = −K 2 x + ke

∫

(r − Cx)dτ. Write the augmented state and output equations for this system

with integral control, for example, z˙ = A 2 z + B 2 r, where A 2 is 3 × 3. Find gains K 2 , ke such that

the system has closed-loop poles at s = − 4 , − 5 , − 10 , and evaluate the steady state error for a step

input r(t).

We are adding a new state variable to represent the integrator. To avoid confusing it with the

original x state variables, we’ll call it w.

z =







w

x 1

x 2







Manipulate the definition of u to get dynamics for w.

u = −K 2 x + ke

∫

(r − C 1 x) dτ

u = −K 2 x + kew

w =

∫

(r − C 1 x) dτ

w ˙ = (r − C 1 x)

Notice that the new form, ˙z = A 2 z + B 2 r, doesn’t have u. We need to make A 2 capture the “inside

loop” of the full state feedback. Manipulate the original dynamics.

x˙ = A 1 x + B 1 u

= A 1 x + B 1 (−K 2 x + kew)

= (A 1 − B 1 K 2 )x + B 1 kew

Now we can write the state and output equations.

z˙ = A 2 z + B 2 r y = C 2 z

z ˙ =

[

0 −C 1

B 1 ke A 1 − B 1 K 2

] [

w

x

]

+







0

0





 r^ y^ =

[

0 C 1

]

[

w

x

]

z ˙ =







0 − 10 0

ke − 12 − k 1 − 7 − k 2













w

x 1

x 2





 +







0

0





 r^ y^ =

[

]







w

x 1

x 2







To find the gains K 2 , ke, find the characteristic equation for A 2 , and make it equal (s + 4)(s +

5)(s + 10):

∆(s) =

∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣







−s − 10 0

0 −s 1

ke − 12 − k 1 −s − 7 − k 2







∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣

= −s

∣ ∣ ∣ ∣ ∣ ∣

[

−s 1

− 12 − k 1 −s − 7 − k 2

]∣∣

∣

∣

∣

∣

− 10

∣ ∣ ∣ ∣ ∣ ∣

[

−s − 7 − k 2 ke

]∣∣

∣

∣

∣

∣

+ 0

∣ ∣ ∣ ∣ ∣ ∣

[

0 −s

ke − 12 − k 1

]∣∣

∣

∣

∣

∣

= −s(s(s + 7 + k 2 ) + 12 + k 1 ) − 10(ke)

= −s

3 − (7 + k 2 )s

2 − (12 + k 1 )s − 10 ke

We can flip the sign of ∆(s) for convenience, since we only care about its roots

= s

3

• (7 + k 2 )s

2

• (12 + k 1 )s + 10ke

= (s + 4)(s + 5)(s + 10) = s

3

• 19s

2

• 110s + 200

7 + k 2 = 19 ⇒ k 2 = 12

12 + k 1 = 110 ⇒ k 1 = 98

10 ke = 200 ⇒ ke = 20

Now we can evaluate the steady state error.

z˙ = A 2 z + B 2 r

0 = A 2 zss + B 2 rss

zss = −A

− 1 2 B^2 rss

yss = −C 2 A

− 1 2 B^2 rss

ess = rss − yss

=

(

I + C 2 A

− 1 2 B^2

)

rss

=









1 +

[

]







0 − 10 0

20 − 110 − 19







− 1 





0

0















= 0

d) Augment the system matrix A 1 to include an integrator (keeping the phase variable form),

and then choose u = K 3 (r · [1 0 0]

T − x). Write the augmented state and output equations for

this system for example, z˙ = A 3 z + B 3 r, where A 3 is 3 × 3. Find gains K 3 such that the system

has closed-loop poles at s = − 4 , − 5 , − 10 , and evaluate the steady state error for a step input r(t).

Here is the signal flow diagram of the original system. Notice that we start in phase variable

form.

u

∫ x^2 ∫ x^1 y

− 12

− 7

0

There is only one way to add an integrator to this system that keeps the phase variable form

and that does not modify the dynamics of the original system: after x 1.

u

∫ x^2 ∫ x^1 ∫ xINT y

− 12

− 7

0

0

10

0

Notice that this new integrator doesn’t actually affect anything (both gains leading away from

xINT are zero). That’s OK; when we do full state feedback, this new state will start affecting the

rest of the system.

We can write the pre-state-feedback system equations based on this diagram:

z =







xINT

x 1

x 2







z^ ˙ = A

′ 3 z^ +^ B

′ 3 u^ =







0 − 12 − 7





 z^ +







0

0





 u

y = C 3 z =

[

]

z

Next, we have to bring the “inside loop” into A 3 , like we did in part (c). To do this, we write

out the dynamics, then get rid of u by substitution.

z˙ = A

′ 3 z^ +^ B

′ 3 u

= A

′ 3 z^ +^ B

′ 3 K^3 (r^ ·^ [1^0 0]

T − x)

= (A

′ 3 −^ B

′ 3 K^3 )z^ +







0

0







[

k 1 k 2 k 3

]







0

0





 r

=







−k 1 − 12 − k 2 − 7 − k 3





 z^ +







0

0

k 1





 r

y =

[

]

z

Now we have A 3 , B 3 , and C 3. Next, we need to find the K 3 values. This time it’s easier, because

A 3 , B 3 are in phase variable form. We want the characteristic equation to equal (s + 4)(s + 5)(s +

10. = s

3

• 19s

2

• 110s + 200, so the bottom row of A 3 must be (− 200 , − 110 , −19).

−k 1 = − 200 ⇒ k 1 = 200

− 12 − k 2 = 1 − 10 ⇒ k 2 = 98

− 7 − k 3 = − 19 ⇒ k 3 = 12

Finally, we evaluate the steady state error. Because the system in (c) and (d) have the same form:

z˙ = Aiz + Bir: we can re-use the derivation of the steady state error from (c).

ess =

(

I + C 3 A

− 1 3 B^3

)

rss

=









1 +

[

]







− 200 − 110 − 19







− 1 





0

0















= 1

What’s going on here? The open loop system (defined by A

′ 3 , B

′ 3 ) clearly has a pole at the origin

(one of the elements of the bottom row of A

′ 3 is zero).^ Doesn’t that mean it’s a Type 1 system,

and therefore zero steady state error when driven by a unit step?

There are a couple of ways of looking at this. As mentioned earlier, that the new integrator isn’t

actually doing anything in the open loop form. It’s integrating x 1 , but the result of the integration

doesn’t affect the other states or the output.

Another way to look at it is to convert the open loop state space system into a transfer function.

(See section 3.6 of Nise). When you do this, you find that the TF representation of the open loop

system doesn’t have a pole at the origin.

2 Linear Quadratic Regulator (35 pts)

STOP

x

x

car 2 car 1 u2 u

−102 −100^0 t = 0 sec

safety zone

STOP

car 1

0

car 2

safetyzone

t = 10 sec

Consider a pair of cars in a “platoon” regulated to a stopping position, as illustrated above. The

dynamics of car 1 are (x 2 = ¨x 1 = u 1 and car 2 has a plant model x˙ 4 = ¨x 3 = 0. 5 u 2 where u 1 and u 2

are the car’s thrust due to engine and braking. (The offset of x 1 from car 1’s rear bumper prevents

the cars from colliding when x 3 = x 1 .) The outputs of the system are y 1 = x 1 and y 2 = x 3 − x 1.

Note that if y 2 > 0 then car 2 has intruded on the safety zone of car 1.

Initial conditions are car 1 at -100 m, 20 m/s, and car 2 at -102 m, 25 m/s.

a) Write the system equations in state space form.











x ˙ 1

x ˙ 2

x ˙ 3

x ˙ 4











=











0 0 0 0

0 0 0 0





















x 1

x 2

x 3

x 4











+











0 0

0 0











[

u 1

u 2

]

[

y 1

y 2

]

=

[

− 1 0 1 0

]











x 1

x 2

x 3

x 4











x(0) =











− 100

− 102











b) Use the LQR method (from Matlab function lqr(sys,Q,R), with Q=diag([2.5,10,2.5,10])

and R=diag([1,1]) to find an optimal K for the state feedback control u = −Kbx. Plot x(t) and

u(t) for the given initial condition (Matlab initial) and state feedback with gain Kb.

−

−

−

−

−

−

−

−

−

−

0

To: Out(1)

−4.5 0 1 2 3 4 5 6 7 8 9 10

−

−3.

−

−2.

−

−1.

−

−0.

0

To: Out(2)

Response to Initial Conditions

Time (seconds)

Amplitude

−100 0 1 2 3 4 5 6 7 8 9 10

−

−

−

−

0

20

t

u

Control effort

u 1 u 2

Kb =

[

1. 5811 3. 6280 0 0

0 0 1. 5811 4. 0404

]

Notice that the safety zone is violated by about 5cm at t = 7s.

c) In part b) both cars get within 1 m of the stop sign at a speed less than 1 m/s at 10 seconds.

However, car 2 lags more than 4 meters behind car 1 during the stopping trajectory. Find a new

cost function Q = diag([q 1 q 2 q 3 q 4 ]) which keeps car 2 within 2 m of car 1 for the whole trajectory,

while maintaining y 2 < 0 to prevent a collision, and distance from stop sign less than 1 m at 10

sec with velocity less than 1 m/s. (Suggestion: let q 1 = 2. 5 and q 2 = 10 and search for new q 3 and

q 4 .) Plot x(t) and u(t) for the given initial condition (Matlab initial) and state feedback with

new gain Kc.

There are many solutions that satisfy the constraints. The gray region in the following graph

shows q 3 , q 4 values that work.

6 7 8 9 10 11 12 13 14

20

25

30

35

40

45

50

55

60

Final

Safety

Distance

q

4

q 3

Taking Q=diag([2.5 10 10 40]) gives the following response:

−

−

−

−

−

−

−

−

−

−

0

To: Out(1)

−2 0 1 2 3 4 5 6 7 8 9 10

−1.

−1.

−1.

−1.

−

−0.

−0.

−0.

−0.

0

To: Out(2)

Response to Initial Conditions

Time (seconds)

Amplitude

−160 0 1 2 3 4 5 6 7 8 9 10

−

−

−

−

−

−

−

0

20

40

t

u

Control effort

u 1 u 2

Kc =

[

1. 5811 3. 6280 0 0

0 0 3. 1623 7. 2560

]

d) Find the solution to the Riccati equation P using Matlab function care(A,B,Q,R) and esti-

mate the cost J = (xT^ P x)(0) for each of b) and c).

Pb =











5. 7363 1. 5811 0 0

1. 5811 3. 6280 0 0

0 0 6. 3884 3. 1623

0 0 3. 1623 8. 0807

     x

T (0)Pbx(0) = 107880

Pc =











5. 7363 1. 5811 0 0

1. 5811 3. 6280 0 0

0 0 22. 9454 6. 3246

0 0 6. 3246 14. 5119

     x

T (0)Pcx(0) = 268030

Answers for Pc will vary depending on the choice of Q. Notice that the top left corner of Pc is

the same, since q 1 and q 2 were not changed.

e) Briefly compare the tradeoffs between control effort and time response between the two cases.

Notice that the time response, control effort, and K entries associated with the first car are

identical in (b) and (c). This will be true as long as the suggestion to keep q 1 and q 2 the same is

followed.

Car 2’s response, control effort, and K entries are different between the two cases. The control

effort used in (c) is bigger in magnitude than in (b). Correspondingly, the time response of y 2 (the

distance between the two cars) is much different as well: in (c), it seems to be decaying to zero

much faster than in (b). y 2 in (b) shows evidence of several modes with similar time constants at

work (the back-and-forth motion). In (c), y 2 seems to be dominated by a single, faster mode.

In general, we get better performance (faster convergence to zero) in (c) compared to (b), at

the cost of greater control effort.

3 Discrete Time Control (35 pts)

Given the following continuous time (CT) system

x ˙ = Ax + Bu =







− 8 − 12 − 5





 x^ +







0

0





 u(t),^ y^ = [2 1 0]^ x

the corresponding discrete time (DT) system is

x[n + 1] = Gx[n] + Hu[n] , y[n] = Cx[n]

which can be found using the Matlab function c2d(sys,T,’zoh’).

a) With initial condition x 0 = [1 0 0]

′ , plot the ZIR using Matlab function initial() for the

CT system and the DT system (with T = 0. 4 sec).

>> A = [0 1 0;0 0 1;-8 -12 -5];

>> B = [0;0;1];

>> C = [2 1 0];

>> D = 0;

ctsys = ss(A,B,C,D);

dtsys = c2d(ctsys,0.4,’zoh’);

initial(ctsys,[1;0;0],[0:.1:4]);

initial(dtsys,[1;0;0],[0:.4:4]);

Continuous time response:

Zero-Order Hold Discrete time response:

b) For output feedback control u = k(r − y), sketch the root locus for the equivalent transfer

function for the continuous time (CT) system.

[NUM,DEN] = ss2tf(A,B,C,D)

NUM =

0 -0.0000 1.0000 2.

DEN =

1.0000 5.0000 12.0000 8.

Root locus for TF:

s + 2

s

3

• 5s

2

• 12s + 8

=

s + 2

(s + 1)(s + 2 + 2j)(s + 2 − 2 j)

c) Determine the closed loop pole locations for the CT system for k = 5 and plot the closed-loop

step response using Matlab.

[z,p,k] = zpkdata(feedback(5*ss(A,B,C,D),1));

p{1}

ans =

-1.7285 + 2.9459i

-1.7285 - 2.9459i

-1.

step(feedback(5*ss(A,B,C,D),1),[0:.1:4]);

Continuous step response for system with feedback:

d) The closed loop DT system has state equation

x[n + 1] = (G − kHC)x[n] + kHr[n] , y[n] = Cx[n]

(which can be found using the Matlab feedback function). Using Matlab, determine the closed

loop pole locations for the DT system for k = 5 and sampling period T = 0. 4 sec and plot the step

response.

[z,p,k] = zpkdata(feedback(5*dtsys,1));

p{1}

ans =

0.

0.2579 + 0.5838i

0.2579 - 0.5838i

step(feedback(5*dtsys,1),[0:.4:4]);

Discrete step response for system with feedback:

e) Use Matlab (iteratively if necessary) to find a sampling period T which gives a closed-loop

step response that is “reasonably close” to the CT closed-loop step response. Determine closed-loop

pole locations, and plot the DT step response.

step(feedback(5*c2d(ctsys,0.4,’zoh’),1))

step(feedback(5*c2d(ctsys,0.3,’zoh’),1))

step(feedback(5*c2d(ctsys,0.2,’zoh’),1))

step(feedback(5*c2d(ctsys,0.1,’zoh’),1),[0:.1:4])

[z,p,k] = zpkdata(feedback(5*c2d(ctsys,0.1,’zoh’),1));

p{1}

ans =

0.

0.8159 + 0.2510i

0.8159 - 0.2510i

With sampling period of 0.1s, there is at least 8 samples for the first “hump” and the tail settles

quickly to the equilibrium value without jittering around.

Discrete step response for system with feedback (Sampling Period changed to T = 0 .1s)

f ) Briefly explain why the CT and DT ZIR responses from a) above are reasonably close, but

the closed loop responses from c) and d) (with T = 0. 4 sec) do not agree at all. (Hint, consider

e

AT .)

We know there will be some error in the open-loop responses due to discritization error. Intu-

itively, controlling the system based on an erroneous output would cause even more error in the

responses.

Lets look at a quantitative analysis of the differences:

For the discrete system:

G = e

A· 0. 4







0. 9491 0. 3180 0. 0392

− 0. 3138 0. 4785 0. 1219

− 0. 9755 − 1. 7771 − 0. 1312







and, in feedback,

G − 5 HC =







0. 8856 0. 2863 0. 0392

− 0. 7060 0. 2824 0. 1219

− 2. 1945 − 2. 3866 − 0. 1312







Suppose we discretized the closed loop continous system instead so that the sampler is outside

the controller:

Gcl = e

(A− 5 BC)· 0. 4







0. 8899 0. 2836 0. 0366

− 0. 6585 0. 2680 0. 1006

− 1. 8117 − 2. 3696 − 0. 2353







Note how these two matrices are similar, but there are a few entries that are off by a few tenths.

In the first case, the controller is acting on observations containing discretization error, wheras

in the second case, the controller is acting continously, and then the system is discretized. The

difference in the matrices show that a controller can, in effect, amplify discretization error when

modelling a continuous system.