Introduction To Control Theory Including Optimal Control

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

1. System Dynamics and Differential Equations 1.1 Introduction

1.4 shows a tank with: inflow rate q i , outflow rate q o , head

- Investigation of dynamics behavior of idealized components, and the collection of interacting components or systems are equivalent to studying differential equations. - Dynamics behavior of systems ⇔ differential equations.

level h , and the cross-section are of the tank is A . inflow

1.2 Some System Equations h (head)

Inductance circuit: In electromagnetic theory it is known that, for a coil, in Fig.1.1, having an inductance L , the electromotive force (e.m.f.) E is proportional to the rate of change of the current I , at the instant considered, that is, E=L

dI dt

or I =

∫

1 E dt L

(1.1) L

valve Fig. 1.4

If q = q o − q i is the net rate of inflow into a tank, over a period δ t , and δ h is the corresponding change in the head level, then q δ t = A δ h and in the limit as δ t → 0 q=A

E Fig. 1.1

Capacitance circuit: Similarly, from electrostatics theory we know, in Fig. 1.2, that the voltage V and current I through a capacitor of capacitance C are related at the instant t by

∫ I dt

dE or I = C dt

or h =

1 A

∫ q dt

(1.4)

dy =x dt

(1.5)

(1.2)

Equation (1.5) is interesting because: - its solution is also the solution to any one of the systems considered. - it shows the direct analogies which can be formulate between quite different types of components and systems. - it has very important implications in mathematical modeling because solving a differential equation leads to the solution of a vast number of problems in different disciplines, all of which are modeled by the same equation.

C

E Fig. 1.2

Dashpot device: Fig. 1.3 illustrates a dashpot device which consists of a piston sliding in an oil filled cylinder. The motion of the piston relate to the cylinder is resisted by the oil, and this viscous drag can be assumed to be proportional to the velocity of the piston. If the applied force is f (t ) and the corresponding displacement is y (t ) , then the Newton’s Law of motion is

f =µ

dh dt

All these equations (1.1)-(1.4) have something in common: they can all be written in the form a

1 E= C

discharge q0

dy dt

(1.3)

where the mass of the piston is considered negligible, and µ is the viscous damping coefficient. y (t )

Over small ranges of the variables involved the loss of accuracy may be very small and the simplification of calculations may be great. It is known that the flow rate through a restriction such as a discharge valve is of the form q =V

P

where P is the pressure across the valve and V is coefficient dependent on the properties of the liquid and the geometry of the valve. Fig.1.5 shows the relation between the pressure and the flow rate. q flow rate q1

assumed linear

f (t )

Fig.1.3 Liquid level: To analyze the control of the level of a liquid in a tank we must consider the input and output regulations. Fig.

P1

P (pressure)

Fig. 1.5

___________________________________________________________________________________________________________ 1 Chapter 1 System Dynamics and Differential Equation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Assume that in the neighborhood of the pressure P = P1

We can also write the set of simultaneous first order equation (1.7) and (1.8) as one differential equation of the second order. On differentiating (1.7), we obtain

change in pressure ≈ R' change in flow rate

h&&2 =

where R' is a constant called the resistance of the valve at the point considered. This type of assumption, which assumes that an inherently nonlinear situation can be approximated by a linear one albeit in a restricted range of the variable, is fundamental to many applications of control theory.

1

A 1 1 ⎛ 1 1 ⎞ 1 ⎜ ⎟ h2 + + q − 2 h&2 − h2 A1 A1 A1 ⎜⎝ R1 R 2 ⎟⎠ A1 R1 A 1 1 = q − 2 h&2 − h2 A1 A1 A1 R 2 Then h&&2 = =

2

A2

A1 h1

R1

h2

⎡ 1 1 q−⎢ A1 A2 R1 ⎢⎣ A2

⎡ 1 h&&2 + ⎢ ⎣⎢ A2

⎛ 1 1 ⎜ ⎜R + R 2 ⎝ 1

dh1 1 = q− (h1 − h2 ) dt R1

⎞ 1 ⎤& 1 ⎟+ ⎟ A R ⎥ h2 − A A R R h2 1 1 ⎥⎦ 1 2 1 2 ⎠

⎛ 1 1 ⎜ ⎜R +R 2 ⎝ 1

Consider a control system: a tank with an inflow q i and a

(1.7)

discharge q o in Fig. 1.7

For tank 2

⎛ 1 1 ⎜ ⎜R + R 2 ⎝ 1

float

lever

dh 1 1 (h1 − h2 ) − h2 A2 2 = dt R1 R2

⇒

⎞& ⎟ h2 ⎟ ⎠

1.3 System Control

dh 1 1 1 h&1 = 1 = q− h1 + h2 dt A1 A1 R1 A1 R1

dh 1 1 h&2 = 2 = h1 − dt A2 R1 A2

⎛ 1 1 ⎜ ⎜R + R 2 ⎝ 1

⎞ 1 ⎤& 1 1 ⎟+ ⎟ A R ⎥ h 2 + A A R R h2 = A A R q 1 1 ⎦⎥ 1 2 1 2 1 2 1 ⎠ (1.11) Equation (1.11) is a differential equation of order 2 of the form &y& + a1 y& + a 2 y = b0 x .

R2

For tank 1

⇒

1 1 & 1 1 q− h2 − h2 − A1 A2 R1 A1 R1 A1 A2 R1 R2 A2

or

Fig. 1.6

A1

(1.10)

=

inflow q

⎞& ⎟ h2 ⎟ ⎠

1 1 1 h&1 = q− h1 + h2 A1 A1 R1 A1 R1

(1.6)

A two tank systems with one inflow (q) and two discharge valves is shown in Fig. 1.6.

⎛ 1 1 ⎜ ⎜R + R 2 ⎝ 1

From (1.7) and (1.8),

At a given point, the pressure P = hρ g , define R = R' /( ρ g ) , we can rewrite the above relation as R = h / qo

1 & 1 h1 − A2 R1 A2

⎞ ⎟ h2 ⎟ ⎠

R

(1.8)

inflow

outflow regulator valve Fig. 1.7

We can write (1.7) and (1.8) in matrix form as follows 1 ⎡ − ⎡ h&1 ⎤ ⎢⎢ A1 R1 ⎢& ⎥ = ⎢ 1 1 ⎣ h2 ⎦ − ⎢ A R A2 ⎣ 2 1

1 ⎤ ⎡ 1 ⎤ ⎥ A1 R1 ⎥ ⎡ h1 ⎤ + ⎢ A ⎥ q ⎢ ⎥ ⎞ ⎛ 1 h 1 ⎥⎣ 2⎦ ⎢ 1⎥ ⎟ ⎜ ⎣ 0 ⎦ ⎜ R + R ⎟⎥ 2 ⎠⎦ ⎝ 1

that is, in the form

h& = A h + B q

(1.9)

where 1 ⎡ ⎢− A R ⎡ h&1 ⎤ 1 1 h= ⎢& ⎥ , A= ⎢ 1 1 ⎢ h ⎣ 2⎦ − ⎢ A R A 2 1 2 ⎣

1 ⎤ ⎡1 ⎤ ⎥ A1 R1 ⎥, B=⎢A ⎥ ⎛ 1 ⎢ 1⎥ 1 ⎞⎥ ⎟ ⎜ ⎣ 0 ⎦ ⎜ R + R ⎟⎥ 2 ⎠⎦ ⎝ 1

The purpose of controller is to maintain the desired level h . The control system works very simply. Suppose for some reason, the level of the liquid in the tank rises above h . The float senses this rise, and communicates it via the lever to the regulating valve which reduces or stops the inflow rate. This situation will continue until the level of the liquid again stabilizes at its steady state h . The above illustration is a simple example of a feedback control system. We can illustrate the situation schematically by use of a block diagram as in Fig. 1.8.

ε

h1

x

regulator

plant

h2

h2

Fig. 1.8 ___________________________________________________________________________________________________________ 2 Chapter 1 System Dynamics and Differential Equation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

theory is to design – in terms of the transfer function – a system which satisfies certain assigned specifications. This objective is primarily achieved by a trial and error approach.

with h1

: desired value of the head

h2

: actual value of head

ε = h1 − h2 : error The regulator receives the signal ε which transforms into a movement x of the lever which in turn influences the position of the regulator valve. To obtain the dynamics characteristics of this system we consider a deviation ε from the desired level h1 over a short time interval dt . If q i and q o are the change in the inflow and outflow, from (1.4),

A

dε = qi − q o dt

where (1.6), q o = ε / R . So that the equation can be written as AR

Modern control theory is not only applicable to linear autonomous systems but also to time-varying systems and it is useful when dealing with nonlinear systems. In particularly it is applicable to MIMO systems – in contrast to classical control theory. The approach is based on the concept of state. It is the consequence of an important characteristic of a dynamical system, namely that its instantaneous behavior is dependent on its past history, so that the behavior of the system at time t > t 0 can be determined given (1) the forcing function (that is, the input), and (2) the state of the system at t = t 0 In contrast to the trial and error approach of classical control, it is often possible in modern control theory to make use of optimal methods.

dε + ε = R qi dt

the change in the inflow, depends on the characteristics of the regulator valve and may be of the form q i = K ε , K is a constant

Another type of control system is known as an open loop control system, in which the output of the system is not involved in its control. Basically the control on the plant is exercised by a controller as shown in Fig. 1.9. input

output controller

plant

Fig. 1.9 1.4 Mathematical Models and Differential Equations

Many dynamic systems are characterized by differential equations. The process involved, that is, the use of physical laws together with various assumptions of linearity, etc., is known as mathematical modeling. Linear differential equation

&y& + 2 y& − 3 y = u → time-invariant (autonomous) system &y& − 2 t y& + y = u → time-varying (non-autonomous) system Nonlinear differential equation

&y& + 2 y& 2 − 2 y = u &y& − y 2 y& + y = u

1.5 The Classical and Modern Control Theory

Classical control theory is based on Laplace transforms and applies to linear autonomous systems with SISO. A function called transfer function relating the input-output relationship of the system is defined. One of the objectives of control

___________________________________________________________________________________________________________ 3 Chapter 1 System Dynamics and Differential Equation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

2. Transfer Functions and Block Diagrams 2.1 Introduction

Consider - Review of Laplace transform - Using Laplace transform to solve a differential equation

then

2.2 Review of Laplace Transforms

F ( s) =

∫

u = f (t ) dv = e − st dt ∞

⇒

du = d [ f (t )]

f (t ) e − st dt

0 ∞

⎡ ⎞⎤ ⎛ 1 = ⎢ f (t ) ⎜ − e − st ⎟⎥ − ⎠⎦ 0 ⎝ s ⎣

Definition: The Laplace transform of f (t ) , a sectionally continuous function of time, denoted by L [ f (t )] , is defined as

∞

(2.1)

1 1 ⎤ ⎡ = − ⎢ f (t ) e − st ⎥ + s ⎦0 s ⎣

where s = σ + iω , σ and ω are real variables. f (t ) is called

f ( 0) 1 ⎡ d ⎤ = + L ⎢ f (t ) ⎥ s s ⎣ dt ⎦

L [ f (t )] =

∫

∞

f (t ) e

− st

dt = F ( s )

0

inverse transform of F (s) , such as L

−1

[ F ( s )] = f (t ) .

Example 2.1 _______________________________________

⎧0 for t < 0 , c is a Consider the step functions f (t ) = ⎨ ⎩c for t ≥ 0 constant as shown in Fig. 2.1. Find the Laplace transform f (t ) c

∞

0

0

=−

dt ∞

⎞⎛ 1 ⎞⎤ f ( x)dx ⎟⎟ ⎜ − e − st ⎟⎥ s 0 ⎠⎦⎥ 0 ⎝ ⎠

∫

t

∫

∫

⎤ 1 − st ⎞ ⎡ t ⎜ − e ⎟ d ⎢ f ( x)dx ⎥ ⎝ s ⎠ ⎣ 0 ⎦

∞⎛

∫

⎤ 1 − st ⎞ d ⎡ t ⎜ − e ⎟ ⎢ f ( x )dx ⎥ dt ⎠ dt ⎣ 0 ⎝ s ⎦

∞⎛

0

∫

1 ∞ f (t ) e − st dt s 0 F (s) = s =

Some properties of Laplace transforms A is a constant and L [ f (t )] = F ( s )

∫

(5) The delay function

(1) L [ A f (t )] = A L [ f (t )] = A F ( s )

L [ f (t − α )] = e −α s F ( s )

(2) L [ f1 (t ) + f 2 (t )] = L [ f1 ( s )] + L [ f 2 ( s )] = F1 ( s ) + F2 ( s ) (3) Differential

(2.2)

(2.6)

Proof Given the Laplace transform F (s ) of a function f (t ) such that f (t ) = 0 for t < 0 . We wish to find the Laplace transform of f (t − α ) where f (t − α ) = 0 for t < α .(Fig. 2.2) f (t )

Proof d (uv) = udv + vdu 0

− st

0

__________________________________________________________________________________________

0 ∞

⎞

−

c (so long as σ > 0 ) s

∞

d [ f (t )] e − st dt dt

0

∞

∞

0

∫ ⎜⎜⎝ ∫ f ( x)dx ⎟⎟⎠ e

c c ⎡1 ⎤ c e − st dt = −c ⎢ e −st ⎥ = − lim e − st + s t →∞ s ⎣s ⎦0

t →∞

⇒

∞

(2.5)

∞⎛ t

⎡⎛ = ⎢⎜⎜ ⎣⎢⎝

t

⎡d ⎤ L ⎢ f (t )⎥ = sF ( s ) − f (0) dt ⎣ ⎦

∫

1 − st ⎞ ⎜ − e ⎟ d [ f (t )] ⎠ ⎝ s

∫

So long as σ (=real part of s ) is positive, lim e − st = 0 , hence F ( s) =

0

⎤ F (s) ⎡ t L ⎢ f ( x )dx ⎥ = s 0 ⎦ ⎣

Fig. 2.1

∫

∞⎛

(4) Integration

Proof ⎡ t ⎤ L ⎢ f ( x)dx ⎥ = 0 ⎣ ⎦ 0

∫

⎤ ⎡d Hence L ⎢ f (t )⎥ = sF ( s) − f (0) dt ⎦ ⎣

∫

F (s) =

v = − 1s e − st

f (t − α )

∞

∫ ∫ ∫ ∫ udv = [uv] − ∫ vdu

d (uv) = udv + vdu 0 ∞

t 0 α t Fig. 2.2 ___________________________________________________________________________________________________________ 4 Chapter 2 Transfer Functions and Block Diagrams

⇒

0

∞ 0

0

0

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Consider the function F ( s) =

∫

∞

(8) Initial value theorem This theorem makes it possible to find f (0) directly from F (s ) , without actually evaluating the inverse Laplace transform. The theorem states that provided lim sF ( s ) exists, then

f ( x) e − sx dx .

0

Let x = t − α , then F ( s) =

∫

∞

f (t − α ) e − s (t −α ) dt

s→∞

f (0) = lim sF ( s )

0

=e

αs

∫

∞

f (t − α ) e

− st

(2.9)

s→∞

dt

0

= eα s L [ f (t − α )]

Since L [ f ' (t )] = sF ( x) − f (0) and lim

∫

∞

s→∞ 0

f ' (t ) e −st dt = 0 ,

hence 0 = lim sF ( s ) − f (0) , and (2.9) follows.

hence L [ f (t − α )] = e −α s F ( s )

t →∞

(9) Final value theorem This theorem makes it possible to obtain f (∞) directly from F (s ) . Provided the indicated limits exist

(6) The unit impulse function Consider the function shown in Fig. 2.3 ⎧1 0≤t ≤c ⎪ f (t ) = ⎨ c ⎪0 t>c ⎩

lim f (t ) = f (∞) = lim sF ( s )

s→∞

(2.10)

s→0

f (t )

Indeed, as above

1 c

L [ f ' (t )] =

∫

∞

f ' (t ) e − st dt = sF ( s ) − f (0)

0

Now lim

0

t c Fig. 2.3 We can express f (t ) in terms of step functions as

f (t ) =

f ' (t ) e − st dt =

Let denote Dy =

Taking Laplace transform

s e −cs =1 c→0 s

=

∫

⇒ [ s 2Y ( s ) − sy (0) − y ' (0)] − 2[ sY ( s ) − y (0)] + Y ( s ) =

⎛t ⎞ f ⎜ ⎟ e − s t dt ⎝α ⎠ f ( x) e

0

= α F (α s )

L [ D 2 y ] − L [2 Dy ] + L [ y ] = L [e 3t ]

(2.7)

(7) Time scale change Suppose that, L [ f (t )] = F ( s ) , find L [ f (t / α )] , α > 0

∞

d2y dy , D2 y = 2 , … dt dt

Solve the equation D 2 y − 2 Dy + y = e 3t subject to the initial conditions y (0) = 1 and y ' (0) = 2 .

= lim

0

= f ( ∞ ) − f ( 0)

Example 2.6 _______________________________________

1 (1 − e −c s ) c→0 c s d (1 − e −c s ) dc = lim d c→0 (c s ) dc

L [δ (t )] = lim

∫

∞ 0

By taking the Laplace transform of a differential equation it is transformed into an algebraic one in the variable s . This equation is rearranged so that all the terms involving the dependent variable are on one side and the solution is obtained by taking the inverse transform.

Using the results of Example 2.1 and (2.6)

∞

0

2.3 Applications to differential equations

u (t ) − u (t − c) c

⎡ ⎛ t ⎞⎤ L ⎢ f ⎜ ⎟⎥ = ⎣ ⎝ α ⎠⎦

∞

∫ f ' (t )dt = f (t )

s→0

that is, a step function beginning at t = 0 minus a step function beginning at t = c . This leads to the definition of the unit impulse function or the Dirac delta function denoted by δ (t ) defined by

c→0

∞

Hence f (∞) − f (0) = lim sF ( s) − f (0) , and (2.10) follows.

1 1 u (t ) − u (t − c) c c

δ (t ) = lim

∫

s→0 0

− sα x

α dx,

x=

t

α (2.8)

⇒

( s 2 − 2 s + 1) Y ( s ) − s =

⇒

Y (s) =

⇒

Y (s) =

1 s−3

1 s −3

s 2 − 3s + 1 ( s − 3)( s − 1) 2

3/ 4 1/ 2 1/ 4 + + s − 1 ( s − 1) 2 s − 3

Taking inverse transform

___________________________________________________________________________________________________________ 5 Chapter 2 Transfer Functions and Block Diagrams

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

L

⇒

−1

[Y ( s)] =

3 L 4

−1 ⎡

1 ⎤ 1 ⎢ s − 1⎥ + 2 L ⎣ ⎦

−1

⎡ 1 ⎤ 1 ⎢ ⎥+ L 2 ⎣⎢ ( s − 1) ⎦⎥ 4

−1 ⎡

1 ⎤ ⎢ s − 3⎥ ⎣ ⎦

3 1 1 y (t ) = e t + te t + e 3t 4 2 4

__________________________________________________________________________________________

2.4 Transfer functions

Hence G(s) =

H (s) R = Q( s) 1 + A R s

Note that G (s) is dependent only on the system characteristics such as A and R . __________________________________________________________________________________________

u (t )

system

y (t )

Example 2.8 _______________________________________

Find the transfer function of the system described by the

Fig. 2.6 In general form the n th order SISO system (Fig. 2.6) may have an associated differential equation a0 D n y + a1 D n−1 y + K + a n y = b0 D m u + b1 D m−1u + K + bm u (2.11) where u (t ) : input y (t ) : output. n ≥ m : proper system n < m : improper system

equation D 2 y + 3Dy + 2 y = u where u (t ) is the input and y (t ) is the output. Taking Laplace transform both sides with zero-initial conditions, we obtain ( s 2 + 3s + 2) Y ( s ) = U ( s )

Hence G ( s ) =

H (s) 1 = Q( s ) ( s + 1)( s + 2)

__________________________________________________________________________________________

2.5 Block Diagrams

If all initial conditions are zero, taking the Laplace transform of (2.11) yields

Block diagrams The equation Y ( s) = G ( s ) U ( s) is represented by the Fig.2.7

(a0 s n + a1s n−1 + K + a n ) Y ( s ) = (b0 s m + b1s m−1 + K + bm )U ( s )

U (s ) Definition The transfer function G (s ) of a linear autonomous system is the ratio of the Laplace transform of the output to the Laplace transform of the input, with the restriction that all initial conditions are zero.

From the above equation, it follows that G( s) =

Y ( s ) b0 s m + b1s m−1 + K + bm = U ( s ) a0 s n + a1s n−1 + K + a n

E (s) = R( s) − C ( s)

R (s ) (2.12)

(2.12) is also written as

C (s ) Fig. 2.8 Identity diagram To show the relation Y ( s ) = Y1 ( s ) = Y2 ( s)

(2.13)

Obtain the transfer function for a fluid tank having a crosssectional area A with one inflow q (t ) , a head h(t ) and with an outflow valve offering a resistance R . dh 1 1 1 q− h1 + h2 , Using equation (1.7): h&1 = 1 = dt A1 A1 R1 A1 R1 the system can be described as follows dh dh 1 A = q − h or A R +h=qR dt R dt

Y2 ( s ) = Y ( s )

Fig. 2.9 Reduction of block diagram

Y (s) = G ( s ) = B ( s ) A( s ) U ( s) U (s )

On taking the Laplace transform with zero initial conditions

Y1 ( s ) = Y ( s )

Y (s )

Example 2.7 _______________________________________

A R s H ( s) + H ( s) = Q(s) R

Y (s )

Fig. 2.7 Sensing diagram Using for representation of addition and/or subtraction of signals.

The transfer function G (s ) depends only on the system parameters and is independent of the input.

Y ( s) = G ( s)U ( s)

G(s)

A(s)

X (s )

B(s)

(2.15)

Y (s )

U (s )

B(s)A(s)

Y (s )

Fig. 2.10

___________________________________________________________________________________________________________ 6 Chapter 2 Transfer Functions and Block Diagrams

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Block in parallel C ( s ) = A( s ) R ( s ) − A( s ) B ( s)C ( s ) C (s) A( s ) = ⇒ G(s) = R( s ) 1 + A( s ) B( s ) R (s ) E (s )

(2.16)

C (s )

A(s)

R (s ) G ( s )=

A( s ) 1 + A( s ) B ( s )

C (s )

B(s)

Fig. 2.11 Example 2.9 _______________________________________

Obtain the system transfer function for the positive feedback system shown in Fig. 2.12 R (s )

E (s )

A(s)

C (s )

B(s)

Fig. 2.12 E ( s) = R( s) + B( s) C ( s) C ( s ) = A( s ) E ( s )

⇒

C ( s ) [1 − A( s ) B ( s )] = A( s ) R ( s )

and G ( s ) =

C (s) A( s ) = R( s ) 1 − A( s ) B ( s )

A special case and an important case of a feedback system is the unity feedback system shown in Fig. 2.1.3 R (s )

E (s )

A(s)

C (s )

Fig. 2.13 The unity feedback system transfer function is G ( s) =

A( s ) 1 + A( s )

(2.17)

__________________________________________________________________________________________

Example 2.10_______________________________________

Consider a system having the transfer function A( s ) =

K . s−2

Time response of this system to a unit impulse Y ( s ) = A( s ) U ( s ), U ( s ) = 1

⇒

y (t ) = K e 2t

For a unity feedback closed loop system G(s) =

⇒

A( s ) K = 1 + A( s) s + ( K − 2)

y (t ) = K e −( K − 2) t

__________________________________________________________________________________________

___________________________________________________________________________________________________________ 7 Chapter 2 Transfer Functions and Block Diagrams

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

3. State Space Formulation 3.1 Introduction (Review) 3.1.1 Eigenvalues and eigenvectors Consider a matrix A of order n × n . If there exists a vector x ≠ 0 and a scalar λ such that Ax = λ x

⎡λ1 0 ⎢ Λ = ⎢ 0 λ2 M M ⎢0 0 ⎣

L L O L

⎤ ⎥ ⎥ = diag{λ1 , λ2 ,L, λn } λn ⎥⎦ 0 0 M

is called eigenvalue matrix of A . Hence

then x is called an eigenvector of A . λ is called an eigenvalue of A . The above equation can be written in the form [ λ I − A]x = 0

Λ = X −1 A X

3.1.2 The Cayley-Hamilton theorem Given a square matrix A (of order n × n ) and integers r and s , then Ar A s = Ar + s = A s+r = A s Ar

where I is the unit matrix (of order n × n ). It is known that the above homogeneous equation has non-trivial (that is, nonzero) solution only if the matrix [ λ I − A ] is singular, that is if det(λ I − A) = 0

This is an equation in λ of great importance. We denoted it by c(λ ) , so that c(λ ) = det(λ I − A) = 0

It is called the characteristic equation of A . Written out in full, this equation has the form a12 ⎡λ − a11 ⎢ a 21 λ − a 22 c (λ ) = ⎢ M ⎢ M an 2 ⎣ a n1

L a1n ⎤ L a 2n ⎥ = 0 O M ⎥⎥ L λ − a nn ⎦

On expansion of determinant, c(λ ) is seen to be a polynomial of degree n in λ , having the form c(λ ) = λn + b1λn−1 + K + bn−1λ + bn = (λ − λ1 )(λ − λ2 ) L (λ − λn ) =0

λ1 , λ2 ,L, λn , the roots of c(λ ) = 0 , are the eigenvalues of A . Assuming that A has distinct eigenvalue λ1 , λ2 ,L, λn , the corresponding eigenvectors x1 , x 2 ,L, x n are linearly independent. The (partitioned) matrix X = [x1 x 2 L x n ]

is called the modal matrix of A. Since A xi = λ xi

where

(i = 1,2,L, n)

This property leads to the following definition. Definition Corresponding to a polynomial in a scalar variable λ f (λ ) = λk + b1λk −1 + K + bk −1λ + bk

define the (square) matrix f ( A) , called a matrix polynomial, by f ( A) = A k + b1 A k −1 + K + bk −1 A + bk I

where I is the unit matrix of order n × n . For example, corresponding to

f (λ ) = λ2 − 2λ + 3 and

A = ⎡ 1 1⎤ , we have ⎢⎣− 1 3⎥⎦ f ( A) = ⎡⎢ 1 1⎤⎥ ⎡⎢ 1 1⎤⎥ − 2 ⎡⎢ 1 1⎤⎥ + 3⎡⎢1 0⎤⎥ = ⎡⎢ 1 2⎤⎥ ⎣− 1 3⎦ ⎣− 1 3⎦ ⎣− 1 3⎦ ⎣0 1⎦ ⎣− 2 5⎦

Of particularly interest to us are polynomials f having the property that f ( A) = 0 . For every matrix A one such polynomial can be found by the Cayley-Hamilton theorem which states that: Every square matrix A satisfies its own characteristic equations. For the above example, the characteristic equation of A is

c(λ ) = (λ − λ1 )(λ − λ2 ) = λ2 − 4λ + 4 where λ1 , λ2 are eigenvalues of A . So that c( A) = A 2 − 4 A + 4 I = ⎡ 1 1⎤ ⎡ 1 1⎤ − 4 ⎡ 1 1⎤ + 4 ⎡1 0⎤ ⎢⎣− 1 3⎥⎦ ⎢⎣− 1 3⎥⎦ ⎣⎢− 1 3⎥⎦ ⎣⎢0 1⎥⎦ = ⎡⎢0 0⎤⎥ ⎣0 0 ⎦

it follows that AX = ΛX

In fact the Cayley-Hamilton theorem guarantees the existence of a polynomial c(λ ) of degree n such that c( A) = 0 .

___________________________________________________________________________________________________________ 8 Chapter 3 State Space Formulation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

3.2 State Space Forms

⎡ x1 ⎤ y = [1 0 0] ⎢ x 2 ⎥ ⎢x ⎥ ⎣ 3⎦

Consider the system equation in the form y ( n ) + a1 y ( n −1) + K + a n −1 y& + a n y = u

(3.1)

It is assumed that y (0), y& (0), L , y ( n−1) (0) are known. If we

define x1 = y , x 2 = y& , L , x n = y ( n −1) , (3.1) is written as x&1 = x 2 x& 2 = x 3

Let x1 =

1 y + 1 &y& 5 5 1 ( −1 + 2i ) &y& x 2 = 15 (2 + i ) y − 12 iy& + 10 1 (1 + 2i ) &y& x 3 = 15 (2 − i ) y + 12 iy& − 10

(3.6)

where i 2 = −1 . Then

M

x& n −1 = x n x& n = − a n x1 − a n−1 x 2 − K − − a1 x n + u

which can be written as a vector-matrix differential equation ⎡ x&1 ⎤ ⎡ 0 1 ⎢ x& 2 ⎥ ⎢ 0 0 ⎢ M ⎥=⎢ M M ⎢ x& ⎥ ⎢ 0 0 − 1 n ⎢ ⎥ ⎢ ⎢⎣ x& n ⎥⎦ ⎣− a n − a n−1

(2) Diagonal form

L 0 0 ⎤ ⎡ x1 ⎤ ⎡0⎤ L 0 0 ⎥ ⎢ x 2 ⎥ ⎢0 ⎥ O M M ⎥ ⎢ M ⎥ + ⎢M⎥ u L 0 1 ⎥ ⎢ x n −1 ⎥ ⎢0⎥ L − a 2 − a1 ⎥⎦ ⎢⎢ x ⎥⎥ ⎢⎣1 ⎥⎦ ⎣ n ⎦

(3.2)

(3.7)

__________________________________________________________________________________________

In general form, a MIMO system has state equations of the form

that is, as x& = A x + B u , where x , A and B are defined in equation (3.2). The output of the system

⎡ x1 ⎤ ⎢ ⎥ y = [1 0 L 0] ⎢ x 2 ⎥ M ⎢x ⎥ ⎣ n⎦

⎡ x&1 ⎤ ⎡2 0 0 ⎤ ⎡ x1 ⎤ 1 ⎡ 2 ⎤ ⎢ x& 2 ⎥ = ⎢0 i 0 ⎥ ⎢ x 2 ⎥ + ⎢− 1 + 2i ⎥ u ⎢ x& ⎥ ⎢0 0 − i ⎥ ⎢ x ⎥ 10 ⎢ − 1 − 2i ⎥ ⎦⎣ 3⎦ ⎣ ⎦ ⎣ 3⎦ ⎣ and ⎡ x1 ⎤ y = [1 1 1] ⎢ x 2 ⎥ ⎢x ⎥ ⎣ 3⎦

(3.3)

⎧ x& = A x + B u ⎨ ⎩ y = C x + Du

(3.14)

and a SISO system has state equations of the form as (3.4). 3.3 Using the transfer function to define state variables

that is, as, where C = [1 0 L 0] . The combination of equations (3.2) and (3.3) in the form

It is sometimes possible to define suitable state variables by considering the partial fraction expansion of the transfer function. For example, given the system differential equation

⎧ x& = A x + B u ⎨ ⎩y = C x

&y& + 3 y& + 2 y = u& + 3u

(3.4)

The corresponding transfer function is are known as the state equations of the system considered. The matrix A in (3.2) is said to be in companion form. The components of x are called the state variables x1 , x 2 , L , x n . The corresponding n-dimensional space called the state space. Any state of the system is represented by a point in the state space.

G(s) = Hence

Y (s) = X 1 (s) + X 2 (s)

Example 3.1 _______________________________________

2U ( s ) s +1 U (s) X 2 (s) = − s+2 X 1 (s) =

Obtain two forms of state equations of the system defined by

&y&& − 2 &y& + y& − 2 y = u where matrix A corresponding to one of the forms should be diagonal.

On taking the inverse Laplace transforms, we obtain x&1 = − x1 + 2u x& 2 = −2 x 2 − u

(a) A in companion form Let the state variable as x1 = y , x 2 = y& , x 3 = &y& . Then

⎡ x&1 ⎤ ⎡0 1 0⎤ ⎡ x1 ⎤ ⎡0⎤ ⎢ x& 2 ⎥ = ⎢0 0 1 ⎥ ⎢ x 2 ⎥ + ⎢0⎥ u ⎢ x& ⎥ ⎢2 − 1 2⎥ ⎢ x ⎥ ⎢1 ⎥ ⎦⎣ 3⎦ ⎣ ⎦ ⎣ 3⎦ ⎣ and

Y ( s) 2 1 s+3 = = − U ( s ) ( s + 1)( s + 2) s + 1 s + 2

or

(3.5)

⎡ x&1 ⎤ ⎡− 1 0 ⎤ ⎡ x1 ⎤ ⎡ 2 ⎤ ⎢⎣ x& 2 ⎥⎦ = ⎢⎣ 0 − 2⎦⎥ ⎢⎣ x 2 ⎥⎦ + ⎢⎣− 1⎥⎦ u ⎡x ⎤ y = [1 1] ⎢ 1 ⎥ ⎣x2 ⎦

___________________________________________________________________________________________________________ 9 Chapter 3 State Space Formulation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

We can of course make a different choice of state variables, for example

x = 2e −3t + 4e −3t

=

U (s) s +1 U (s) X 2 (s) = s+2 then X 1 (s) =

3τ

0

= 2e −3t + 4e −3t

Y (s) = 2 X 1 (s) − X 2 (s)

t

∫ e τ dτ ( te − e

22 e −3t 9

3t

1 3

+ 19

1 3t 9

)

+ 43 t − 94

__________________________________________________________________________________________

To use this method to solve the vector matrix equation (3.15), we must first define a matrix function which is the analogue of the exponential, and we must also define the derivative and the integral of a matrix.

x&1 = − x1 + u &x 2 = −2 x 2 + u

∞

Definition Since e z =

∑ n! z 1

n

(all z ), we define

0

∞

Now the state equations are eA =

∑ n! A 1

n

= A0 + A +

1 2 1 A + K = I + A + A2 + K 2! 2!

⎡ x&1 ⎤ ⎡− 1 0 ⎤ ⎡ x1 ⎤ ⎡1⎤ ⎢⎣ x& 2 ⎥⎦ = ⎢⎣ 0 − 2⎦⎥ ⎢⎣ x 2 ⎥⎦ + ⎢⎣1⎥⎦ u

(where A0 ≡ I ) for every square matrix A .

⎡x ⎤ y = [2 − 1] ⎢ 1 ⎥ ⎣ x2 ⎦

Example 3.4 _______________________________________

0

⎡ 1 0 0⎤ Evaluate e At given A = ⎢− 1 0 0⎥ ⎢⎣ 1 1 1 ⎥⎦

3.4 Direct solution of the state equation By a solution to the state equation x& = A x + B u

(3.15)

We must calculate A 2 , A3 ,… ⎡ 1 0 0⎤ A 2 = ⎢ − 1 0 0 ⎥ = A = A3 = K ⎢⎣ 1 1 1⎥⎦

we mean finding x at any time t given u (t ) and the value of x at initial time t 0 = 0 , that is, given x(t 0 ) = x 0 . It is instructive to consider first the solution of the corresponding scalar differential equation

It follows that

x& = a x + b u

e At = I + A t +

(3.16)

1 2 2 1 33 A t + A t +K 2! 3 1 2 1 3 ⎛ ⎞ = I + A ⎜ t + t + t + K⎟ 2 ! 3 ⎝ ⎠

given x = x0 at t = 0 . The solution of (3.15) is found by an analogous method. Multiplying the equation by the integrating factor e − at , we obtain

= I + A(e t − 1) t 0 0 ⎤ ⎡1 0 0⎤ ⎡⎢ e − 1 ⎥ 0 ⎥ = ⎢0 1 0⎥ + ⎢ − e t + 1 0 t t t ⎣⎢0 0 1⎥⎦ ⎢⎣ e − 1 e − 1 e − 1⎦⎥ ⎡ et 0 0⎤ ⎥ ⎢ = ⎢− e t + 1 1 0⎥ t t t e e e − − 1 1 ⎦⎥ ⎣⎢

e − at ( x& − ax ) = e − at b u

or

d −at [e x] = e −at b u dt

Integrating the result between 0 and t gives

__________________________________________________________________________________________

e

−at

x − x0 =

t

∫e

− aτ

b u (τ ) dτ

0

that is x=

e at x0 123

complementary function

+

t

∫ e b u(τ ) dτ 144424443 −a( t −τ )

(3.17)

0

particular int ergral

Example 3.3 _______________________________________ Solve the equation x& + 3 x = 4t given that x = 2 when t = 0 .

In fact, the evaluation of the matrix exponential is not quite as difficult as it may appear at first sight. We can make use of the Cayley-Hamilton theorem which assures us that every square matrix satisfies its own characteristic equation. One direct and useful method is the following. We know that if A has distinct eigenvalues, say λ1 , λ2 ,L, λn , there exists a non-singular matrix P such that ⎡λ1 0 ⎢ P A P = ⎢ 0 λ2 M M ⎢0 0 ⎣ −1

L L O L

0⎤ 0 ⎥ = diag{λ , λ ,L, λ } = Λ 1 2 n M ⎥ λn ⎥⎦

We then have A = P Λ P −1 , so that Since a = −3 , b = 4 and u = t , on substituting into equation (3.17) we obtain A 2 = ( PΛP −1 )( PΛP −1 ) = PΛ ( P −1P )ΛP −1 = PΛ2 P −1 ___________________________________________________________________________________________________________ 10 Chapter 3 State Space Formulation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

e At = ⎡⎢ 2 2 ⎤⎥ e −t + ⎡⎢−1 −2⎤⎥ e −2t ⎣1 1⎦ ⎣− 1 − 1⎦

and in general A r = PΛr P −1

(r = 1,2,L)

If we consider a matrix polynomial, say f ( A) = A − 2 A + I , we can write it as 2

= P f (Λ ) P −1

Since (for the above example)

f ( λ ) = λ 2 − 2λ + I ⎤ ⎡λ1 0 ⎥ ⎢ ⎥ − 2 ⎢ 0 λ2 M M ⎥ ⎢ λ2n ⎥⎦ ⎣ 0 0 0 0 M

⎡λ12 − 2 0 ⎢ λ22 − 2λ2 + 1 =⎢ 0 ⎢ M M ⎢ 0 0 ⎣

called

the

poles

and

Definition Let A(t ) = [aij (t )], then

= P diag{ f (λ1 ), f (λ2 ),L, f (λn )}P −1

L L O L

λ1 , λ2 ,L, λn are

e λ1t , e λ2t ,L, e λnt are called the modes of the system. We now define the derivative and the integral of a matrix A(t ) whose elements are functions of t .

= P ( Λ2 − 2Λ + I ) P −1

0 ⎢ 2 0 λ 2 =⎢ ⎢M M ⎢0 0 ⎣

eigenvalues

__________________________________________________________________________________________

f ( A) = PΛ2 P −1 − 2 PΛP −1 + P I P −1

⎡λ12

In the solution to the unforced system equation x = e At x 0 , the

L L O L

0 ⎤ ⎡1 0 0 ⎥ + ⎢0 1 M ⎥ ⎢M M λn ⎥⎦ ⎢⎣0 0

L L O L

0⎤ 0⎥ M⎥ 1⎥⎦

⎤ L 0 ⎥ L 0 ⎥ ⎥ O M L λ2n − 2λn + 1⎥⎦

(1)

d d A(t ) = A& (t ) = [ (aij )], and dt dt

(2)

∫ A(t )dt = [∫ a

that is, each element of the matrix is differentiated (or integrated).

⎡ 6t sin 2t ⎤ For example, if A = ⎢ 2 , then 3 ⎥⎦ ⎣t + 2 A& = ⎡⎢ 6 2 cos 2t ⎤⎥ , and 0 ⎦ ⎣2t

∫

= diag{ f (λ1 ), f (λ2 ),L, f (λn )} The results holds for the general case when f (x) is a polynomial of degree n . In generalizing the above, taking f ( A) = e At , we obtain e At = P diag{e λ1t , e λ2t ,L, e λnt } P −1

ij (t ) dt ],

(3.18)

⎡ 3t 2 − 12 cos 2t ⎤ ⎥ + C ( C is a constant matrix) Adt = ⎢ 3 3t ⎥⎦ ⎢⎣ 13 t + 2t

From the definition a number of rules follows: if α and β are constants, and A and B are matrices, d (i) (α A + β B ) = α A& + β B& dt (ii)

b

b

b

a

a

a

∫ (α A + β B)dt = α ∫ Adt + β ∫ Bdt

Example 3.5 _______________________________________

d (iii) ( A B ) = A B& + A& B dt

Find e At given A = ⎡ 0 2 ⎤ ⎢⎣− 1 − 3⎥⎦

(iv) If A is a constant matrix, then

The eigenvalues of A are λ1 = −1 and λ2 = −2 . It can be verified that P = ⎡ 2 1 ⎤ and that P −1 = ⎡ 1 1 ⎤ . ( P can ⎢⎣− 1 − 1⎥⎦ ⎢⎣− 1 − 2⎥⎦

⊗ Note that although the above rules are analogous o the rules for scalar functions, we must not be dulled into accepting for matrix functions all the rules valid for scalar functions. For example, although d ( x n ) = n x n−1 x& in general, dt

be found from A = P Λ P −1 ). Using (3.18), we obtain

it is not true that

⎡ −t 0 ⎤ ⎡ 1 1 ⎤ e At = ⎡ 2 1 ⎤ ⎢e ⎢⎣− 1 − 1⎥⎦ 0 e −2t ⎥ ⎢⎣− 1 − 2⎥⎦ ⎣ ⎦ ⎡ 2e −t − e −2t 2e −t − 2e −2t ⎤ = ⎢ −t ⎥ − 2t − e −t + 2e −2t ⎦ ⎣− e + e

then

e At in its spectral or modal form as

and so that

(3.19)

where λ1 , λ2 ,L, λn are the eigenvalues of A and A1 , A2 ,L , An are matrices of the same order as the matrix A . In the above example, we can write

d dt

( A n ) = n A n−1 A& .For example,

⎡ 2 ⎤ A = ⎢t 2t ⎥ , ⎣0 3⎦ 2 d ⎡ 3 ⎤ ( A 2 ) = ⎢4t 6t + 6⎥ 0 0 dt ⎣ ⎦ ⎡4t 3 4t 2 ⎤ & 2 AA = ⎢ 0 ⎥⎦ ⎣ 0

when

From the above discussion it is clear that we can also write

e At = A1e λ1t + A2 e λ2t + K + An e λnt

d At (e ) = Ae At dt

d 2 A ≠ 2 AA& . dt

We now return to our original problem, to solve equation (3.15): x& = A x + B u . Rewrite the equation in the form x& − A x = B u

⇒

e − At (x& − A x) = e − At ( B u )

___________________________________________________________________________________________________________ 11 Chapter 3 State Space Formulation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

or

With this notation equation (3.20) becomes

d − At (e x) = e − At B u dt

∫

e − At x(t ) − x(0) =

t

∫e

− Aτ

3.5 Solution of the state equation by Laplace transforms

B u (τ ) dτ

Since the state equation is in a vector form we must first define the Laplace transform of a vector.

so that t

∫e

A(t −τ )

B u (τ ) dτ

(3.20) Definition

0

Example 3.6 _______________________________________

A system is characterized by the state equation ⎡ x&1 ⎤ ⎡ 0 2 ⎤ ⎡ x1 ⎤ ⎡0⎤ ⎢⎣ x& 2 ⎥⎦ = ⎢⎣− 1 − 3⎥⎦ ⎢⎣ x2 ⎥⎦ + ⎢⎣1⎥⎦ u If the forcing function u (t ) = 1 for t ≥ 0 and x(0) = [1 − 1]T . Find the state x of the system at time t . We have already evaluated e At in Example 3.5. It follows that −2t −t −2t ⎤ ⎡ −t e At x(0) = ⎢ 2e −t − e −2t 2e −t − 2e −2t ⎥ ⎡ 1 ⎤ − e + 2e ⎦ ⎢⎣− 1⎥⎦ ⎣− e + e

= e −2t ⎡⎢ 1 ⎤⎥ ⎣− 1⎦

∫

t

∫ ⎡ = ⎢2e ∫ ⎣ 2e t

0

−(t −τ )

⎡ L [ x1 (t )] ⎤ ⎡ X 1 ( s ) ⎤ ⎢ L [ x2 (t )]⎥ ⎢ X 2 ( s) ⎥ We define L [x(t )] = ⎢ ⎥ = ⎢ M ⎥ = X( s ) M ⎢ L [ x (t )]⎥ ⎢ X ( s)⎥ n ⎣ ⎦ ⎣ n ⎦ From this definition, we can find all the results we need. For example, ⎡ L [ x&1 (t )] ⎤ ⎡ sX 1 ( s) − x1 (0) ⎤ ⎥ ⎢ ⎥ ⎢ & L [x& (t )] = ⎢ L [ x 2 (t )]⎥ = ⎢ sX 2 ( s) − x2 (0) ⎥ M M ⎢ L [ x& (t )]⎥ ⎢ sX ( s) − x (0)⎥ n n ⎦ ⎣ ⎦ ⎣ n

s X( s ) − x(0) = A X( s ) + B U ( s ) where U ( s ) = L [u (t )]

− 2e −2(t −τ ) ⎤ dτ ⎥ − e −(t −τ ) ⎦

−2(t −τ )

or

−t − 2t ⎤ ⎡ = ⎢1 − 2−et +−e2t ⎥ e e − ⎣ ⎦

−t −2t ⎤ ⎡ −t −2t ⎤ ⎡ −2t ⎤ ⎡ x(t ) = ⎢ e −2t ⎥ + ⎢1 − 2−et +−e2t ⎥ = ⎢1 − 2−et + 2−e2t ⎥ ⎣− e ⎦ ⎣ e − e ⎦ ⎣ e − 2e ⎦

X( s) = ( sI − A) −1 x(0) + ( sI − A) −1 B U ( s)

__________________________________________________________________________________________

At

The matrix e in the solution equation (3.20) is of special interest to control engineers; they call it the state-transition matrix and denote it by Φ(t ) , that is, Φ(t ) = e At

(3.21)

For the unforced system (such as u (t ) = 0 ) the solution (Eq. (3.20)) becomes x(t ) = Φ(t ) x(0)

so that Φ(t ) transforms the system from its state x(0) at some initial time t = 0 to the state x(t ) at some subsequent time t hence the name given to the matrix. Since e At e − At = I . It follows that [e At ]−1 = e − At .

Φ −1 (t ) = e − At = Φ(−t ) At − Aτ

Φ(t ) Φ(−τ ) = e e

=e

A(t −τ )

( sI − A) X( s ) = x(0) + B U ( s )

Unless s happens to be an eigenvalue of A , the matrix ( sI − A) is non-singular, so that the above equation can be solved giving

Hence the state of the system at time t is

Also

⎡ x1 ⎤ ⎢x ⎥ Let x(t ) = ⎢ 2 ⎥ M ⎢x ⎥ ⎣ n⎦

Now we can solve the state equation (3.15). Taking the transform of x& = A x + B u , we obtain

t

e A(t −τ ) B u (τ ) dτ = e A(t −τ ) ⎡⎢0⎤⎥ dτ ⎣1⎦ 0 0

Hence

(3.22)

0

0

x(t ) = e At x(0) =

t

x(t ) = Φ(t ) x(0) = Φ(t − τ ) B u (τ ) dτ

On integration, this becomes

= Φ (t − τ ) .

(3.23)

and the solution x(t ) is found by taking the inverse Laplace transform of equation (3.23).

Definition ( sI − A) −1 is called the resolvent matrix of the system. On comparing equation (3.22) and (3.23) we find that Φ (t ) = L

−1

{( sI − A) −1}

Not only the use of transforms a relatively simple method for evaluating the transition matrix, but indeed it allows us to calculate the state x(t ) without having to evaluate integrals.

Example 3.7 _______________________________________ Use Laplace transform to evaluate the state x(t ) of the system describe in Example 3.6. ⎡ x&1 ⎤ ⎡ 0 2 ⎤ ⎡ x1 ⎤ ⎡0⎤ ⎢⎣ x& 2 ⎥⎦ = ⎢⎣− 1 − 3⎥⎦ ⎢⎣ x2 ⎥⎦ + ⎢⎣1⎥⎦ u

___________________________________________________________________________________________________________ 12 Chapter 3 State Space Formulation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

T z& = A T z + B u y = CT z

For this system ( sI − A) = ⎡⎢ s −2 ⎤⎥ , so that ⎣1 s + 3⎦

( sI − A) −1 =

or as

1 ⎡ s + 3 2⎤ s ( s + 3) + 2 ⎢⎣ − 1 s ⎥⎦

z& = A1 z + B1 u

(3.26)

y = C1 z

s+3 2 ⎡ ⎤ ⎢ ( s + 1)( s + 2) ( s + 1)( s + 2) ⎥ ⎥ =⎢ −1 s ⎢ ⎥ ⎢⎣ ( s + 1)( s + 2) ( s + 1)( s + 2) ⎥⎦

where A1 = T −1 A T , B1 = T −1 B, C1 = C T . The transformation (3.25) is called state-transformation, and the matrices A and A1 are similar. We are particular interested in the

1 2 2 ⎤ ⎡ 2 − ⎢ s +1 − s + 2 + + s 1 s 2 ⎥ =⎢ 1 1 1 2 ⎥ ⎢− ⎥ + − + s +1 s + 2⎦ ⎣ s +1 s + 2

transformation when A1 is diagonal (usually denoted by Λ ) and A is in the companion form, such as (3.2)

So that −2t −t −2t ⎤ ⎡ −t L −1{( SI − A) −1} = ⎢ 2e −t − e −2t 2e −t − 2e −2t ⎥ = Φ(t ) − e + e − e + 2 e ⎣ ⎦

Hence the complementary function is as in Example 3.6. For the particular integral, we note that since

⎡ x&1 ⎤ ⎡ 0 1 ⎢ x& 2 ⎥ ⎢ 0 0 ⎢ M ⎥=⎢ M M ⎢ x& ⎥ ⎢ 0 0 1 n − ⎥ ⎢ ⎢ ⎢⎣ x& n ⎥⎦ ⎣− a n − a n −1

It is assumed that the matrix A has distinct eigenvalues λ1 , λ2 ,L, λn . Corresponding to the eigenvalue λi there is the

1 L {u (t )} = , then s

eigenvector x i such that

1 1 ⎡ s + 3 2⎤ ⎡0⎤ ( sI − A) −1 BU ( s ) = ( s + 1)( s + 2) s ⎢⎣ − 1 s ⎥⎦ ⎣⎢1 ⎥⎦

A x i = λi xi

2 ⎡ ⎤ ⎢ s ( s + 1)( s + 2) ⎥ ⎥ =⎢ 1 ⎢ ⎥ ⎣⎢ ( s + 1)( s + 2) ⎦⎥

V is called the modal matrix, it is non-singular and can be used as the transformation matrix T above. We can write the n equations defined by equation (3.27) as AV = Λ V

−t −2 t ⎤ ⎡ {( sI − A) BU ( s )} = ⎢1 − 2−et +−e2t ⎥ − e e ⎣ ⎦ −1

which is the particular integral part of the solution obtained in Example 3.6. __________________________________________________________________________________________

3.6 The transformation from the companion to the diagonal state form Note that the choice of the state vector is not unique. We now assume that with one choice of the state vector the state equations are

x& = A x + B u y =Cx

(3.24)

where ⎡λ1 0 ⎢ Λ = ⎢ 0 λ2 M M ⎢0 0 ⎣

(3.28)

L L O L

0⎤ 0 ⎥ = diag{λ , λ ,L, λ } 1 2 n M ⎥ λn ⎥⎦

From equation (3.28) we obtain Λ = V −1 AV

(3.29)

The matrix A has the companion form 1 0 ⎡ 0 ⎢ 0 0 1 M M A=⎢ M ⎢ 0 0 0 ⎢⎣− a0 − a1 − a 2

L 0 ⎤ L 0 ⎥ O M ⎥ L 1 ⎥ L − a n−1 ⎥⎦

where A is the matrix of order n × n and B and C are matrices of appropriate order.

so that its characteristic equation is

Consider any non-singular matrix T of order n × n . Let

det( sI − A) = λn + a n−1λn−1 + K + a1λ + a0 = 0 ]

x =T z

(3.27)

V = [x1 x 2 L x n ]

On taking the inverse Laplace transform, we obtain L

(i = 1,2,L, n)

We define the matrix V whose columns are the eigenvectors

2 1 ⎤ ⎡1 + ⎢ − ⎥ = ⎢ s s +1 s + 2 ⎥ 1 1 ⎢ ⎥ − ⎣ s +1 s + 2 ⎦

−1

0 ⎤ ⎡ x1 ⎤ ⎡0⎤ L 0 0 ⎥ ⎢ x 2 ⎥ ⎢0 ⎥ L 0 O M M ⎥ ⎢ M ⎥ + ⎢M⎥ u 1 ⎥ ⎢ x n −1 ⎥ ⎢0⎥ L 0 L − a 2 − a1 ⎥⎦ ⎢⎢ x ⎥⎥ ⎢⎣1 ⎥⎦ ⎣ n ⎦

(3.25)

By solving this equation we obtain the eigenvalues λ1 , λ2 ,L, λn . The corresponding eigenvectors have an

then z is also a state vector and equation (3.24) can be written interesting form. Consider one of the eigenvalues λ and the as ___________________________________________________________________________________________________________ 13 Chapter 3 State Space Formulation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

corresponding eigenvector x = [α1 α 2 L α n ]T . Then the equation A x = λ x , corresponds to the system of equations

⎧α 2 ⎪ ⎪α 3 ⎨ ⎪ ⎪α n ⎩

so the transformation V −1x gives the new choice of the state vector ⎡ z1 ⎤ 0 4 ⎤ ⎡ y⎤ ⎡ 4 ⎢ z 2 ⎥ = 1 ⎢8 + 4i − 10i − 2 + 4i ⎥ ⎢ y& ⎥ ⎢ z ⎥ 20 ⎢8 − 4i 10i − 2 − 4i ⎥ ⎢ &y&⎥ ⎣ ⎦⎣ ⎦ ⎣ 3⎦

= λ α1 = λα2 M = λ α n−1

[

Setting α1 = 1 , we obtain x = 1 λ λ2 L λ the modal matrix in this case takes the form 1 ⎡ 1 ⎢ λ1 λ1 V =⎢ M M ⎢ n−1 n−1 λ1 ⎣λ1

]

n −1 T

L 1 ⎤ L λ1 ⎥ O M ⎥ ⎥ L λ1n−1 ⎦

The state equations are now in the form of equation (3.26), that is . Hence

z& = A1z + B1u y = C1z where A1 = V −1 AV = diag{2, i,−i}

(3.30)

1 ⎡ 2 ⎤ ⎢− 1 + 2i ⎥ 10 ⎢ − 1 − 2i ⎥ ⎣ ⎦ C1 = C V = [1 1 1] B1 = V −1B =

In this form V is called Vandermonde matrix; it is nonsingular since it has benn assumed that λ1 , λ2 ,L, λn are distinct. We now consider the problem of obtaining the transformation from the companion form to diagonal form.

__________________________________________________________________________________________

Example 3.8 _______________________________________

Consider the system in state space form of Eq.(3.24)

Having chosen the state vector so that

x& = A x + B u

&y&& − 2 &y& _ y& − 2 y = u

y =Cx

is written as the equation (in companion form)

Taking the Laplace transforms we obtain

⎡ x&1 ⎤ ⎡0 1 0⎤ ⎡ x1 ⎤ ⎡0⎤ ⎢ x& 2 ⎥ = ⎢0 0 1 ⎥ ⎢ x 2 ⎥ + ⎢0⎥ u ⎢ x& ⎥ ⎢2 − 1 2⎥ ⎢ x ⎥ ⎢1⎥ ⎦⎣ 3⎦ ⎣ ⎦ ⎣ 3⎦ ⎣ ⎡ x1 ⎤ y = [1 0 0]⎢ x 2 ⎥ ⎢x ⎥ ⎣ 3⎦

s X (s) = A X (s) + B U (s)

Find the transformation which will transform this into a state equation with A in diagonal form. The characteristic equation of A is det(λ I − A) = λ3 − 2λ2 + λ − 2 = (λ − 2)(λ2 + 1) = 0

that is λ1 = 2, λ2 = i and λ3 = −i . From (3.30) the modal matrix is ⎡1 1 1 ⎤ V = ⎢2 i − i ⎥ ⎢⎣4 − 1 − 1⎥⎦ and its inverse can be shown to be V −1 =

0 4 ⎤ 1 ⎡ 4 ⎢8 + 4i − 10i − 2 + 4i ⎥ 20 ⎢8 − 4i 10i − 2 − 4i ⎥ ⎣ ⎦

The transformation is defined by equation (3.25), that is x = V z or z = V −1x . The original choice for x was

x = [ y y& &y&]

T

3.7 The transfer function from the state equation

(3.24)

Y (s) = C X (s) ⇒

X ( s ) = ( sI − A) −1 B U ( s )

and G ( s ) =

Y (s) = C ( sI − A) −1 B U (s)

(3.31)

Example 3.9 _______________________________________ Calculate the transfer function from the system whose state equations are x& = ⎡1 −2 ⎤ x + ⎡1 ⎤ u ⎢⎣3 − 4⎥⎦ ⎢⎣2⎥⎦ y = [1 − 2]x

G ( s ) = [1 − 2]⎡⎢ s − 1 2 ⎤⎥ ⎣ − 3 s + 4⎦

−1

⎡1 ⎤ ⎢⎣2⎥⎦ −1

s+4 −2 ⎤ ⎡ ⎢ ( s + 1)( s + 2) ( s + 1)( s + 2) ⎥ 1 ⎥ ⎡ ⎤ = [1 − 2]⎢ s−2 3 ⎥ ⎢⎣2⎥⎦ ⎢ ⎢⎣ ( s + 1)( s + 2) ( s + 1)( s + 2) ⎥⎦ =−

3s + 2 ( s + 1)( s + 2)

__________________________________________________________________________________________

Example 3.10_______________________________________ Given that the system

x& = A x + B u y =Cx

___________________________________________________________________________________________________________ 14 Chapter 3 State Space Formulation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

has the transfer function G (s) , find the transfer function of the system

z& = T −1 AT z + T −1 B u y = CT z If the transfer function is G1 ( s ) G1 ( s ) = C T ( sI − T −1 AT ) −1T −1 B = C T {T −1 ( sI − A)T }−1T −1 B = C ( sI − A) −1 B = G(s) so that G1 ( s ) = G ( s ) . __________________________________________________________________________________________

___________________________________________________________________________________________________________ 15 Chapter 3 State Space Formulation

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.4 Transient and Steady State Response Analysis 4.1 Introduction

Notice that when t = a , then y (t ) = y (a) = 1 − e −1 = 0.63 . The

Many applications of control theory are to servomechanisms which are systems using the feedback principle designed so that the output will follow the input. Hence there is a need for studying the time response of the system.

response is in two-parts, the transient part e −t / a , which approaches zero as t → ∞ and the steady-state part 1, which is the output when t → ∞ .

The time response of a system may be considered in two parts: • Transient response: this part reduces to zero as t → ∞ • Steady-state response: response of the system as t → ∞

Consider the output of a linear system in the form (4.1)

K1 K2 z− p z , where K1 = K and K 2 = K − s s+ p p p

Y ( s) = Hence,

K1 {

steady − state part

For a transient response analysis it is customary to use a reference unit step function u (t ) for which 1 s

− K 2 e − pt 1 424 3

(4.5)

transient part

With the assumption that z > p > 0 , this response is shown in Fig. 4.2. K1

y = K1 − K 2e − pt

K2

It then follows that

1 1 1 = − (a s + 1) s s s + 1 / a

Y ( s) =

where K =b/a z = 1 / b : the zero of the system p = 1 / a : the pole of the system

y (t ) =

1 U (s) as + 1

U (s) =

(4.4)

When U ( s ) = 1 / s , Eq. (4.4) can be written as

Consider the first order system of the form a y& + y = u , its transfer function is Y ( s) =

bs + 1 s+z U (s) = K U (s) as + 1 s+ p

Y ( s) =

4.2 Response of the first order systems

Y ( s) = G ( s)U ( s) where Y (s ) : Laplace transform of the output G (s ) : transfer function of the system U (s ) : Laplace transform of the input

If the derivative of the input are involved in the differential equation of the system, that is, if a y& + y = b u& + u , then its transfer function is

(4.2)

K1 − K 2

K 2e − pt

On taking the inverse Laplace of equation (4.2), we obtain

y (t ) =

−

1 {

steady − state part

−t / a e1 23

(t ≥ 0)

t

Fig. 4.2

(4.3)

transient part

Both of the input and the output to the system are shown in Fig. 4.1. The response has an exponential form. The constant a is called the time constant of the system.

We note that the responses to the systems (Fig. 4.1 and Fig. 4.2) have the same form, except for the constant terms K1 and

K 2 . It appears that the role of the numerator of the transfer function is to determine these constants, that is, the size of y (t ) , but its form is determined by the denominator.

u (t )

1.00

4.3 Response of second order systems An example of a second order system is a spring-dashpot arrangement shown in Fig. 4.3. Applying Newton’s law, we find

0.63

M &y& = − µ y& − k y + u (t ) t

a

Fig. 4.1

where k is spring constant, µ is damping coefficient, y is the distance of the system from its position of equilibrium point, and it is assumed that y (0) = y& (0) = 0 .

___________________________________________________________________________________________________________ 16 Chapter 4 Transient and Steady State Response Analysis

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

y

Fig. 4.3

the form p1, 2 = α ± i β where α = a1 / 2 and β =

U (s) =

K s 2 + a1 s + a 2

where

K s ( s + p1 )( s + p 2 )

where p1, 2 =

(4.6)

a1 ± a12 − 4a 2 2

transfer function G ( s ) =

K1 =

K

α2 + β2

, K2 =

K (− β − iα ) 2 β (α 2 + β 2 )

, K3 =

K (− β + iα ) 2 β (α 2 + β 2 )

(Notice that K 2 and K 3 are complex conjugates) It follows that y (t ) = K1 + K 2 e −(α +iβ )t + K 3 e −(α −iβ )t

, p1 and p 2 are the poles of the K

2

s + a1 s + a 2

= K1 + e −αt [( K 2 + K 3 ) cos β t + ( K 3 − K 2 )i sin β t ] (using the relation e iβ t = cos β t + i sin β t )

, that is, the zeros of the =

denominator of G(s). There are there cases to be considered:

=

Case 1: a12 > 4a 2 → over-damped system

In this case p1 and p 2 are both real and unequal. Eq. (4.6) can be written as Y ( s) =

4a 2 − a12 .

K3 K1 K2 + + , s s + p1 s + p 2

Y ( s) =

U (s)

where K = 1 / M , a1 = µ / M , a 2 = k / M . Applying a unit step input, we obtain Y (s) =

1 2

Hence

On taking Laplace transforms, we obtain

M s2 + µ s + k

(4.10)

In this case, the poles p1 and p 2 are complex conjugate having

u (t ) = M &y& + µ y& + k y

1

(1 − e − pt − p t e − pt )

Case 3: a12 < 4a 2 → under-damped system

u (t )

Y ( s) =

p2

µ

M

Hence

K

y (t ) =

k

K3 K1 K2 + + s s + p1 s + p 2

K

α2 + β2 K 2

α +β

2

+ −

K

α2 + β2 K 2

α +β

2

e −αt (− cos β t − 1+

α2 β2

α sin β t β

(4.11)

e −αt sin( β t + ε ) (4.12)

where tan ε = β / α Notice that when t = 0 , y (t ) = 0 . The there cases discussed above are plotted in Fig. 4.4.

(4.7)

y (t )

where

1

under damped system u (t )

K K K K = , K2 = , K3 = K1 = p1 p 2 a 2 p1 ( p1 − p 2 ) p 2 ( p 2 − p1 )

critically damped system over damped system

(notice that K1 + K 2 + K 3 = 0 ). On taking Laplace transform of Eq.(4.7), we obtain y (t ) = K1 + K 2 e − p1t + K 3 e − p2t

The

transient

part

of

the

t

Fig. 4.4

(4.8) solution

is

seen

to

be

expect that when the damping is 0 (that is, a1 = 0 ) the system

K 2 e − p1t + K 3 e − p2t .

should oscillate indefinitely. Indeed when a1 = 0 , then

Case 2: a12 = 4a 2 → critically damped system

In this case, the poles are equal: p1 = p 2 = a1 / 2 = p , and K3 K K2 Y ( s) = = 1 + + s s + p ( s + p) 2 s (s + p) 2 K

From Fig. 4.4, we see that the importance of damping (note that a1 = µ / M , µ being the damping factor). We would

(4.9)

α = 0 , and β = a 2 and since sin ε = 1 and cos ε = 0 , then ε = π / 2 , Eq. (4.12) becomes y (t ) =

Hence y (t ) = K1 + K 2 e − pt + K 3 t e − pt , where K1 = K / p 2 , K 2 = − K / p 2 and K 3 = − K / p so that

[

π ⎞⎤ K K ⎡ ⎛ 1 − cos a 2 t ⎢1 − sin ⎜ a 2 t + ⎟⎥ = a2 ⎣ 2 ⎠⎦ a 2 ⎝

]

This response of the undamped system is shown in Fig.4.5.

___________________________________________________________________________________________________________ 17 Chapter 4 Transient and Steady State Response Analysis

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

(1) Overshoot defined as maximum overshoot × 100% final desired value

y (t )

(2) Time delay t d , the time required for a system response to reach 50% of its final value. 2π

4π

a2

a2

0

(3) Rise time, the time required for the system response to rise from 10% to 90% of its final value.

t

Fig. 4.5 There are two important constants associated with each second order system. • The undamped natural frequency ω n of the system is the frequency of the response shown in Fig. 4.5: ω n = a 2 • The damping ratio ξ of the system is the ratio of the actual damping µ ( = a1M ) to the value of the damping

µ c , which results in the system being critically damped a1 µ = . (that is, when a1 = 2 a 2 ). Hence ξ = µ c 2 a2 We can write equation (4.12) in terms of these constants. We note that a1 = 2ω nξ and a 2 = ω n2 . Hence

1+α 2 / β 2 = 1+

a12 4a 2 − a12

2 a2

=

4a 2 − a12

=

1

(4.13)

In general, the quality of a system may be judged by various standards. Since the purpose of a servomechanism is to make the output follow the input signal, we may define expressions which will describe dynamics accuracy in the transient state. Such expression are based on the system error, which in its simplest form is the difference between the input and the output and is denoted by e(t ) , that is, e(t ) = y (t ) − u (t ) , where y (t ) is the actual output and u (t ) is the desired output ( u (t ) is the input).

1−ξ 2

ξ

e (t )dt

0

.

the damping ratio ξ . There typical graphs are shown in Fig. 4.6. Some definitions

steady-state error ess

Having chosen an appropriate performance index, the system which minimizes the integral is called optimal. The object of modern control theory is to design a system so that it is optimal with respect to a performance index and will be discussed in the part II of this course. 4.4 Response of higher order systems

We can write the transfer function of an n th - order system in the form

G (s) =

0.5

0

e 2 (t )dt

∞

the ‘normalised’ response y (t ) against ω t for various values of

1.0 0.9

∞

0 ∞

∫ t e (t )dt

It is conventional to choose K / a 2 = 1 and then plot graphs of

maximum overshoot

(2) integral of absolute error (IAS)

∫ ∫

(3) integral of time multiplied absolute error criterion (ITAE)

where

y (t )

In fact, one can often improve one of the parameters but at the expense of the other. For example, the overshoot can be decreased at the expense of the time delay.

(1) integral of error squared (IES)

⎛ ⎞ K ⎜ 1 ⎟ 1− e −ω nξ t sin(ω t + ε ) ⎟ 2 ⎜ 2 ωn ⎜ ⎟ 1−ξ ⎝ ⎠

ω = ω n 1 − ξ 2 and tan ε =

(5) Steady-state error ess , the difference between the steady state response and the input.

The expression called the performance index can take on various forms, typically among them are:

1−ξ 2

Eq. (4.12) becomes

y (t ) =

(4) Settling time, the time required for the eventual settling down of the system response to be within (normally) 5% of its final value.

K ( s m + b1 s m −1 + K + b m ) s n + a1 s n −1 + K + a n

(4.14)

Example 4.1________________________________________

rise time

With reference to Fig. 2.11, calculate the close loop transfer 1 function G ( s ) given the transfer functions A( s ) = and td t s+3 B( s) = 2 / s Fig. 4.7 ___________________________________________________________________________________________________________ 18 Chapter 4 Transient and Steady State Response Analysis 0.1

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

R (s ) E (s )

C (s ) R (s )

A(s)

G ( s )=

C (s )

A( s ) 1 + A( s ) B ( s )

B(s)

Y ( s) = K (c1 + c 2 s + K + c s n−m−1 ) +

We obtain s s 2 + 3s + 2

=

s ( s + 1)( s + 2)

__________________________________________________________________________________________

The response of the system having the transfer function (4.14) to a unit step input can be written in the form

Y ( s) =

K ( s + z1 )( s + z 2 ) L ( s + z m ) s ( s + p1 )( s + p2 )L ( s + p n )

(4.15)

where z1 , z 2 ,L, z m : the zeros of the numerator

Example 4.2________________________________________ Find the response of the system having a transfer function 5( s 2 + 5s + 6) 3

s + 6 s 2 + 10 s + 8

to a unit step input.

We first assume that n ≥ m in equation (4.14); we then have two cases to consider: Case 1: p1 , p 2 ,L, p n are all distinct numbers. The partial fraction expansion of equation (4.15) has the form

K K1 K2 + + L + n+1 s s + p1 s + pn

The inverse Laplace transform of the first right term of (4.20) involves the impulse function and various derivatives of it. The second term of (4.20) is treated as in Case 1 or Case 2 above.

G(s) =

p1 , p 2 ,L, p n : the zeros of the denominator

Y ( s) =

(4.16)

In this case, 5( s 2 + 5s + 6)

Y ( s) =

3

s + 6 s 2 + 10 s + 8

Y ( s) =

K3 K1 K K4 + 2 + + s s + 4 s + (1 + i ) s + (1 − i )

where K1 =

y (t ) = K1 + K 2 e − p1t + K + K n+1e − pnt

Hence

Case 2: p1 , p 2 ,L, p n are not distinct any more. Here at least

one of the roots, say p1 , is of multiplicity r , that is

K ( s + z1 )( s + z 2 ) L ( s + z m ) s( s + p1 ) ( s + p 2 )L ( s + pn )

K n −r + 2 K1 K K 2r + 21 + L + +L+ s s + p1 s + p n−r +1 ( s + p1 ) r

(4.18)

⎧ K ⎫⎪ K t j −1e − pt , ( j = 1,2,L, r ) , the ⎨ ⎬= ⎪⎩ ( s + p ) j ⎪⎭ ( j − 1)! response has the form

Since L

−7 + i −7 − i 15 1 , K4 = . , K 2 = − , K3 = 4 4 4 4

15 1 −4t 1 −t − e + e (−14 cos t + 2 sin t ) 4 4 4 15 1 −4t 2000 −t − e + e sin(t + 352) 4 4 4

__________________________________________________________________________________________

4.5 Steady state error

The partial fraction expansion of equation (4.17) has the form

Y ( s) =

y (t ) = =

(4.17)

r

5( s + 2)( s + 3) s( s + 4)[ s + (1 + i )][s + (1 − i )]

=

The partial fraction expansion as

K1 , K 2 ,L, K n+1 are called the residues of the expansions. The response has the form

Y ( s) =

K1 ( s n + d1s n−1 + K + d n )

s ( s n + a1s n−1 + K + a n ) (4.20) where c s , d s , K and K1 are all constants.

Fig. 2.11

G(s) =

such as resistanceless circuits. We then divide the numerator until we obtain a proper fraction so that when applying a unit step input, we can write Y ( s ) as

Consider a unity feedback system as in Fig. 4.8

R (s )

E (s )

A(s)

C (s )

−1 ⎪

y (t ) = K1 + K 21e − p1t + K 22 e − p1t + K +

K 2 r r −1 − p1t t e + (r − 1)!

K 3e − p2t + L + K n−r + 2 e − pn −r +1t

(4.19)

We now consider that n < m in equation (4.14); which is the case when the system is improper; that is, it can happen when we consider idealized and physically non-realisable systems,

Fig. 4.8 where r (t ) : reference input c(t ) : system output e(t ) : error We define the error function as e (t ) = r (t ) − c (t )

(4.21)

___________________________________________________________________________________________________________ 19 Chapter 4 Transient and Steady State Response Analysis

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

hence, ess = lim e(t ) . Since E ( s ) = R( s ) − A( s ) E ( s ) , it

c(t )

t →∞

steady-state error ess

R( s) and by the final value theorem 1 + A( s ) sR ( s ) = lim sE ( s ) = lim (4.22) s →0 s →0 1 + A( s )

follows that E ( s ) = ess

k u (t )

We now define three error coefficients which indicate the steady state error when the system is subjected to three different standard reference inputs r (s ) .

t

Fig. 4.11

(1) step input, r (t ) = k u (t ) ( k is a constant)

Example 4.3________________________________________

ess = lim

sk /s k = A( s) 1 + lim A( s)

s →0 1 +

Find the (static) error coefficients for the system having a 8 open loop transfer function A( s ) = s ( 4 s + 2)

s →0

lim A( s ) = K p , called the position error constant, then

s→0

ess

K p = lim A( s) = ∞ s→0

k − ess k or K p = = ess 1+ K p

(4.23)

K v = lim sA( s ) = 4 s→0

K a = lim s 2 A( s ) = 0 s→0

c(t )

k u (t )

__________________________________________________________________________________________

steady-state error ess

From the definition of the error coefficients, it is seen that ess

k

depends on the number of poles at s = 0 of the transfer function. This leads to the following classification. A transfer function is said to be of type N if it has N poles at the origin. Thus if

t

Fig. 4.9

k 2

, so that ess =

(4.24)

K (− z1 )L (− z m ) where K1 = (4.25) (− p1 )L (− p n ) sj K1 is called the gain of the transfer function. Hence the steady K1

s →0

k k or K v = , where Kv ess

s K v = lim sA( s ) is called the velocity error constant. s→0

c(t )

K ( s − z1 ) L ( s − z m ) s j ( s − p1 )L ( s − p n )

At s = 0 , A( s ) = lim

(2) Ram input, r (t ) = k t u (t ) ( k is a constant) In this case, R ( s) =

A( s ) =

state error ess depends on j and r (t ) as summarized in Table 4.1 Table 4.1 j

steady-state error ess

0 1 2

k u (t )

System Type 1 Type 2 Type 3

ess r(t)=ku(t)

r(t)=ktu(t)

r(t)=½kt2u(t)

Finite 0 0

∞ finite 0

∞ ∞ finite

4.6 Feedback Control t

Fig. 4.10 1 2 k t u (t ) ( k is a constant) 2 k k k In this case, R ( s ) = 3 , so that ess = or K a = , where Ka ess s

(3) Parabolic input, r (t ) =

K a = lim s 2 A( s) is called the acceleration error constant. s→0

Consider a negative feedback system in Fig. 4.12 R (s )

E (s )

A(s)

C (s )

B(s)

Fig. 4.12 The close loop transfer function is related to the feed-forward transfer function A(s ) and feedback transfer function B (s) by

___________________________________________________________________________________________________________ 20 Chapter 4 Transient and Steady State Response Analysis

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

G(s) =

A( s ) 1 + A( s ) B ( s )

(4.26)

We consider a simple example of a first order system for K and B ( s ) = c , so that which A( s ) = as + 1

K /a Kc + 1 s+ a On taking Laplace inverse transform, we obtain the impulse response of the system, where G(s) =

K = as + Kc + 1

K (1) c = 0 (response of open loop system): g (t ) = e −t / a a a K − K ca+1 t K − αt (2) c ≠ 0 : g (t ) = e = e , where α = K c +1 a a a and α are respectively the time-constants of the open loop and closed loop systems. a is always positive, but α can be either positive or negative.

Fig. 4.13 shows how the time responses vary with different values of K c . g (t )

Kc ≤ −1

Kc = −5

Kc = 0 Kc = 4 t

If the impulse response does not decay to zero as t increase, the system is unstable. From the Fig. 4.13, the instability region is defined by Kc ≤ −1 . In many applications, the control system consists basically of a plant having the transfer function A(s ) and a controller having a transfer function B (s) , as in Fig. 4.14. E (s )

B(s)

Q (s)

controller

A(s)

C (s )

plant

where K p is a constant, called the controller gain. The transfer function of this controller is Q( s) = Kp B( s) = E (s)

(4.28)

(3) Integral Controller

∫

t

In this case q(t ) = K e(t )dt , hence 0

B( s) = K / s

(4.29)

(4) Derivative Controller de In this case q (t ) = K , hence dt (4.30)

⎛ K ⎞ B ( s ) = K p ⎜1 + 1 s ⎟ = K p (1 + K s ) ⎜ Kp ⎟ ⎝ ⎠

(4.30)

With the closed loop transfer function A( s) B ( s ) 1 + A( s ) B( s )

The controllers can be of various types. (1) The on-off Controller The action of such a controller is very simple. ⎧Q if e(t ) > 0 q (t ) = ⎨ 1 ⎩Q2 if e(t ) < 0

∫

t

In this case q(t ) = K p e(t ) + K1 e(t )dt , hence 0

⎛ K 1⎞ ⎛ K⎞ B ( s ) = K p ⎜1 + 1 ⎟ = K p ⎜ 1 + ⎟ ⎜ Kp s⎟ s⎠ ⎝ ⎝ ⎠

(4.31)

(7) Proportional-Derivative-Integral Controller (PID) t de In this case q (t ) = K p e(t ) + K1 + K 2 e(t )dt , hence dt 0

∫

⎛ K K 1⎞ B ( s) = K p ⎜1 + 1 s + 2 ⎟ = K p (1 + k1s + k 2 / s ) ⎜ Kp K p s ⎟⎠ ⎝ Example 4.4________________________________________

Fig. 4.14

G(s) =

(2) Proportional Controller For this control action q (t ) = K p e(t )

(6) Proportional-Integral Controller (PI)

Fig. 4.13

R (s )

The on-off controller is obviously a nonlinear device and it cannot be described by a transfer function.

(5) Proportional-Derivative Controller (PD) de In this case q (t ) = K p e(t ) + K1 , hence dt

Kc = −1

stable region

q (t ) is output signal from the controller Q1 , Q2 are some constants

B( s ) = Ks

Kc = −3

K a

where

(4.27)

Design a controller for a plant having the transfer function A( s ) = 1 /( s + 2) so that the resulting closed loop system has a zero steady state error to a reference ramp input. For zero steady state error to a ramp input, the system must be of type 2. Hence if we choose an integral controller with B ( s ) = K / s then the transfer function of the closed loop system including the plant and the controller is A( s ) B ( s ) K = 1 + A( s ) B ( s ) s 3 + 2 s 2 + K

___________________________________________________________________________________________________________ 21 Chapter 4 Transient and Steady State Response Analysis

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.3

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

The response of this control system depends on the roots of the denominator polynomial s 3 + 2 s 2 + K = 0 . If we use PI controller, B ( s) = K p (1 + K / s ) the system is of type 2 and response of the system depends on the roots of

s 3 + 2 s 2 + K p s + KK p = 0 __________________________________________________________________________________________

___________________________________________________________________________________________________________ 22 Chapter 4 Transient and Steady State Response Analysis

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.5 Stability 5.1 Introduction

G(s) =

5.2 The Concept of Stability

s n + a1s n−1 + K + a n = 0

A system is said to be asymptotically stable if its weighting function response decays to zero as t tends to infinity.

(5.2)

Eq. (5.2) is called characteristic equation of the system. It follows that the system is asymptotically stable if and only if the zeros of the characteristic equation (that is the finite poles of the transfer function) p1 , p 2 ,L, p n are negative or have negative real parts. Fig. 5.1 illustrates graphically the various cases of stability.

If the response of the system is indefinitely oscillatory, the system is said to be marginally stable. Back to the system transfer function has the form

s-plane graph ( s = σ + jω ) ω

Response graph

Remark

y

Asymptotically stable

Real and negative

σ

t

ω

y

Unstable

Real and positive

σ

t

ω

y

Marginally stable

Zero

σ Conjugate complex with negative real part

t

ω

y

σ ω

Conjugate imaginary (multiplicity r=1)

ω

Asymptotically stable

t

Marginally stable

t

Unstable

t

Unstable

y

σ ω

Conjugate with positive real part

t

y

σ

Conjugate imaginary (multiplicity r=2)

y

σ ω

Roots of multiplicity r=2 at the origin

(5.1)

s n + a1s n−1 + K + a n

and the system response is determined, except for the residues, by the poles of G (s ) , that is by the solution of

The properties of a system are characterized by its weighting function g (t ) , that is, its response to a unit impulse or equivalently by the Laplace transform of the weighting function – the transfer function.

Type of root

K ( s m + b1s m−1 + K + bm )

y

y=t

σ

Unstable t

Fig. 5.2 5.3 Routh Stability Criterion 5.4 Introduction to Lyapunov’s Method

Lyapunov’s method is the most general method we know for system stability analysis. It is applicable for both linear and nonlinear systems of any order. For linear system, it provides both the necessary and sufficient conditions, whereas for

___________________________________________________________________________________________________________ 23 Chapter 5 Stability

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

nonlinear systems it provides only the sufficient conditions for asymptotically stable. Lyapunov considered the solution x(t ) of the state equations to trace out a trajectory (or an orbit or a path) which, for an asymptotically stable system, terminates at the origin of the n-dimensional state-space, known as the phase plane when n=2.

x12 +

1 g /l

x 22 = A 2

(5.9)

The trajectories (for various values of A and ε ) are ellipses as shown in Fig. 5.3 x2

Take a simple example: consider the motion of an undamped pendulum. The equation of motion is

θ&& +

g sin θ = 0 l

x1

(5.4) Fig. 5.3

which is nonlinear system.

θ

Every trajectory shown in Fig. 5.3 is closed, showing that the pendulum, assumed to have no resistance to its motion, will swing indefinitely. A closed trajectories is typical of a periodic motion, that is, one for which the solution x(t ) has the property x(t + T ) = x(t ) , where T is the period of the motion. l

When the oscillations are damped, the trajectories would terminate at the point x = 0 so that it may have the form shown in Fig. 5.4. x2

mg

Fig. 5.2

x0 x1

Assume that the displacements are small, we have sin θ ≈ θ and (5.4) can be rewritten as Fig. 5.4 g l

θ&& + θ = 0

(5.5)

This is well known linear equation of simple harmonic motion, and has the solution

θ = A sin( g / l t + ε )

In general autonomous system, the state equation correspondding to equation (5.7) takes the form

⎧ x&1 = P ( x1 , x2 ) ⎨ ⎩ x& 2 = Q( x1 , x2 )

(5.9)

(5.6)

The solutions corresponding to Eq. (5.8) are where A and ε are constants determined by the initial conditions. Let define the state variables: x1 = θ , x 2 = θ& . (5.5) can be written as

⎧ x&1 = x 2 ⎪ ⎨ g ⎪ x& 2 = − l x1 ⎩

and obtain the trajectory by eliminating t in (5.8)

(5.10)

and can be plotted to give the trajectory. For an n-dimensional state equation, the algebra is similar.

(5.7)

Definition

A point ( x10 , x20 ) , for which P ( x10 , x20 ) = 0 = Q( x10 , x20 ) is called a critical point.

and the solution

⎧⎪ x1 = A sin( g / l t + ε ) ⎨ ⎪⎩ x 2 = A g / l cos( g / l t + ε )

⎧ x1 = φ (t ) ⎨ ⎩ x2 = ψ (t )

(5.8)

For a linear system, there is only one equilibrium point and it is the point x = 0 . A system is asymptotically stable when it is stable and the trajectory eventually approaches x = 0 as t → ∞ . In mathematical terminology, these definitions take the form shown in Fig. 5.5

___________________________________________________________________________________________________________ 24 Chapter 5 Stability

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

The equilibrium point x = 0 (Fig. 5.5) is said to be (1) stable if , given any t > 0 , there exists δ > 0 so that every solution of (5.10) which at t = 0 is such that x12 + x22 = ϕ 2 (0) + ψ 2 (0) < δ 2 implies

The potential energy U of the system U = mgl (1 − cosθ ) ≈ 12 mglθ 2 (for small θ )

ϕ 2 (0) + ψ 2 (0) < ε 2 for t ≥ 0

The kinetic energy of the system

(2) asymptotically stable if, (a) it is stable, and ⎧ lim ψ (t ) = 0 ⎪t →∞ (b) ⎨ or lim [ψ 2 (t ) + ψ 2 (t )] = 0 ψ (t ) = 0 t →∞ ⎪⎩tlim →∞

K = 12 ml 2θ& 2 Hence the system energy is V = 12 mglθ 2 + 12 ml 2θ& 2

(3) unstable if it is not stable

= 12 mglx12 + 12 ml 2 x22

x2

δ

Notice that V (0,0) = 0 . It follows that

ε

x12 x1

2

+

a V

x ( 0)

x22 b 2V

= 1 , where a 2 =

2 2 and b 2 = 2 mgl ml

The trajectory is (again) seen to be an ellipse (Fig. 5.6) with major and minor axes having lengths 2a 2V and 2b 2V respectively.

unstable

x2

Fig. 5.5

b 2V

If we consider an unforced system, having the state equation

14243

In the definition, starting ‘near’ the equilibrium point is defined by the ‘ δ -neighborhood’, whereas the ‘reasonable’ distance is the ‘ ε -neighborhood’.

x& = A x

(5.13)

= V ( x1 , x2 )

asymptotically stable

0

stable

Lyapunov’s method depends on the state trajectory – but not directly. We illustrate his method using the pendulum example.

0

1442443 a 2V

x1

(5.11)

we have found the solution x = e At x(0)

= ( A1e λ1t + A2 e λ2t + K + An e λnt ) x(0)

(5.12)

Fig. 5.6 The rate of change of the energy along a trajectory is

For asymptotic stability, lim x(t ) = 0 , where x(t ) is the t →∞

norm of the (solution) vector x defined by

x = (x, x) = x12 + x22 + K + xn2 Hence for asymptotic stability all the state variables x1 , x2 ,L, xn must decrease to zero as t → ∞ . Form (5.12),

we see that this is the case when Re λi < 0 (i = 1,2,L, n) , where λi is the eigenvalues of A , that is, the solution of the characteristic equation det( sI − A) = 0 .

dV (x) dV dx1 dV dx2 = + V& = dt dx1 dt dx2 dt = mglx1

(5.14)

dx1 dx + ml 2 x2 2 dt dt

= mglx1 x2 + ml 2 x2 (− g / l ) x1 =0 Hence the total energy is a constant along any trajectory. If dV / dt < 0 , the energy is decreasing monotonically along every trajectory. This means that every trajectory moves toward the point of minimum energy, that is the equilibrium point (0,0). This is an essence the Lyapunov’s criterion for asymptotic stability.

___________________________________________________________________________________________________________ 25 Chapter 5 Stability

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Lyapunov’s function V (x) is similar to the energy function, it has the following properties: (1) V (x) and its partial derivatives ∂V / ∂xi (i = 1,2,L, n) , are continuous. (2) V (x) > 0 for x ≠ 0 in the neighborhood of x = 0 , and V ( 0) = 0 (3) V& (x) < for x ≠ 0 (in the neighborhood of x = 0 ), and

V& (0) = 0 V (x) satisfying the above properties is called a Lyapunov function. Definition V (x) is • positive-definite if V (x) > 0 for x ≠ 0 , and V (0) = 0 • positive-semi-definite if V (x) ≥ 0 for all x • negative-definite if −V (x) is positive-definite • negative-semi-definite if −V (x) is positive-semi-definite

V ( x) = x ′ P x ⎡ p11 p12 ⎢p p = [x1 x 2 L x n ]⎢ 21 22 M M ⎢p ⎣ n1 p n 2 n

=

L p12 ⎤ ⎡ x1 ⎤ L p 2n ⎥ ⎢ x2 ⎥ O M ⎥⎢ M ⎥ L p nn ⎥⎦ ⎢⎣ xn ⎥⎦

(5.15)

n

∑∑ p x x

ij i j

i =1 j =1

The matrix P is of order n × n and is symmetric. Sylvester’s Criterion: the necessary and sufficient conditions for a quadratic form V ( x) = x′ P x to be positivedefinite is that all the successive principal minors of the symmetric matrix P are positive, that is

p11 > 0 ,

p11 p12 p13 p11 p12 > 0 , p 21 p 22 p 23 > 0 , etc. p 21 p 22 p31 p32 p33

Example 5.6_______________________________________ Example 5.5_______________________________________

Consider the function

Classify the Lyapunov function

V (x) = 4 x12 + 3 x 22 + x32 + 4 x1 x 2 − 2 x1 x3 − 2 x 2 x3

a. V ( x) = x12 + 3 x22

Rewrite V (x) in the form

V (x) > 0, ∀x ≠ 0 and V (0) = 0 ⇒ V (x) is positive-definite b. V (x) = −( x1 + x2 ) 2 − 2 x12 V (x) < 0, ∀x ≠ 0 and V (0) = 0 ⇒ V (x) is negative-definite c. V ( x) = ( 2 x1 + x 2 ) 2 Whenever x2 = −2x1 , V (0) = 0 , otherwise V (0) > 0 . Hence V (x) ≥ 0 ⇒ V (x) is positive-semi-definite

V (x) = 4 x12 + 3 x 22 + x32 + 4 x1 x 2 − 2 x1 x3 − 2 x 2 x3 ⎡ 4 2 − 1⎤ ⎡ x1 ⎤ = [x1 x 2 x3 ]⎢ 2 3 − 1⎥ ⎢ x2 ⎥ ⎢ ⎥ ⎣⎢− 1 − 1 1 ⎥⎦ ⎣ x3 ⎦ 4 2 −1 Since 4>0, 4 2 = 8 > 0 , 2 3 − 1 = 5 > 0 , Sylvester’s 2 3 −1 −1 1 criterion assures us that V (x) is positive-definite. _________________________________________________________________________________________

d. V (x) = x12 + 2 x 22 − 4 x1 x 2 V (1,1) < 0 and V (1,3) > 0 . V (x) assumes both positive and negative values ⇒ V (x) is indefinite _________________________________________________________________________________________

Lyapunov’s Theorem The origin (that is the equilibrium point) is asymptotically stable if there exists a Lyapunov function in the neighborhood of the origin. If V& (x) ≤ 0 , then the origin is a stable point.

This is sometimes referred to as stability in the sense of Lyapunov. 5.5 Quadratic Forms Definition A quadratic form is a scalar function V (x) of variables x = [x1 x 2 L x n ]T defined by

5.6 Determination of Lyapunov’ Functions Definition The matrix A is said to be stable, if the corresponding system defined by equation x& = A x is stable. As explained above, to determine whether A is stable we seek a symmetric matrix P so that V ( x) = x′ P x is a positivedefinite function, and

V& ( x) = x& ′ P x + x′ P x& = ( A x)′ P x + x′ PA x = x′[ A′ P + PA] x

(5.16)

is negative. Since P is symmetric, [ A′ P + PA] is symmetric. Hence V& ( x) is a quadratic form. Let −Q = A′P + PA

(5.17)

for V& ( x ) to be negative-definite, x′ Q x must be positivedefinite.

___________________________________________________________________________________________________________ 26 Chapter 5 Stability

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Lyapunov’s theorem can now be expressed in the following form: Given a positive-definite matrix Q , the matrix A is stable if the solution to the equation (5.17) results in a positivedefinite matrix P .

Example 5.7_______________________________________ Use Lyapunov’s direct method to determine whether the following systems are asymptotically stable: ⎧ x& = x 2 a. ⎨ 1 ⎩ x& 2 = 2 x1 − x 2 Solving A′P + PA = − I , we obtain 1 P = ⎡ −3 −1⎤ 4 ⎢⎣− 1 1 ⎥⎦ which, by Sylvester’s criterion, is not positive definite. Hence the system is not stable. ⎧ x&1 = x 2 b. ⎨ ⎩ x& 2 = −5 x1 − 2 x 2 1 We obtain P = ⎡17 1⎤ which is positive definite. Hence 10 ⎢⎣ 1 3⎥⎦ the system is asymptotically stable.

5.7 The Nyquist Stability Criterion

If the transfer function of a system is not known, it is often possible to obtain it experimentally by exciting the system with a range of sinusoidal inputs of fixed amplitude by varying frequencies, and observing the resulting changes in amplitude and phase shift at the output. The information obtained in this way, known as the system frequency response, can be used in various ways for system analysis and design. 5.8 The Frequency Response

Consider the system Y ( s ) = G ( s ) U ( s ) . For a sinusoidal input of unit amplitude u (t ) = sin ω t so that U ( s) = ω /( s 2 + ω 2 ) and Y ( s0 = G (s)

Determine using Lyapunov’s method, the range of K for the system in Fig. 5.8 to be asymptotically stable U (s )

K s +1

X 2 (s)

3 s+2

X 1 (s)

Y ( s ) can be expanded in partial fractions as n

Y (s) =

From Fig. 5.8 3 X1 = X2 s+2 K X2 = (− X 1 + U ) s +1

⇒

x&1 = −2 x1 + 3 x 2

⇒

x& 2 = − Kx1 − x 2 + Ku

Kj

j =1

Solving A′P + PA = − I

⎡− 2 − K ⎤ ⎡ p11 p12 ⎤ + ⎡ p11 p12 ⎤ ⎡ − 2 3 ⎤ = ⎡− 1 0 ⎤ ⎢⎣ 3 − 1 ⎥⎦ ⎢⎣ p 21 p 22 ⎥⎦ ⎢⎣ p 21 p 22 ⎥⎦ ⎣⎢− K − 1⎥⎦ ⎣⎢ 0 − 1⎥⎦ 1 ⎡ K 2 + 3K + 3 3 − 2 K ⎤ 3K + 15⎥⎦ 18 K + 12 ⎢⎣ 3 − 2 K

Sylvester’s criterion leads us to conclude that K > −2 / 3 .

+ j

Q1 Q2 − s + iω s − iω

(5.18)

and Kj =

Q1 =

(s + p j ) ω G(s) s2 + ω 2

s =− p j

(s + p j ) ω G(s)

s2 + ω 2

=− s = −i ω

(s + p j ) ω G(s)

s2 + ω 2

= s =i ω

1 G (−i ω ) 2i

1 G (i ω ) 2i

On taking the inverse Laplace transform of equation (5.18), we obtain n

The state equation is ⎡ x&1 ⎤ ⎡ − 2 3 ⎤ ⎡ x1 ⎤ ⎡ 0 ⎤ ⎢⎣ x& 2 ⎥⎦ = ⎢⎣− K − 1⎥⎦ ⎢⎣ x 2 ⎥⎦ + ⎢⎣ K ⎥⎦ u

We obtain P =

∑s+ p

Q2 =

Fig. 5.8

s +ω2

Assuming that G (s ) has poles at s = − p1 ,− p 2 , L ,− p m ,

_________________________________________________________________________________________

Example 5.8_______________________________________

ω 2

y (t ) =

∑K e j

− p jt

j =1 1 42 4 43 4 transient form

+ Q1e −iω t + Q2 e i ω t

(5.19)

144424443 steady − state form

The steady-state term can be written as 1 [G (iω ) e iω t − G (−iω ) e −iω t ] = Im{G (iω ) e iω t } 2i

(5.20)

In polar form we can express G (iω ) = G (iω ) e iφ , where

φ = tan −1

Im{G (iω )} , and Eq. (5.12) becomes Re{G (iω )}

Im{ G (iω ) e i (ω t +φ ) } = G (iω ) sin(ωt + φ )

(5.21)

_________________________________________________________________________________________

___________________________________________________________________________________________________________ 27 Chapter 5 Stability

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

(5.21)

Eq. (5.21) shows that the steady-state output is also sinusoidal, and relative to the input the amplitude is multiplied by G (iω ) . The frequency diagram is the plot of G (iω ) in the complex plane as the frequency varies form −∞ to +∞ . Since for the systems considered the transfer function G (s ) are the ratios of two polynomials with real coefficients, it follows that G (iω ) = the complex conjugate of G (−iω ) and it is sufficient to plot G (iω ) for 0 ≤ ω < ∞ ; the remainder of the locus is symmetrical about the real axis to this plot. Example 5.9_______________________________________

Determine the frequency diagrams for the following transfer functions

where

φ j (k ) = tan −1

Im{G j (iω )}

( j = 1,2, L , k )

Re{G j (iω )}

Example 5.10______________________________________

Construct the frequency diagram for the system having the transfer function

G(s) =

K ( K is constant) s ( s + 1)

⎛ 1 ⎞⎛ K ⎞ Since G (iω ) = ⎜ ⎟ ⎜ ⎟ = G1 (iω ) G2 (iω ) , we obtain ⎝ iω ⎠ ⎝ iω + 1 ⎠ 1 1 G1 (iω ) = and φ = −π / 2 − tan −1 ω . The plot is ω 1+ω2 given in Fig. 5.10 Imaginary axis

=0

Im{ G (iω ) e i (ω t +φ ) } = G (iω ) sin(ωt + φ )

ω

a. G ( s ) = Ks

ω = −∞

G (iω ) = iKω , hence G (iω ) = Kω and ϕ (ω ) = 90 0 ⇒ the plot is the positive imaginary axis (for ω ≥ 0 )

Real axis

b. G ( s ) = K / s

ϕ

G

G (iω ) = K /(iω ) , hence G (iω ) = K / ω and ϕ (ω ) = −90

ω = +∞

0

ω =0

⇒ the plot is the negative imaginary axis (for ω ≥ 0 ) c. G ( s ) = K /(a + s )

Fig. 5.10 ________________________________________________________________________________________

G (iω ) = K /( a + iω ) = K (a − iω ) /(a 2 + ω 2 ) G (iω ) = K / a 2 + ω 2

and

,

hence

ϕ (ω ) = − tan −1 (ω / a)

⇒ the plot for ω ≥ 0 is given in Fig. 5.9

To use the Nyquyist criterion we must plot the frequency response of the system’s open-loop transfer function. Consider the closed loop system in the Fig. 5.11

Locus corresponding to negative frequency

Imaginary axis

R (s )

E (s )

G(s)

C (s )

C1 ( s ) = H ( s ) C ( s )

ω = −∞ K 2a

0

ϕ ω = +∞

K a

G (iω

ω =0

H(s)

Real axis

)

increasing ω

Fig. 5.11 The closed loop transfer function C (s) G ( s) = G( s) 1 + G ( s) H ( s)

(5.22)

Fig. 5.9 _________________________________________________________________________________________

If the transfer function involves several factors, the above procedure can be used for each factor individually, and the overall result is obtained by the following rules If G ( s ) = G1 ( s ) G 2 ( s ) L G k ( s ) , then

G (iω ) = G1 (iω ) G2 (iω ) L Gk (iω ) , and

φ (ω ) = φ1 (ω ) + φ 2 (ω ) + K + φ k (ω ) ,

The ratio of the feedback signal C1 ( s) to the actuating error

C1 ( s ) = G ( s ) H ( s ) and is called the system E (s) open-loop transfer function. signal E ( s ) is

If we consider the unity negative feedback system shown in C (s) = G(s) H (s) Fig. 5.12 we again find that 1 E (s)

___________________________________________________________________________________________________________ 28 Chapter 5 Stability

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

R (s )

E (s )

s = 1 + iω

C (s )

G(s)H(s)

and corresponding along D ' A'

C1( s )

F ( s) =

Fig. 5.12 And both of the systems in Fig. 5.11 and Fig. 5.12 have the same open-loop transfer functions. Assume that if there are any factors which cancel in the product G ( s) H ( s) they involve only poles and zeros in the left haft s-plane. We will consider the system to be stable so long as the close loop transfer function equation (5.22) has no poles in the RHP. We note that s = s1 is a pole of equation (5.22) if (1) 1 + G ( s1 ) H ( s1 ) = 0 and G ( s1 ) ≠ 0 (2)

G(s) =

A( s) ( s − s1 ) n

and 1 + G ( s ) H ( s ) =

B( s) ( s − s1 ) m

, where

n > m A( s1 ) ≠ 0 and B ( s1 ) ≠ 0 . From (1) and (2), we conclude that the system having the transfer function equation (5.22) is stable if and only if [1 + G ( s ) H ( s )] has no zeros in RHP.

5.9 An Introduction to Conformal Mappings To establish the Nyquist criterion we regard the open loop transfer function G ( s ) H ( s ) = F ( s ) of a feedback system as a mapping from the s-plane to the F ( s) -plane. As an example, consider F ( s ) = 1 /( s + 4) . We consider two planes: (1) the s-plane, having a real axis, Re{s} = σ and imaginary axis Im{s} = ω , where s = σ + iω (2) the F ( s) -plane which has similarly a real axis Re{F ( s)} and imaginary axis Im{F ( s )} Consider corresponding paths in the two planes as in Fig. 5.13 s − plane

F ( s ) − plane

iω

B

Im{F ( s )}

i

A

C’

0 .1

D’ 1 -1

0

σ

0 .1

0 .2

0 .3

Re{F ( s )}

A’ C

D

-i

γ − path

−1 ≤ ω ≤ 1

− 0 .1 Γ − path

B’

Fig. 5.13 In Fig. 5.13 it is shown that F ( s ) traces out circular arcs as s moves along the insides of the square. This can of course be proved analytically. For example, along DA

1 5 − iω = = u + iv 5 + iω 25 + ω 2

where u =

5 25 + ω

2

and v = −

ω . They relationship is 25 + ω 2

(u − 1 / 10) 2 + v 2 = (1 / 10) 2 This is the equation of a circle, center (1 / 10,0) and radius 1 / 10 . In Fig. 5.13 we note that (1) the closed contour γ in the s -plane is mapped into a closed contour Γ in the F ( s) -plane (2) Any point inside γ maps into a point inside Γ (3) We can define a positive direction in which s traces out its path. It does not matter which direction is chosen to be positive. It is interesting to consider the direction in which F ( s ) sweeps out its path, as s moves along a contour which encloses the pole ( s = −4) of F ( s ) . This is shown in Fig. 5.14. s − plane

B

iω

Im{F ( s )}

A

C’

i

F ( s ) − plane

D’

σ -5

-4

-3

-2 -i

C

D

Re{F ( s )}

0

B’

γ − path

A’ Γ − path

Fig. 5.14 In this case we note that as s traces the γ -contour in the positive direction, F ( s) traces the corresponding Γ -contour in the opposite, that is, negative direction..

Cauchy’s Theorem Let γ be a closed contour in the s -plane enclosing P poles and Z zeros of F ( s ) , but not passing through any poles or zeros of F ( s ) . As s traces out γ in the positive direction, F ( s) traces out a contour Γ in the F ( s ) − plane, which encircles the origin ( Z − P ) times in the positive direction. In Fig. 5.13, γ does not enclose any poles or zeros of F ( s ) , so that P = 0 = Z ⇒ the origin of the F ( s ) -plane is not encircled by the Γ -contour. In Fig. 5.14, γ encircles a pole of F ( s ) , so that P = 1 , Z = 0 , hence Z − P = −1 ⇒ Γ -contour encircles the origin (-1) times in the positive direction, that is, Γ encircles the origin once in the negative direction.

___________________________________________________________________________________________________________ 29 Chapter 5 Stability

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

5.10 Application of Frequency Response

Conformal

Mappings

to

the

In section 5.8 we discussed the frequency response, which required s to trace out the imaginary axis in the s-plane. In section 5.9, Cauchy’s theorem calls for a closed contour in the s-plane avoiding the poles and zeros of F ( s ) . We now combine these two contours and define the Nyquist path. Nyquist path consist of the entire imaginary axis, and a semicircular path of large radius R in the right half of s-plane, as shown in Fig. 5.15 s − plane

F ( s) =

K s ( s + 1)( s + 2)

when (1) K = 1 and (2) K = 10 Since F ( s) has a pole at the origin, the Nyquist contour has a semicircular indentation as shown in Fig. 5.18(a)

F ( s ) − plane

Nyquist path

R→

Use the Nyquist criterion to determine the stability of the closed-loop system having an open-loop transfer function

s − plane iω i∞

Im{F ( s )}

∞

i 0+

σ

Re{F ( s )}

0

r

θ Γ

γ

i0

The small semicircular indentations are used to avoid poles and zeros of F ( s) which lie along the imaginary axis. The semi circle in the s-plane is assumed large enough to enclose all the poles and zeros of F ( s) which may lie in the RHP.

The Nyquist Criterion To determine the stability of a unity feedback system, we consider the map Γ in the F ( s) -plane of the Nyquist contour traversed in a positive direction. The system is stable if (1) Γ does not encircle the point (-1,0) when the open-loop transfer function G ( s ) H ( s ) has no poles in the RHP. (2) Γ encircles the point (-1,0) P times in the positive direction when the open-loop transfer function G ( s ) H ( s ) has P poles in the RHP. Otherwise, the system is unstable.

ω = 2 ω = +∞

α σ

0.2

0.1

−

−i∞

Fig. 5.15

F ( s ) − plane ω = −0 F (i 0 − )

R

iω

Example 5.11______________________________________

0

ω = −∞

F (i 0 + )

ω = +0

(a)

(b)

Fig. 5.18 The F ( s) contour for K = 1 corresponding to the Nyquist path is shown in Fig. 5.18(b). (1) For K = 1 , the point (−1,+∞) is not encircled by the F ( s ) contour. Since F ( s ) has no poles in the RHP, we can conclude that the system is stable. (2) For K = 10 , the magnitude F (iω ) is 10 times larger than when K = 1 for each finite value of ω while the phase φ is unchanged. This means that Fig. 5.18 is valid for K = 10 if we change the coordinate scales by a factor of 10. This time the critical point (−1,+∞) is not encircled twice. We can conclude that there are two closed-loop poles in the RHP, and the system is unstable. ________________________________________________________________________________________

Fig. 5.7 shows examples of a stable and an unstable system for the case where P = 0 . F ( s ) − plane

-1

F ( s ) − plane

0

-1

stable system

0

unstable system

Fig. 5.17

___________________________________________________________________________________________________________ 30 Chapter 5 Stability

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.6 Controllability and Observability 6.1 Introduction It is seen from equation (6.6) that if b′i , the i th row of B1 , has all zero components, then z&i = λi z i + 0 , and the input

6.2 Controllability We consider a system described by the state equations ⎧ x& = A x + B u ⎨ ⎩y = C x

(6.1)

where A is n × n , B is n × m and C is r × n .

By definition, the system equation (6.6) is controllable if only if a control function u(t ) can be found to transform the initial

With the transformation

x = Pz

u(t ) has no influence on the i th mode (that is e λit ) of the system. The mode is said to be uncontrollable, and a system having one or more such modes is uncontrollable. Otherwise, where all the modes are controllable the system is said to be completely state controllable, or just controllable.

state (6.2)

z 0 = [z 01 z 02 L z 0 n ]T

to

a

specifies

state

z1 = [z11 z12 L z1n ] . T

we can transform equation (6.1) into the form Example 6.1_______________________________________ ⎧ x& = A1z + B1u ⎨ ⎩y = C1x

(6.3)

where A1 = P −1 A P , B1 = P −1 B and C1 = C P . Assuming that A has distinct eigenvalues λ1 , λ2 ,L, λn we can choose

Check whether representation

the

system

having

the

state-space

1 2⎤ ⎡ 4⎤ x& = ⎡− ⎢⎣ 3 4⎥⎦ x + ⎢⎣6⎥⎦ u y = [1 − 2]x

P so that A1 is a diagonal matrix, that is,

is controllable. A1 = diag{λ1 , λ2 ,L, λn }

The eigenvalues λ1 = 1 , λ2 = 2 and the corresponding eigenvectors x1 = [1 2]T , x 2 = [2 3]T , so that the modal matrix is

If n = m = r = 2 the first order of (6.3) has the form ⎡ z&1 ⎤ ⎡λ1 0 ⎤ ⎡ z1 ⎤ ⎡ b11 b12 ⎤ ⎡ u1 ⎤ ⎢⎣ z& 2 ⎥⎦ = ⎢⎣ 0 λ2 ⎥⎦ ⎢⎣ z 2 ⎥⎦ + ⎢⎣b21 b22 ⎥⎦ ⎢⎣u 2 ⎥⎦

P = ⎡1 2⎤ and P −1 = ⎡ 3 −2⎤ ⎢⎣1 3⎥⎦ ⎢⎣− 1 1 ⎥⎦

which is written as ⎧ z& 1 = λ1z1 + b1′ u ⎨ ⎩z& 2 = λ2 z 2 + b ′2 u

Using transformation x = P z , the sate equation becomes (6.4)

where b′1 and b′2 are the row vectors of the matrix B1 . The output equation is

⎡ c 21 ⎤ ⎡ y1 ⎤ ⎡ c11 c12 ⎤ ⎡ z1 ⎤ ⎡ c11 ⎤ ⎢⎣ y 2 ⎥⎦ = ⎢⎣c21 c22 ⎥⎦ ⎢⎣ z 2 ⎥⎦ = ⎢⎣c21 ⎥⎦ z1 + ⎢⎣c 22 ⎥⎦ z 2 or y = c1z1 + c 2 z 2

This equation shows that the first mode is uncontrollable and so the system is uncontrollable. ________________________________________________________________________________________

Definition (6.5)

where c1 ,c 2 are the vectors of C1 . So in general, equation (6.3) can be written in the form

⎧z& i = λi z i + b ′i u (i = 1,2,L, n) ⎪⎪ n ⎨ ci z i y = ⎪ ⎪⎩ i =1

∑

z& = ⎡⎢1 0⎤⎥ z + ⎡⎢0⎤⎥ u ⎣0 2 ⎦ ⎣ 2⎦ y = [− 1 − 4]z

The matrix Q = [ B | A B | L | A n−1 B ] is called the system controllability matrix. The Controllability Criterion rank{Q} = n ⇔ the system is controllable Example 6.2_______________________________________

(6.6)

Using the controllability criterion verify that the system in example 6.1 is uncontrollable.

___________________________________________________________________________________________________________ Chapter 6 Controllability and Observability

31

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

For the system b = ⎡⎢4⎤⎥ and Ab = ⎡⎢ −1 2⎤⎥ ⎡⎢4⎤⎥ = ⎡⎢ 8 ⎤⎥ , ⎣6⎦ ⎣− 3 4⎦ ⎣6⎦ ⎣12⎦

so that r{Q} = r {[b Ab]} = r{⎡4 8 ⎤} = 1 ⎢⎣6 12⎥⎦ Since the rank of Q is less than 2, the system is uncontrollable. ________________________________________________________________________________________

The above example illustrates the importance of the observability concept. In this case we have a non-stable system, whose instability is not observed in the output measurement.

Definition The matrix R = [ C | C T A | L | C T A n−1 ]T is called the system observability matrix.

The Observability Criterion rank{R} = n ⇔ the system is observable

Example 6.3_______________________________________

Example 6.5_______________________________________

Determine whether the system, governed by the equation

Using observability criterion, verify that the system examined in example 6.4 is unobservable

⎡ x&1 ⎤ ⎡ 0 6 ⎤ ⎡ x1 ⎤ ⎡ 3 6 ⎤ ⎡ u1 ⎤ ⎢⎣ x& 2 ⎥⎦ = ⎢⎣− 1 − 5⎥⎦ ⎢⎣ x 2 ⎥⎦ + ⎢⎣− 1 − 2⎥⎦ ⎢⎣u 2 ⎥⎦ is controllable. Using the controllability criterion ⎡ ⎤ r{Q} = r {[B A B ]} = r{⎢ 3 6 − 6 − 12⎥} 4 ⎦ ⎣− 1 − 2 2 It is obvious that the rank of this matrix is 1. The system is therefore uncontrollable.

For this system c ′ = [3 − 2] and A = ⎡ −5 4⎤ , ⎢⎣− 6 5⎥⎦

hence R = ⎡ 3 −2⎤ , so that r {R} = 1 ⎢⎣− 3 2 ⎦⎥

________________________________________________________________________________________

6.3 Observability

⎧ z& = A1z + B1u . If a row of Consider the system in the form ⎨ ⎩y = C1x the matrix C1 is zero, the corresponding mode of the system will not appear in the output y . In this case the system is unobservable, since we cannot determine the state variable corresponding to the row of zeros in C1 from y . Example 6.4_______________________________________

Determine whether the system have the state equations ⎡ x&1 ⎤ ⎡ − 5 4⎤ ⎡ x1 ⎤ ⎡1 ⎤ ⎢⎣ x& 2 ⎥⎦ = ⎢⎣− 6 5⎥⎦ ⎢⎣ x 2 ⎥⎦ + ⎢⎣2⎥⎦ u ⎡x ⎤ y = [3 − 2]⎢ 1 ⎥ ⎣ x2 ⎦ is observable. Using the modal matrix P = ⎡1 2⎤ , the transform matrix ⎢⎣1 3⎥⎦ x = P z , transform the state-equations into ⎡ z&1 ⎤ ⎡− 1 0⎤ ⎡ z1 ⎤ ⎡− 1⎤ ⎢⎣ z& 2 ⎥⎦ = ⎢⎣ 0 1⎥⎦ ⎢⎣ z 2 ⎥⎦ + ⎢⎣ 1 ⎥⎦ u ⎡z ⎤ y = [1 0]⎢ 1 ⎥ ⎣z2 ⎦

This results shows that the system is unobservable because the output y does not influenced by the state variable z 2 .

Since the rank of R is less than 2, the system is unobservable. ________________________________________________________________________________________

6.4 Decomposition of the system state From the discussion in the two previous sections, the general state variables of a linear system can be classified into four exclusive groups as follows: • controllable and observable • controllable and unobservable • uncontrollable and observable • uncontrollable and unobservable Assuming that the system has distinct eigenvalues, by appropriate transforms the state equations can be reduced to the form below ⎡ x& 1 ⎤ ⎡ A1 0 0 0 ⎤ ⎡ x1 ⎤ ⎡ B1 ⎤ ⎢x& 2 ⎥ ⎢ 0 A2 0 0 ⎥ ⎢x 2 ⎥ ⎢ B ⎥ ⎢ x& ⎥ = ⎢ 0 0 A 0 ⎥ ⎢ x ⎥ + ⎢ 2 ⎥ u 3 ⎥⎢ 3⎥ ⎢ 0 ⎥ ⎢ 3⎥ ⎢ ⎣⎢x& 4 ⎦⎥ ⎣⎢ 0 0 0 A4 ⎥⎦ ⎣⎢x 4 ⎦⎥ ⎣ 0 ⎦ ⎡ x1 ⎤ ⎢x ⎥ y = [C1 0 C3 0] ⎢ 2 ⎥ x ⎢ 3⎥ ⎣⎢x 4 ⎦⎥

(6.11)

The (transformed) systems matrix A has been put into a “block-diagonal” form where each Ai (i = 1,2,3,4) is in diagonal form. The suffix i of the state-variable vector xi implies that the elements of this vector are the state variables corresponding to the i th group defined above.

________________________________________________________________________________________

___________________________________________________________________________________________________________ Chapter 6 Controllability and Observability

32

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Example 6.6______________________________________

Example 6.7_______________________________________

Consider the system

Illustrate the above discussion with the systems considered in examples 6.1 and 6.4

0 0 ⎤ ⎡ x1 ⎤ ⎡1 0 0 ⎥ ⎢⎢ x 2 ⎥⎥ ⎢0 0 0 ⎥ ⎢ x3 ⎥ + ⎢1 0 0 ⎥ ⎢ x 4 ⎥ ⎢0 0 0 ⎥ ⎢ x ⎥ ⎢1 5 0 − 4⎥⎦ ⎢ x ⎥ ⎢⎣0 ⎣ 6⎦ ⎡ x1 ⎤ ⎢ x2 ⎥ ⎡ y1 ⎤ ⎡0 0 − 1 2 0 1⎤ ⎢⎢ x3 ⎥⎥ ⎢⎣ y 2 ⎥⎦ = ⎢⎣1 0 1 − 1 0 1⎥⎦ ⎢ x4 ⎥ ⎢ x5 ⎥ ⎢⎣ x6 ⎥⎦ ⎡ x&1 ⎤ ⎡− 2 ⎢ x& 2 ⎥ ⎢ 0 ⎢ x& ⎥ ⎢ ⎢ 3⎥ = ⎢ 0 ⎢ x& 4 ⎥ ⎢ 0 ⎢ x&5 ⎥ ⎢ 0 ⎢⎣ x& 6 ⎥⎦ ⎣ 0

0 −1 0 0 0 0

0 0 0 0 1 0 0 −3 0 0 0 0

− 1⎤ 0⎥ 1 ⎥ ⎡ u1 ⎤ 0 ⎥ ⎢⎣u 2 ⎥⎦ 0⎥ 2 ⎥⎦

0 0 1 0 0 −4 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 −3 0 0

0 ⎤ ⎡ x1 ⎤ ⎡1 0 ⎥ ⎢⎢ x3 ⎥⎥ ⎢1 0 ⎥ ⎢ x 6 ⎥ + ⎢0 0 ⎥ ⎢ x5 ⎥ ⎢1 0 ⎥ ⎢ x ⎥ ⎢0 − 1⎥⎦ ⎢ x 4 ⎥ ⎢⎣0 ⎣ 2⎦

⎡ z&1 ⎤ ⎡1 0⎤ ⎡ z1 ⎤ ⎡0⎤ ⎢⎣ z& 2 ⎥⎦ = ⎢⎣0 2⎦⎥ ⎢⎣ z 2 ⎥⎦ + ⎢⎣2⎥⎦ u ⎡z ⎤ y = [− 1 − 4] ⎢ 1 ⎥ ⎣z2 ⎦

and can be represented by Fig.6.2

By inspection we can rewritten the above system in the following form ⎡ x&1 ⎤ ⎡− 2 ⎢ x&3 ⎥ ⎢ 0 ⎢ x& ⎥ ⎢ ⎢ 6⎥ = ⎢ 0 ⎢ x&5 ⎥ ⎢ 0 ⎢ x& 4 ⎥ ⎢ 0 ⎢⎣ x& 2 ⎥⎦ ⎣ 0

In example 6.1, the state equations were transferred into the diagonal form

− 1⎤ 1⎥ 2 ⎥ ⎡ u1 ⎤ 0 ⎥ ⎢⎣u 2 ⎥⎦ 0⎥ 0 ⎥⎦

u

z&1 = z1

-1

z&2 = 2 z2 + 2u

-4

y

Fig. 6.2 The system has two modes corresponding to the poles λ = 1 and λ = 2 . The transfer function

⎡ 1 ⎤ ⎢ s − 2 0 ⎥ ⎡0⎤ −8 G ( s) = C [ sI − A] B = [− 1 − 4] ⎢ = 1 ⎥ ⎢⎣2⎥⎦ s − 2 ⎢ 0 ⎥ s − 1⎦ ⎣ −1

⎡ x1 ⎤ ⎢ x3 ⎥ ⎡ y1 ⎤ ⎡0 − 1 1 0 2 0⎤ ⎢⎢ x6 ⎥⎥ ⎢⎣ y 2 ⎥⎦ = ⎢⎣1 1 1 0 − 1 0⎥⎦ ⎢ x5 ⎥ ⎢ x4 ⎥ ⎣⎢ x2 ⎥⎦

It is seen that the transfer function involves only the mode corresponding to λ = 2 - the controllable and observable one. The uncontrollable mode λ = 1 does not appear in the transfer function.

hence • controllable and observable

x1, x 3 , x 6

In example 6.4 the transformed state equations are

• controllable and unobservable

x5

• uncontrollable and observable

x4

• uncontrollable and unobservable

x2

⎡ z&1 ⎤ ⎡− 1 0⎤ ⎡ z1 ⎤ ⎡− 1⎤ ⎢⎣ z& 2 ⎥⎦ = ⎢⎣ 0 1⎥⎦ ⎢⎣ z 2 ⎥⎦ + ⎢⎣ 1 ⎥⎦ u ⎡z ⎤ y = [1 0] ⎢ 1 ⎥ ⎣z2 ⎦

________________________________________________________________________________________

We can represent the decompositions of the state variables into four groups by a diagram

and can be represented in Fig.6.3

S1

u

u

z&1 = z1 + u

y

y

z&2 = 2 z2 + 2u

S2

Fig. 6.3 S3 S4 Fig. 6.1 In general, a transfer function G ( s ) represents only the subsystem S1 of the system considered, and indeed on adding to S1 the subsystem S 2 , S 3 , S 4 will not effect G ( s ) .

In this case ⎡ 1 ⎤ ⎢ s + 1 0 ⎥ ⎡− 1⎤ −1 G ( s) = C [ sI − A] B = [1 0] ⎢ = 1 ⎥ ⎢⎣ 1 ⎥⎦ s + 1 ⎢ 0 ⎥ s − 1⎦ ⎣ −1

This result shows that the unobservable mode λ = 1 does not involved in the transfer function. ________________________________________________________________________________________

___________________________________________________________________________________________________________ Chapter 6 Controllability and Observability

33

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Definition ⎧ x& = A x + B u The state equation ⎨ are said to be a realization ⎩y = C x −1

of the transfer function G (s ) if G ( s ) = C [ sI − A] B .

Definition

⎡ 0 ⎢ ⎢ 0 C=⎢ M ⎢ ⎢ 0 ⎢− a ⎣ 0

1

0

0

1

M

M

0 − a1

0 − a2

⎤ ⎡0 ⎤ ⎥ ⎢ ⎥ 0 0 ⎥ ⎢0 ⎥ O M ⎥ and d = T −1b = ⎢0⎥ ⎥ ⎢ ⎥ L 1 ⎥ ⎢M⎥ ⎢1⎥ L − a n −1 ⎥⎦ ⎣ ⎦ (6.15)

L

0

⎧ x& = A x + B u A realization equation ⎨ of a transfer function ⎩y = C x G ( s) is said to be minimal if it is both controllable and observable.

Also the characteristic equation of A yields

Example 6.8_______________________________________

We construct the transformation matrix T in the following way:

Obtain a minimal realization of a system having the transfer function G ( s) =

1 ⎡ − s −( s + 1)⎤ ( s + 1)( s + 2) ⎢⎣3s + 2 2( s + 1) ⎥⎦

1 ⎡ − s −( s + 1)⎤ ( s + 1)( s + 2) ⎢⎣3s + 2 2( s + 1) ⎥⎦

⎡ 1 ⎤ 0 ⎥ ⎡ c11 c12 ⎤ ⎢ s + 1 ⎡ b11 b12 ⎤ =⎢ ⎢ 1 ⎥ ⎢⎣b21 b22 ⎥⎦ ⎣c 21 c 22 ⎥⎦ ⎢ 0 ⎥ s + 2⎦ ⎣

p ′q 2 ⎡ p ′q1 ⎢ ′ ′Aq 2 p A q p 1 ⎢ ⎢ M M ⎢ n−1 n −1 ′ ′ ⎣⎢p A q1 p A q 2

________________________________________________________________________________________

6.5 A Transformation into the Companion Form There exists a transformation (6.12)

transforming the system

is in a companion form, that is

L p ′q n ⎤ ⎡1 ⎥ ⎢ L p ′Aq n ⎥ ⎢0 = ⎥ ⎢M O M ⎥ ⎢ L p ′A n−1q n ⎦⎥ ⎣⎢0

0 L 0⎤ ⎥ 1 L 0⎥ M O M⎥ ⎥ 0 L 1 ⎦⎥ (6.18)

On taking the product ⎡ p′ ⎤ ⎢ ′ ⎥ pA ⎥ T −1 AT = ⎢ A [q1 q1 L q n ] ⎢ M ⎥ ⎢ n−1 ⎥ ⎢⎣p ′A ⎥⎦

we obtain the matrix (6.13)

where A is n × n matrix assumed to have n distinct eigenvalues, into

where C = T −1 AT

(6.17)

so that T −1T = I , hence

⎡ x&1 ⎤ ⎡− 1 0 ⎤ ⎡ x1 ⎤ ⎡1 0⎤ ⎢⎣ x& 2 ⎥⎦ = ⎢⎣ 0 − 2⎥⎦ ⎢⎣ x 2 ⎥⎦ + ⎢⎣2 1⎥⎦ u ⎡x ⎤ y = ⎡⎢ 1 − 1⎤⎥ ⎢ 1 ⎥ ⎣− 1 2 ⎦ ⎣ x 2 ⎦

z& = C z + d u

(6.16)

T = [q1 q1 L q n ]

And the minimal realization is therefore

x& = A x + b u

T −1

⎡ p′ ⎤ ⎢ ′ ⎥ pA ⎥ =⎢ ⎢ M ⎥ ⎢ n −1 ⎥ ⎢⎣p ′A ⎥⎦

Assume for the moment that T −1 is a non-singular matrix, so that T exists and has the portioned form

= C [ sI − A] −1 B

x =Tz

Let p be a vector of order n × 1 . Then p ′A, p ′A 2 , L , p ′A n −1 are n row vectors, which we collect to make up the n rows of the matrix T −1 , that is

The system has two modes λ = −1 and λ = −2 . Hence we can write G ( s) =

λ I − A = λn + a n−1λn−1 + L + a 0 = λ I − C = 0

(6.14)

p ′Aq 2 ⎡ p ′Aq1 ⎢ ′ 2 ′A 2 q 2 p A q p 1 ⎢ ⎢ M M ⎢ n n ⎣⎢p ′A q1 p ′A q 2

L p ′Aq n ⎤ ⎥ L p ′A 2 q n ⎥ O M ⎥ ⎥ L p ′A n q n ⎦⎥

(6.19)

Note that the first (n − 1) rows of equations (6.19) are the same (row by row) as the last (n − 1) rows of the matrix equation (6.18) ⇒ the matrix equation (6.19) is in the companion form, and

___________________________________________________________________________________________________________ Chapter 6 Controllability and Observability

34

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

⎧ − a0 = p ′A n q1 ⎪ ⎪ − a1 = p ′A n q 2 ⎨ ⎪ ⎪ n ⎩− a n−1 = p ′A q n

(6.20)

We must also have p ′b = 0 ⎧ ⎡0 ⎤ ⎡ p′ ⎤ ⎪ ⎢ ⎥ ⎢ ′ ⎥ ′ p Ab = 0 ⎢ p A ⎥ b = ⎢0⎥ or ⎪⎨ ⎢M⎥ ⎢ M ⎥ ⎪ ⎢ ⎥ ⎢ n−1 ⎥ ⎪p ′A n −1b = 0 ⎢⎣1 ⎥⎦ ⎢⎣p ′A ⎥⎦ ⎩

that is

[

]

Ab L A n −1b = d ′

p′ b

(6.21)

Eq. (6.21) is a system of n simultaneous equations in n unknowns (the components of p ) and will have a unique solution if and only if the rank of the system controllability matrix is n , that is if

[

rank b

]

Ab L A n −1b = n

(6.22)

(then the system is controllable). Example 6.9_______________________________________

The state equation of a system is ⎡ 13 35 ⎤ ⎡−2⎤ x& = ⎢ ⎥x + ⎢ ⎥u ⎣− 6 − 16⎦ ⎣1⎦

Find the transformation so that the state matrix of the transformed system is in a companion form. ⎧ p ′b = 0 , we obtain Using ⎨ ⎩p ′Ab = 0

[ p1

⎡−2⎤ p 2 ] ⎢ ⎥ = 0 and [ p1 ⎣1⎦

⎡ 13 35 ⎤ ⎡−2⎤ p2 ]⎢ ⎥⎢ ⎥ =1 ⎣− 6 − 16⎦ ⎣ 1 ⎦

that is ⎧−2 p1 + p 2 = 0 ⎨ ⎩9 p1 − 4 p 2 = 1

⇒

⎧ p1 = 1 ⎨ ⎩ p2 = 2

and the transformation matrix ⎡ 3 −2⎤ ⎡1 2⎤ −1 T =⎢ ⎥ and T = ⎢ ⎥ − 1 1 ⎣ ⎦ ⎣1 3⎦

The required transformation is ⎡ 3 −2⎤ x& = ⎢ ⎥z ⎣− 1 1 ⎦ ________________________________________________________________________________________

___________________________________________________________________________________________________________

Chapter 6 Controllability and Observability

35

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.7 Multivariable Feedback and Pole Location 7.1 Introduction

Example 7.1_______________________________________

Given the system x& = ⎡ −13 8 ⎤ x + ⎡1 ⎤ u , find the feedback ⎢⎣− 25 15⎥⎦ ⎢⎣2⎥⎦ control law so that the resulting closed loop system has poles at λ = −2 and λ = −3 .

7.2 State Feedback of a SISO System Consider a SISO system described by the state equations

x& = A x + b u

(7.1)

The dynamics of this system are determined by the poles, that is, the eigenvalues λ1 , λ 2 , L , λ n (assumed distinct) of the

With

the

transformation

matrix

T = ⎡⎢1 1 ⎤⎥ ⎣1 2⎦

and

matrix A .

T −1 = ⎡ 2 −1⎤ the state equation becomes ⎢⎣− 1 1 ⎥⎦

If system is controllable, there exists a matrix T such that x = T z transforms the system defined by equation (7.1) to

z& = T −1 A T z + T −1b u = ⎡ 0 1 ⎤ z + ⎡0⎤ u ⎢⎣1⎥⎦ ⎢⎣− 5 2⎥⎦

z& = C z + d u where

the characteristic equation in this case is λ12 − 2λ + 5 = 0 and

(7.2)

1 0 ⎡ 0 ⎢ 0 0 1 C = T −1 A T = ⎢ M M M ⎢ 0 0 0 ⎢⎣− a1 − a 2 − a3

L 0 ⎤ ⎡0 ⎤ ⎢0 ⎥ L 0 ⎥ O M ⎥ , d = T −1b = ⎢ M ⎥ ⎢0 ⎥ L 1 ⎥ ⎢⎣1 ⎥⎦ L − a n ⎥⎦

the poles are at λ = 1 ± 2 i so that the system is unstable. ⎡z ⎤ Apply feedback law in the form u = [k1 k 2 ] ⎢ 1 ⎥ and using ⎣z2 ⎦ equation (7.5) k1 = a1 − r1 = 5 − 6 = −1

Now apply the feedback control of the form

k 2 = a 2 − r2 = −2 − 5 = −7

u (t ) = k′z , k ′ = [k1 , k 2 , L , k n ]

(7.3) and the control law is

yields z& = [C + d k ′] z

(7.4)

u (t ) = k ′ z = k ′ T −1 x = [− 1 − 7 ] ⎡⎢ 2 −1⎤⎥ x = [5 − 6] x ⎣− 1 1 ⎦

which is the system whose dynamics are determined by the eigenvalues of [C + d k′] , that is by the eigenvalues of

• Checking the poles of the closed loop system After assigned poles

1 0 ⎡ 0 ⎢ 0 0 1 ⎢ M M M ⎢ 0 0 0 ⎢⎣− a1 − a 2 − a3

x& = Ax + bu = ( A + bk ′T −1 ) x

L 0 ⎤ ⎡0 0 0 L 0 ⎥ ⎢0 0 0 O M ⎥+⎢M M M L 1 ⎥ ⎢0 0 0 L − a n ⎥⎦ ⎢⎣k1 k 2 k 3

L L O L L

1 0 ⎡ 0 ⎢ 0 0 1 ⎢ M M M ⎢ 0 0 0 ⎢⎣− r1 − r2 − r3

0⎤ 0⎥ M⎥= 0⎥ k n ⎥⎦

and

L 0 ⎤ L 0 ⎥ O M ⎥ L 1 ⎥ L − rn ⎥⎦

(7.5) where − ri = k i − ai , (i = 1,2, L , n)

7.3 Multivariable Systems

By appropriate choice of the components k1 , k 2 , L , k n of k in equation (7.5) we can assign arbitrary values to the coefficients r1 , r2 , L , rn and so achieve any desired pole configurations. Since z = T −1 x it follows that x

λ 2 + 5λ + 6 = (λ + 2)(λ + 3) = 0

________________________________________________________________________________________

s n + rn s n −1 + K + r2 s + r1 = 0

u (t ) = k ′ z = k ′ T

The corresponding characteristic equation is

hence the closed loop poles were assigned at λ = −2 and λ = −3 as desired.

The characteristic equation of [C + d k ′] is

−1

⎡ −13 8 ⎤ ⎡1 ⎤ ⎡ −8 2⎤ A + bk ′T −1 = ⎢ ⎥ + ⎢ ⎥ [5 − 6] = ⎢ ⎥ − 25 15 2 ⎣ ⎦ ⎣ ⎦ ⎣− 15 3⎦

(7.6)

The pole assignment for a multivariable system is more complicated than for a scalar system since a multitude of inputs are to be controlled, and not just one ! We will discuss a relatively simple method dealing with this problem.

___________________________________________________________________________________________________________ 36 Chapter 7 Multivariable Feedback and Pole Location

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

• Matrix property If P and Q are matrices of orders m × n and n × m respectively, then | I m − PQ | = | I n − Q P |

(7.7)

Now we consider a multivariable system having n state variables and m inputs. To simplify the mathematical manipulations we shall initially assume that the state matrix is in a diagonal form, that is the state equation has the form z& = Λz + B1u

(7.10)

where the notation | A | stands for the determinant of A . I m and I n are the unit matrices of order m × n and n × m respectively.

• Proof Firstly, take a notice about some matrix properties

| sI − Λ | = ( s − λ1 )( s − λ 2 ) L ( s − λ n )

0

I

=

I

0

A I

(7.11)

where the system poles are assumed to be distinct. To correspond to equation (7.2) we assume a feedback law of the form

For any square matrix A I 0 A 0 = =| A | 0 A 0 I A

B1 ∈ Rn×m , z ∈ Rn×1 , u ∈ Rm×1

The system open loop characteristic polynomial is

For any A, B : | AB | = | BA |

I

where Λ = diag{λ1 , λ 2 , L , λ n }

u = K1 z

(7.12)

where K1 ∈ Rm×n is a proportional gain matrix.

=1

The closed loop system dynamics

When n = m

z& = (Λ + B1 K1 ) z

| I − PQ | = | P( P −1 − Q) | = | ( P −1 − Q) P | = | I − QP |

(7.13)

The objective of this procedure is to choose K1 so that the When n ≠ m Im P I = m Q In Q =

P In

closed loop system poles ρ1 , ρ 2 , L , ρ n have prescribed values. These poles are the solution of the closed loop characteristic equation, that is the roots of

Im 0 − Q In

I m − PQ P 0 In

(7.9)

= | I m − PQ |

Im Q

P I 0 = m In − Q In =

Im 0

| sI − Λ + B1 K1 | = ( s − ρ1 )( s − ρ 2 ) L ( s − ρ n ) = 0

(7.14)

Select K1 to have the so-called dyadic form Im Q

P In

P I n − QP

= | I n − QP | ⇒ | I m − PQ | = | I n − Q P |

K1 = f d ′

(7.8)

(7.15)

where f = [ f1

f2 L

f m ]′ and d ′ = [d1

d2 L dn ]

with this choice, (7.14) becomes (7.7)

| sI − Λ + B1f d ′ | = | sI − Λ | | I − ( sI − Λ ) −1 B1f d ′ |

Example 7.2_______________________________________

Given P = [ p1

Taking the determinant of both sides and using (7.11) and (7.14) we obtain

p 2 ] and

show that | I m − P Q | = | I n − Q P | P = [ p1

p2 ]

⎡q ⎤ and Q = ⎢ 1 ⎥ ⇒ m = 1, n = 2 ⎣q 2 ⎦

Hence | I m − PQ | = 1 − ( p1q1 + p 2 q 2 ) Also ⎡q p | I n − PQ | = I 2 − ⎢ 1 1 ⎣q 2 p1

q1 p 2 ⎤ ⎡1 − q1 p1 − q1 p 2 ⎤ ⎥=⎢ ⎥ q 2 p 2 ⎦ ⎣ − q 2 p1 1 − q 2 p 2 ⎦

Hence | I n − PQ | = (1 − q1 p1 )(1 − q 2 p 2 ) − q 2 p1 q1 p 2

= 1 − ( p1q1 + p 2 q 2 )

n

n

i =1

i =1

π ( s − ρ i ) = π ( s − λi ) | I − ( sI − λ ) −1 B1f d ′ |

Let

r2 L rn ]′

B1f = r = [r1 P = ( sI − Λ )

−1

(7.17)

(7.19)

B1f ∈ Rn×1

Q = d ′ ∈ R1×n ,

It follows that | sI − Λ + B1f d ′ | = | 1 − d ′( sI − Λ ) −1 B1f |

(7.18)

we then have ___________________________________________________________________________________________________________ 37 Chapter 7 Multivariable Feedback and Pole Location ________________________________________________________________________________________

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

d ′( sI − Λ ) −1 B1f

= [d1

=

⎡ 1 ⎢s − λ 1 ⎢ ⎢ 0 L d n ]⎢ ⎢ M ⎢ ⎢ 0 ⎣⎢

d2

0 1 s − λ2 M 0

⎤ ⎥ ⎥ ⎡ r1 ⎤ ⎥ ⎢r2 ⎥ L ⎥⎢ ⎥ ⎢M⎥ O M ⎥⎢ ⎥ 1 ⎥ ⎣⎢rn ⎦⎥ ⎥ L s − λ n ⎦⎥ L

∑s−λ

⎡−1 0 ⎤ ⎡2⎤ x& = ⎢ ⎥x + ⎢ ⎥u 0 − 2 ⎣ ⎦ ⎣1 ⎦

find a proportional gain matrix K1 such that the closed loop system poles are both at λ = −3 . In this case n = 2, m = 1 so that K1 = f d ′ is a one-row matrix, and f is one-element matrix

d i ri

i =1

Given a system defined by the state equation

0

d r d1r1 d r + 2 2 +L+ n n s − λ1 s − λ 2 s − λn n

=

Example 7.3_______________________________________

(7.20)

The open loop characteristic polynomial is

i

( s − λ1 )( s − λ 2 ) = ( s + 1)( s + 2) = s 2 + 3s + 2

Substituting Eqs. (7.18) and (7.20) into Eq. (7.17) we obtain ⎛ π ( s − ρ i ) = π ( s − λi )⎜⎜1 − i =1 i =1 ⎝ n

n

n

∑ i =1

The closed loop characteristic polynomial is

d i ri ⎞⎟ s − λi ⎟ ⎠

( s − ρ1 )( s − ρ 2 ) = ( s + 3) 2 = s 2 + 6 s + 9

The residuals are

Hence n

R1 = lim

n

π ( s − λi ) − π ( s − ρ i )

i =1

i =1

n

n

=

π ( s − λi )

∑s−λ d i ri

i =1

i =1

s →−1

(7.21) i

and R2 = lim

s →−2

On taking partial fraction n

i =1

i =1

n

n

=

π ( s − λi )

∑s−λ Ri

i =1

i =1

(7.22) i

where the residues Ri at the eigenvalue λ j ( j = 1,2, L , n) can be evaluated as ⎡ ⎤ ( s − λ j ) ⎢ π ( s − λi ) − π ( s − ρ i )⎥ i i = 1 = 1 ⎣ ⎦ n

R j = lim

s →λ j

( s + 2)(−3s − 7) =1 ( s + 1)( s + 2)

Hence (Eq. (7.24)): r1d1 = −4 and r2 d 2 = 1

n

π ( s − λi ) − π ( s − ρ i )

( s + 1)(−3s − 7) = −4 ( s + 1)( s + 2)

⎡ r1 ⎤ ⎡ 2⎤ ⎡2 f ⎤ ⎡ r1 ⎤ ⎢ ⎥ f = ⎢ ⎥ so that ⎢ ⎥ = ⎢ ⎥ r 1 ⎣ ⎦ ⎣ f ⎦ ⎣r2 ⎦ ⎣ 2⎦ ⎧r1 = 2 and Choose f = 1 , ⇒ ⎨ ⎩r2 = 1

n

n

From Eq. (7.19)

(7.23)

π ( s − λi )

i =1

The procedure for the pole assignment can be summarized as follows:

Check The closed loop characteristic polynomial is ⎡ s 0⎤ ⎡− 1 0 ⎤ ⎡2⎤ | sI − Λ − B1 K1 | = ⎢ ⎥−⎢ ⎥ − ⎢ ⎥ [− 2 1] ⎣0 s ⎦ ⎣ 0 − 2 ⎦ ⎣1 ⎦

⎡ s 0⎤ ⎡ − 1 0 ⎤ ⎡− 4 2⎤ = ⎢ ⎥ ⎥−⎢ ⎥−⎢ ⎣0 s ⎦ ⎣ 0 − 2 ⎦ ⎣ − 2 1 ⎦

• Calculate the residues R j ( j = 1,2, L , n) from Eq. (7.23) n ⎡n ⎤ ( s − λ j ) ⎢ π ( s − λi ) − π ( s − ρ i ) ⎥ i i = 1 = 1 ⎣ ⎦ R j = lim n s →λ j π ( s − λi )

=

r2 L rn ]′

• Calculate d i from the relation between (7.21) and (7.22)

• Calculate gain K1 from Eq. (7.15) K1 = f d ′

−2

2

s +1

________________________________________________________________________________________

• Calculate ri after choosing f i using Eq. (19)

d jrj = R j

s+5

= ( s + 3) 2

i =1

B1f = r = [r1

⎧d1 = −2 and K1 = [−2 1] . ⎨ ⎩d 2 = 1

(7.24)

In general case, if the state matrix is not necessarily in a diagonal form. Assume that x& = A x + B u

(7.25)

Using the transformation x = T z , Eq. (7.25) becomes

___________________________________________________________________________________________________________ 38 Chapter 7 Multivariable Feedback and Pole Location

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

z& = T −1 AT z + T −1 B u

= Λ z + B1u

(7.26)

And the closed loop system poles are the roots of the characteristic equation s + 59 98 = s 2 + 6s + 9 = 0 − 32 s − 53

where Λ = T −1 AT and B1 = T −1 B . Applying the inverse transformation, z = T −1x to the resulting state equation, we obtain

as required.

x& = T ΛT −1x + TB1 K1T −1x

7.4 Observers

= A x + BK x

(7.27)

where K = K1T −1

(7.28)

is the proportional gain matrix associated with the state of the system specified by Eq. (7.25). Example 7.4_______________________________________

Given ⎡ 4 10 ⎤ ⎡9 ⎤ x& = ⎢ ⎥x + ⎢ ⎥u − 3 − 7 ⎣ ⎦ ⎣ − 5⎦ Find the proportional gain matrix K such that the closed loop poles are both at λ = −3 . On applying the transformation 5⎤ ⎡2 x =Tz = ⎢ ⎥z − 1 − 3⎦ ⎣

for which 5 ⎤ ⎡3 T −1 = ⎢ ⎥ − − 1 2⎦ ⎣

________________________________________________________________________________________

In may control system it is not possible to measure entire state vector, although the measurement of some of the state variables is practical. Fortunately, Luenberger has shown that a satisfactory estimate of the state vector can be obtained by using the available system inputs and outputs to derive a second linear system known as an observer. The observer is a dynamic system whose state vector, say z (t ) is an approximation to a linear transformation of the state x(t ) , that is z (t ) ≈ T x(t )

(7.29)

The entire state vector z is of course, available for feed back. Consider the system x& = Ax + Bu y = c ′x

(7.30)

where A ∈ Rn×n , B ∈ Rn×m and c ′ ∈ R1×n . To simplify the mathematical manipulation, the system has been assumed to have a single output. But the results are applicable to multi-output systems. u

system x& = Ax + Bu

x

y output

It is found that the system is the same as the one considered in example 7.3. Hence K1 = [−2 1]

observer z& = Fz + hy + Nu

and 5 ⎤ ⎡3 K = K1T −1 = [− 2 1] ⎢ ⎥ = [− 7 − 12] − − 1 2⎦ ⎣

so the feedback control u = [−7 12] x

z

Fig. 7.1 If the dynamics of the observer when regarded as a free system (that is when the input is zero) are defined by z& = F z then, when the inputs are applied, using equation z& = F z + h y + N u

(7.31)

which achieve the desired pole assignment.

where h is an arbitrary vector (of order n × 1 ). Using (7.30) and (7.31), we form the difference

Check The closed loop system yields

z& - T x& = F z + h y + N u - TAx-TBu = F (z − Tx) + (hc ′ − TA + FT )x + ( N − TB )u

⎡ 4 10 ⎤ ⎡−63 −108⎤ ⎡−59 −98⎤ A + BK = ⎢ ⎥+⎢ ⎥=⎢ ⎥ 60 ⎦ ⎣ 32 53 ⎦ ⎣− 3 − 7 ⎦ ⎣ 35

If we choose hc ′ = TA − FT N = TB

(7.32)

(7.33)

___________________________________________________________________________________________________________ 39 Chapter 7 Multivariable Feedback and Pole Location

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

that is

Eq. (7.32) becomes

λ 2 + (h1 + 3)λ + (2h1 + h2 + 2) = 0

d ( z − Tx) = F ( z − Tx) dt

which has the solution

Since the observer eigenvalues are at λ = −3 , the above equation is identical to

z − Tx = e Ft [z (0) − T x(0)]

λ2 + 6λ + 9 = 0

or

hence

z = Tx + e Ft [z (0) − T x(0)]

(7.34)

h1 + 3 = 6

If it is possible to chose z (0) = T x(0) then z = T x for all

2h1 + h2 + 2 = 9

t ≥ 0 . Generally it is more realistic to assume that z (0) is only approximately equal to T x(0) , in which case it is

so that

required that e Ft approaches the zeros matrix as rapidly as possible. This can be achieved by choosing the eigenvalues of F to have enough large negative real parts.

⎡ h ⎤ ⎡3⎤ h = ⎢ 1⎥ = ⎢ ⎥ ⎣h2 ⎦ ⎣1⎦

Having obtained an estimate of T x , now we construct a

It follows that the state matrix of the required observer is

system having the transfer function T −1 to obtain x . This can be avoided if we choose T to be the identity transformation I . The resulting system is called an identity observer.

⎡ −4 1 ⎤ F =⎢ ⎥ ⎣ − 1 − 2⎦

In this case equations (7.33) becomes hc ′ = A − F N=B

(7.35)

From Eq. (7.35) we obtain

F = A − hc ′

(7.36)

so that the observer eigenvalues are seen to be determined by the choice of h .

Check On computing the difference z& − x& , we find ⎡ −4 1 ⎤ ⎡3 0⎤ ⎡1 −1⎤ z& − x& = ⎢ ⎥z + ⎢ ⎥x+ ⎢ ⎥u ⎣− 1 − 2⎦ ⎣1 0⎦ ⎣0 1 ⎦ ⎡− 1 1 ⎤ ⎡1 − 1⎤ −⎢ ⎥x− ⎢ ⎥u − 0 2 ⎣ ⎦ ⎣0 1 ⎦ ⎡− 4 1 ⎤ =⎢ ⎥ ( z − x) ⎣ − 1 − 2⎦ as required. ________________________________________________________________________________________

Example 7.7_______________________________________

Design an identity observer having both eigenvalues at λ = −3 for the system defined by the following state equations ⎡ −1 1 ⎤ ⎡1 −1⎤ x& = ⎢ ⎥x+ ⎢ ⎥u − 0 2 ⎣ ⎦ ⎣0 1 ⎦ y = [1 0] x y = [1 0] x ⇒ only the state variable x1 is available at the output.

F = A − hc ′ ⎡− 1 1 ⎤ ⎡ h1 ⎤ =⎢ ⎥ − ⎢ ⎥ [1 0] ⎣ 0 − 2 ⎦ ⎣ h2 ⎦ ⎡− 1 − h1 1 ⎤ =⎢ ⎥ − 2⎦ ⎣ − h2 The observer characteristic equation is | λI − F | = 0

___________________________________________________________________________________________________________ 40 Chapter 7 Multivariable Feedback and Pole Location

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.8 Introduction to Optimal Control Here v = (v1 , v 2 ) and the external force is the gravity force

8.1 Control and Optimal Control Take an example: the problem of rocket launching a satellite into an orbit about the earth. • A control problem would be that of choosing the thrust angle and rate of emission of the exhaust gases so that the rocket takes the satellite into its prescribed orbit. • An associated optimal control problem is to choose the controls to affect the transfer with, for example, minimum expenditure of fuel, or in minimum time. 8.2 Examples

only Fext = (0,−mg ) . The minimum fuel problem is to choose the controls, β and φ , so as to take the rocket from initial position to a prescribed height, say, y , in such a way as to minimize the fuel used. The fuel consumed is T

∫ β dt

(8.25)

o

8.2.1 Economic Growth 8.2.2 Resource Depletion 8.2.3 Exploited Population 8.2.4 Advertising Policies 8.2.5 Rocket Trajectories The governing equation of rocket motion is

dv = m F + Fext dt where m : rocket’s mass v : rocket’s velocity F : thrust produced by the rocket motor Fext : external force m

where T is the time at which y is reached.

8.2.6 Servo Problem A control surface on an aircraft is to be kept at rest at a position. Disturbances move the surface and if not corrected it would behave as a damped harmonic oscillator, for example (8.22)

θ&& + a θ& + w 2θ = 0

where θ is the angle from the desired position (that is, θ = 0 ). The disturbance gives initial values θ = θ 0 , θ& = θ 0′ , but a servomechanism applies a restoring torque, so that (8.26) is modified to

It can be seen that

θ&& + a θ& + w 2θ = u

c F = β m where : relative exhaust speed c β = − dm / dt : burning rate

(8.23)

(8.27)

The problem is to choose u (t ) , which will be bounded by u (t ) ≤ c , so as to bring the system to θ = 0, θ& = 0 in minimum time. We can write this time as

Let φ be the thrust attitude angle, that is the angle between the rocket axis and the horizontal, as in the Fig. 8.4, then the equations of motion are dv1 cβ = cos φ dt m dv 2 cβ = sin φ − g dt m

(8.26)

T

∫ 1.dt

(8.28)

0

and so this integral is required to be minimized.

8.3 Functionals (8.24)

v

All the above examples involve finding extremum values of integrals, subject to varying constraints. The integrals are all of the form

J=

path of rocket

∫

t1

t0

(8.29)

where F is a given function of function x(t ) , its derivative x& (t ) and the independent variable t . The path x(t ) is defined for t 0 ≤ t ≤ t1

r v

• given a path x = x1 (t ) ⇒ give the value J = J 1

φ

• given a path x = x 2 (t ) ⇒ give the value J = J 2 x

Fig. 8.4 Rocket Flight Path

F ( x, x& , t )dt

in general, J 1 ≠ J 2 , and we call integrals of the form (8.29) functional.

___________________________________________________________________________________________________________ Chapter 8 Introduction to Optimal Control

41

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

8.4 The Basic Optimal Control Problem

Consider the system of the general form x& = f (x, u, t ) , where f = ( f1 , f 2 , L , f n ) ′ and if the system is linear x& = A x + B u

(8.33)

The basic control is to choose the control vector u ∈ U such that the state vector x is transfer from x 0 to a terminal point at time T where some of the state variables are specified. The region U is called the admissible control region. If the transfer can be accomplished, the problem in optimal control is to effect the transfer so that the functional J=

T

∫f 0

0 ( x, u, t ) dt

(8.34)

is maximized (or minimized).

___________________________________________________________________________________________________________ Chapter 8 Introduction to Optimal Control

42

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.9 Variational Calculus 9.1 The Brachistochrone Problem Consider the Johann Bernoulli problem: Find the trajectory from point A to point B in a vertical plane such that a body slides down the trajectory (under constant gravity and no friction) in minimum time. y b A

Suppose we are looking for the path y = y (x) between A and B (Fig.9.3) which gives an extremum value for the functional

∫

t1

J = F ( y, dy / dx, x)dx

(9.5)

t0

y

r v

r g

B

β Cε

εη

a

B ( a , b)

A

α

x

Fig.9.1 The Brachistochrone Problem

a

Introduce coordinates as in Fig.9.1, so that A is the point A(0,0) and B is B (a, b) . Assuming that the particle is released from rest at A , conservation of energy gives 1 2

m v2 − m g x = 0

C0

(9.1)

The particle speed is v = x& 2 + y& 2 . From Fig.9.2, we can see that an element of arc length, δ s , is given approximately by δ s ≈ δ x 2 + δ y 2

b Fig.9.3 Possible paths between A and B

Let C 0 be the required curve, say y 0 ( x) , and let Cε be a neighboring curves defined by yε ( x) = y 0 ( x) + ε η ( x) where, ε : small parameter η (x) : arbitrary differential function of x , η (a) = η (b) = 0 .

gives the optimal curve C 0 .

δy

The value of J for the path Cε for the general problem is given by

δx Fig.9.2 Arc Length Element

J=

Hence, 2

v = lim

From (9.1) and (9.2), dt / ds = (2 gx) − the time of descent is given by t=

b

∫ F(y a

2

∫ (2 gx) ds

1 2

=

1 1

(2 g ) 2

∫

a

0

(9.2)

1

2

so that on integrating,

1 + (dy / dx) 2 dx x

0

+ ε η , y 0′ + ε η ′, x)dx ,

y′ =

dy dη , η′ = dx dx

For a given function y (x) , J = J (ε ) . For extremum values of J , we must have

dJ =0 dε ε =0

(9.7)

Taking the differentiation yields

(9.3)

Our problem is to find the path y = y (x) such that y (0) = 0, y (a ) = b and which minimizes (9.3). 9.2 Euler Equation

(9.6)

Thus equation (9.6) generates a whole class of neighboring curve Cε , depending on the value of ε . The value ε = 0

δs

δ s ds ⎛ dy ⎞ ⎛ dx ⎞ = = ⎜ ⎟ +⎜ ⎟ δ t →0 δ t dt ⎝ dt ⎠ ⎝ dt ⎠

x

dJ = dε

∫ [η F ( y b

a

y

0

+ ε η , y 0′ + ε η ′, x) +

]

y ′F y′ ( y 0 + ε η , y 0′ + ε η ′, x ) dx where Fy =

dF dF , F y′ = dy dy ′

___________________________________________________________________________________________________________ Chapter 9 Variational Calculus

43

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

dF dF dy dF dy ′ = + = η F y + η ′F y′ dε dy dε dy ′ dε

Example 9.1_______________________________________ Find the curve y (x) which gives an extremum value to the functional

y = y0 + εη y ′ = y 0′ + ε η ′

1

∫

J = ( y ′ 2 + 1)dx

Applying condition (9.7) we obtain

0

∫ [η F ( y , y′ , x) +η′F b

a

y

0

0

y ′ ( y0 , y0′ , x )

]dx = 0

with y (0) = 1, y (1) = 2

From now on all the partial derivations of F are evaluated on the optimal path C 0 , so we will leave off the ( y 0 , y 0′ , x) dependence. With this notation, on the optimal path, b

∫ (η F a

y

and on integrating, F y′ = A constant, that is

b

b a

F y′ = 2 y ′ = A

d ⎡ ⎤ + ⎢η F y − η ( F y′ )⎥ dx = 0 dx ⎣ ⎦

∫ a

Integrating again,

or b

η (b) F y′

x =b

(9.10) becomes d ( F y′ ) = 0 dx

+ η ′F y′ )dx = 0

and the integrating the second term by parts gives

η F y′

Here, F = y ′ 2 + 1 ⇒ F y = 0 . Hence the Euler equation

− η ( a ) F y′

x =a

d ⎡ ⎤ + ⎢η F y − η ( F y′ )⎥ dx = 0 dx ⎣ ⎦

∫ a

But η (a ) = η (b) = 0 , so we are left with

(9.8)

y = Ax / 2 + B From boundary conditions A = 2, B = 1 , and so the optimal curve is the straight line y = x +1

b

∫ a

d ⎡ ⎤ ⎢η F y − η dx ( F y′ )⎥ dx = 0 ⎣ ⎦

(9.9)

as illustrated in Fig. 9.5

y

This must hold for all differentiable functions η (x) such that

2

η (a) = η (b) = 0 We now require the following Lemma.

1

Lemma Let η (x) be a differentiable function for a ≤ x ≤ b such that η (a) = η (b) = 0 . If f is continuous on a ≤ x ≤ b and 0 1 Fig. 9.5 Optimal Path

b

∫ f ( x)η ( x) dx = 0

The corresponding extremum value of J is

a

1

for all such η (x) , then f = 0 on a ≤ x ≤ b .

d ( F y′ ) = 0 dx

∫

1

∫

J = ( y ′ 2 + 1)dx = 2 dx = 2

Applying the result of this Lemma to (9.9), immediately gives Fy −

x

0

0

_________________________________________________________________________________________

Example 9.2_______________________________________ (9.10)

on the optimal path. This is known as Euler’s equation (or the Euler-Lagrange equation) and must be satisfied by paths y (x ) which yield extremum values of the functional J .

Find the curve y (x) which has minimum length between (0,0) and (1,1) Let y = y ( x), 0 ≤ x ≤ 1 be any curve between (0,0) and (1,1) . The length, L , is given by

___________________________________________________________________________________________________________ Chapter 9 Variational Calculus

44

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

This is a question of a cycloid, the constant c being chosen so that the curve passes through the point (a, b) . A typical example of such a curve is shown in Fig. 9.6

1

L=

∫

1 + y ′ 2 dx

0

The integrand F = 1 + y ′ 2 is independent of y , that is,

x

A(0,0)

F y = 0 . Hence (9.10) becomes d ( F y′ ) = 0 dx

B ( a , b)

and on the integrating F y′ = A , constant. Thus y′ 1 + y′

y

Fig. 9.6 (a,b) Cycloid from A and B

=A 2

which implies that y ′ = B = constant. Hence y = B x + C , where C is a constant. To pass through (0,0) and (1,1), we obtain y = x . ________________________________________________________________________________________

Example 9.3 The Brachristochrone problem_____________

Determine the minimum value of the functional (9.3), that is a

1 + y′2 dx x

2g ∫ 1

t=

0

1 + y′2 x d ⇒ ( F y′ ) = F y = 0 dx y′ = c , constant. ⇒ F y′ = y (1 + y ′ 2 )

Let define x ≡

=0

d d Using (9.10): F y − ( F y′ ) = 0 ⇒ F y = ( F y′ ) , give dx dx dF d d ( y ′F y′ ) , therefore = y ′ ( F y′ ) + y ′′F y = dx dx dx

dx

1

(1 − cos θ ) ⇒ dx / dθ = sin θ /( 2c 2 ) ,

therefore

y=

1 2 c2

y=

1 2 c2 1 2c2

A plane curve, y = y ( x ) passes through the points (0, b) and ( a, c), (b ≠ c ) in the x − y plane ( a, b, c are positive constants) and the curve lies above the x-axis as shown in Fig. 9.7. y c

(θ − sin θ ) + A , A is a constant

Initial condition: ( x, y ) = (0,0) ≡ θ = 0 ⇒ A = 0 . Hence x=

(9.13)

Example 9.4 Minimum Surface Area___________________

x − c2 x2

2c 2

(9.12)

F − y ′F y′ = constant

1− c2x x

∫

1) F independent of y d d From (9.10): F y − ( F y′ ) = 0 ⇒ ( F y′ ) = 0 , hence dx dx

dF ∂F dy ∂F dy ′ ∂F = + + = y ′F y + y ′′F y ∂y dx ∂y ′ dx { dx ∂x

c2x

y=c

and

The Euler equation, (9.9), takes special forms when the integrand F is independent of either y or x :

2) F independent of x In general

F=

y′ =

________________________________________________________________________________________

F y′ = constant

such that y (0) = 0 and y (a ) = b .

⇒

⊗ Note that: Our theory has only shown us that we have a stationary value of t on this path. The path found does in fact corresponding to minimum t .

b

(a, c )

y = y (x )

(1 − xosθ )

(9.11) (θ − sin θ )

a Fig. 9.7 Curve rotates about x-axis

x

___________________________________________________________________________________________________________ Chapter 9 Variational Calculus

45

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Determine the curve y ( x) which, when rotated about the xaxis, yields a surface of revolution with minimum surface area.

y

B C0

The area of the surface of revolution is given by A=

∫

2πy ds = 2π

∫

a

0

A

2

⎛ dy ⎞ y 1 + ⎜ ⎟ dx ⎝ dx ⎠

so we wish to find the curve y ( x) to minimizes this integral. b a Fig. 9.8 Optimal and Neighboring Curves

The integrand F = y 1 + y ′ 2 . This is independent of x so we can apply (9.13)

Fy −

F − y ′F y′ = A , constant y 1 + y′2 −

⇒

y′2 y 1 + y′

Eq. (9.6) becomes η (b) F y′

=A

⇔

2

2

∫

⇔

(9.17)

x =b

− η ( a ) F y′

x=a

= 0 . This result

must be true for all differentiable functions η (x) . In fact for the class of all differentiable functions η (x) with the restriction

2

y = A (1 + y ′ ) dy 1 y′ = = y 2 − A2 dx A dy 1 = dx A 2 2 y −A

⇔

2

d ( F y′ ) = 0 dx

x

η (a) = 0

∫

which means that

⇔ cosh −1 ( y / A) = x / A + B , B: constant e x + e−x Hyperbolic cosine of x = cosh x ≡ 2

η (b) F y′

x =b

=0

Since the value of η (b) is arbitrary, we conclude that

so that y = A cosh( x / A + B)

(9.14)

The constants are chosen so that the curve passes through (0,b) and (a,c).

∂F = 0 at x = b ∂y ′

(9.18)

Similarly

_________________________________________________________________________________________

∂F = 0 at x = a ∂y ′

9.3 Free End Conditions

This section generalizes the result obtained in section 9.2 for finding extremum values of functionals. Find the path y = y ( x), a ≤ x ≤ b which gives extremum values to the functional J=

∫

b

F ( y , y ′, x)dx

(9.15)

(9.19)

and if only free at one end point, then it can be shown that ∂F = 0 at that end point. We can summarize these results by ∂y ′ ∂F = 0 at end point at which y is not specified. ∂y ′

(9.20)

a

Example 9.5 ______________________________________

with no restriction on the value of y at either end point.

Find the curve y = y (x) which minimizes the functional

The analysis used in 9.2 follows through directly, except that we no longer have the restriction η (a ) = η (b) = 0 . The optimal curve, C∈ , are illustrated in Fig.9.8.

J =

b

x =b

− η ( a ) F y′

2

+ 1)dx

0

where y (0) = 1 and y is free at x = 1 .

Form Eq. (9.8)

η (b) F y′

1

∫ ( y′

x=a

d ⎛ ⎞ ( F y′ ) ⎟ dx = 0 + η ⎜ Fy − dx ⎝ ⎠

∫

We have already seen, from the example 9.1, that the Euler equation gives

a

x (9.16) y= A +B The free end point also satisfy the Euler equation 2 ___________________________________________________________________________________________________________

Chapter 9 Variational Calculus

46

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Hence x 2 + &x&1 = 0 and x1 + &x&2 = 0 . Eliminating x 2 gives

Apply Eq. (9.20) to give ∂F = 0 at x = 1 ∂y ′

that is, 2 y ′ = 0 at x = 1 .

d 4 x1

⎧ y (0) = 1 ⎧ A = 0 ⇒⎨ The boundary conditions ⎨ . Hence the ⎩ y ′(1) = 0 ⎩ B = 1 optimal path is y ( x) = 1, 0 ≤ x ≤ 1

and the optimum value of J is J = 1 , which is clearly a minimum. _________________________________________________________________________________________

9.4 Constrains

The results obtained in the last two sections can be further generalized - firstly to the case where the integrand of more than one independent variable, - secondly to the case of optimization with constraints.

a

1

2

n

1

2

n

(9.21)

where, x1 , x 2 , L , x n are independent functions of t x&1 , x& 2 , L , x& n are differentiation with respect to t We obtain n Euler equations to be satisfied by the optimum path, that is ∂F d ⎛ ∂F − ⎜ ∂xi dt ⎜⎝ ∂x& i

⎞ ⎟ = 0 (i = 1,2, L , n) ⎟ ⎠

(9.22)

∫

0

0 ≤ t ≤π /2

_________________________________________________________________________________________

In many optimization problems, we have constraints on the relevant variables. We first deal with extremum values of a function, say f ( x1 , x 2 ) subject to a constraint of the form

multiplier, λ , and form the augmented function f * = f ( x1 , x 2 ) + λ g ( x1 , x 2 )

(9.24)

We now look for extremum values of f * , that is ∂f * ∂f * = =0 ∂x1 ∂x 2

(9.25)

Solving these equations, x1 , x 2 will be functions of the

f * ≡ f , we have the condition for the extremum values for f .

(9.23) Example 9.7 ______________________________________

Find the extremals of the functional J =

1 − e −π

λ is found by satisfying g ( x1 (λ ), x 2 (λ )) = 0 . Since

Example 9.6 ______________________________________

π /2

e t − e −t

parameter λ , that is, x1 (λ ), x 2 (λ ) . And appropriate value of

Also at any end point where xi is free ∂F =0 ∂x& i

x1 = x 2 = e −π / 2

one of x1 or x 2 using the constraint, we introduce a Lagrange

b

∫ F ( x , x ,L, x , x& , x& ,L, x& , t )dt

⎧ e −π / 2 ⎪A = 1 − e −π ⎪ ⎪ Applying the boundary conditions gives ⎨ B = − A , ⎪C = 0 ⎪ ⎪⎩ D = 0 Hence the optimum path is

g ( x1 , x 2 ) = 0 . Assuming that it is not possible to eliminate

If we wish to find extremum values of the functional J =

dt 4

⎧⎪ x1 = Ae t + Be −t + C sin t + D cos t − x1 = 0 ⇒ ⎨ ⎪⎩ x 2 = Ae t + Be −t − C sin t − D cos t

( x&12 + x& 22 + 2 x1 x 2 )dt

subject to boundary conditions x1 (0) = 0 , x1 (π / 2) = 1 , x 2 (π / 2) = 1 .

Using Eq. (9.22) to give ⎧ ∂F d ⎛ ∂F ⎞ d ⎛ ⎟=0 − ⎜⎜ ⎪ ⎜ 2 x 2 + (2 x&1 ) = 0 ⎟ ⎪ ∂x1 dt ⎝ ∂x&1 ⎠ dt ⎜ ⇒ ⎨ ⎜ d d ⎛ ∂F ⎞ ⎪ ∂F ⎜ ⎟ ⎜ 2 x1 + (2 x& 2 ) = 0 0 − = ⎪ ∂x ⎜ ⎟ dt ⎝ ⎩ 2 dt ⎝ ∂x& 2 ⎠

Find the extremum value of f ( x1 , x 2 ) = x12 + x 22

subject to x1 + 3 x 2 = 4 . We first note that we can easily eliminate x1 from the constraint equation to give f = (4 − 3 x 2 ) 2 + x 22

and we can then solve the problem by evaluating df / dx 2 = 0 . Using a Lagrange multiplier, we consider the augmented function

f * = x12 + x 22 + λ ( x1 + 3 x 2 − 4) Hence ___________________________________________________________________________________________________________

Chapter 9 Variational Calculus

47

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

For extremum value of f * ,

y: λ −

⎧ ∂f * ⎪ ∂x = 0 ⎧ 2 x1 + λ = 0 ⎪ 1 ⇒⎨ ⎨ ∂ * f ⎩2 x 2 + 3λ = 0 ⎪ =0 ⎪⎩ ∂x 2

d ⎛⎜ z& ⎜ dx 2 2 g x (1 + y& z& ) ⎝

z: −λ −

Hence x1 = λ / 2, x 2 = −3λ / 2 and to satisfy the constraint

y& d ⎛⎜ dx ⎜ 2 2 g x (1 + y& z& ) ⎝

⎞ ⎟=0 ⎟ ⎠

Eliminating λ , we obtain y& + z& d ⎛⎜ dx ⎜ 2 2 g x(1 + y& z& ) ⎝

(−λ / 2) + [3(−3λ / 2)] = 4 ⇒ λ = −4 / 5

⎞ ⎟=0 ⎟ ⎠

⎞ ⎟=0 ⎟ ⎠

and the extremum value of f occurs at x1 = 2 / 5,

y& + z&

⇒

x2 = 6 / 5

x(1 + y& z& )

_________________________________________________________________________________________

We can now return to the main problem: Suppose we wish for extremum values of (9.21) J=

b

∫ F ( x , x ,L, x , x& , x& ,L, x& , t )dt = ∫ F (x, x& , t )dt 1

a

2

n

1

2

n

a

(9.21) subject to a constraint between x1 , x 2 , L , x n , x&1 , x& 2 , L , x& n , t which we write as g (x, x& , t ) = 0

(9.26)

then we introduce a Lagrange multiplier, λ , and form the augmented functional

∫

using the constraint z = y − 1 , we obtain for the optimum path y&

b

=A

=

x(1 + y )

&2

A 2

_________________________________________________________________________________________

Example 9.9 ______________________________________

Minimize the cost functional J=

1 2

∫

2

0

&x&2 dt

where x(0) = 1, x& (0) = 1; x(2) = 0, x& (2) = 0

b

J * = ( F + λ g )dt

(9.27)

a

We now find the extremum values of J * in the usual way, regarding the variables as independent. The Lagrange multiplier can then be eliminated to yield the required optimum path.

To use Lagrange multiplier, we introduce the state variable x1 = x, x 2 = x& . The functional is now J=

1 2

∫

2

0

x& 2 dt

and we have the constraint Example 9.8 ______________________________________

Minimize the cost functional

J=

1 2g

∫

a

0

1 + y&z&

The augmented functional is dx

x

where the variables y = y ( x), z = z ( x) are subject to the constraint y = z +1

Introducing the Lagrange multiplier λ , the augmented functional is b⎡

1

a

2g

∫ ⎢⎣⎢

J* =

2⎛ 1

∫ ⎜⎝ 2 x& 0

2 2

⎞ + λ ( x 2 − x&1 ) ⎟ dt ⎠

Euler equation yields d ( −λ ) = 0 dt d x 2 : λ − ( x& ) = 0 dt x1 : 0 −

and where y (0) = 0, y (a ) = b

J* =

x 2 − x&1 = 0

⎤ 1 + y& z& + λ ( y − z − 1)⎥ dx x ⎦⎥

Eliminating λ , we obtain x1 = A

⇒

d 3 x2 dt 3

λ& = 0 λ − &x&2 = 0

= 0 , and

t3 t2 +B +Ct + D 3 2

x2 = A t 2 + B t + C Euler equation yields ___________________________________________________________________________________________________________

Chapter 9 Variational Calculus

48

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.4

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Applying the end conditions implies

1 3 7 2 t + t +t +1 2 4 3 2 7 x2 = t + t + 1 2 2 x1 =

y& + z&

⇒

x(1 + y& z& )

=A

and we obtain for the optimum path J=

1 2

2

∫ (3t − 7 / 2) dt = 13 / 4 2

0

_________________________________________________________________________________________

___________________________________________________________________________________________________________ Chapter 9 Variational Calculus

49

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.10 Optimal Control with Unbounded Continuous Controls and the corresponding optimal control is recovered from (10.2), giving

10.1 Introduction The methods appropriate for variational calculus can readily extended to optimal control problems, where the control vector u is unbounded. Technique for dealing with the general case u ∈ U , a bounded admissible control region, will be discussed in the next chapter. We start with a simple example. Suppose we wish to find extremum values of the functional J=

∫

T

2

(10.1)

sinh 2T

10.2 The Hamiltonian

Consider the general problem of minimizing J=

( x + u )dt 2

sinh 2 (T − t ) − 2 cosh 2 (T − t )

u = x0

T

∫f o

(10.6)

0 ( x, u , t ) dt

o

subject to the differential equation x : the state variable u : the control variable

where

x& = f ( x, u , t )

satisfy the differential equation x& + x = u

(10.2)

(10.7)

and boundary condition x(0) = a . With Lagrange multiplier, λ , we form the augmented functional T

∫ {f

0 ( x, u , t ) + λ [ f

with the boundary conditions x(0) = x 0 , x(T ) = 0 .

J* =

A direct method of solving this problem would be to eliminate u from (10.2) and (10.1), which reduces it to a straightforward calculus of variations problem. In general, it will not able to eliminate the control, so we will develop another method based on the Lagrange multiplier technique for dealing with constraints.

The integrand F = f 0 ( x, u , t ) + λ [ f ( x, u , t ) − x& ] is a function of the two variables x and u, and so we have two Euler equations

We introduce the Lagrange multiplier, λ , and from the augmented functional

(10.8-9)

J* =

∫

T

0

[ x 2 + u 2 + λ ( x& + x − u )] dt

o

(10.3)

The integrand F = x 2 + u 2 + λ ( x& + x − u ) is a function of two variables x and u. We can evaluate Euler’s equation for x and u, that is

∂F d ⎛ ∂F ⎞ − ⎜ ⎟=0 ∂x dt ⎝ ∂x& ⎠

⇒

2 x + λ − λ& = 0

(10.4)

∂F d ⎛ ∂F ⎞ − ⎜ ⎟=0 ∂u dt ⎝ ∂u& ⎠

⇒

2u − λ = 0

(10.5)

( x, u , t ) − x& ]} dt

⎧ ∂F d ⎛ ∂F ⎞ ⎧ ∂f 0 ∂f & − ⎜ ⎟=0 ⎪ ⎪⎪ ∂x + λ ∂x + λ = 0 ⎪ ∂x dt ⎝ ∂x& ⎠ ⇒ ⎨ ⎨ ∂f 0 ∂f ⎪ ⎪ ∂F − d ⎛⎜ ∂F ⎞⎟ = 0 =0 +λ ⎪⎩ ∂u dt ⎝ ∂u& ⎠ ∂u ∂u ⎩⎪ Now defining the Hamintonian

H = f0 + λ f

(10.10)

(10.8) and (10.9) can be rewritten as ∂H ∂x ∂H 0= ∂u

λ& = −

(10.11) (10.12)

From (10.2), (10.4) and (10.5), we obtain

These are the equations, together with (10.2), which govern the optimal paths.

&x& − 2 x = 0

Example 10.1______________________________________

⇒

x = Ae

2t

+ Be

− 2t

or, alternatively x = a cosh 2T + b sinh 2T .With boundary conditions, we obtain

x=

x 0 sinh 2 (T − t )

Using Hamintonian method, find the optimal control u which yields a stationary value for the functional the

1

∫

J = ( x 2 + u 2 )dt

where x& = u

0

with x(0) = 1 and not fixed at t = 1 . sinh 2T ___________________________________________________________________________________________________________ 50 Chapter 10 Optimal Control with Unbounded Continuous Controls

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

with α = 4 , (10.20) has solutions of the form x1 = e mt , where m satisfies

Here f 0 = x 2 + u 2 , f = u and the Hamintonian is H = x2 + u2 + λ u

2

m4 − m2 + ∂H (10.11): λ& = − ⇒ ∂x ∂H (10.12): 0 = ⇒ ∂u

λ& = −2 x

with the boundary conditions, we get 2u + λ = 0

x1 = [a + (b + a / 2 ) t ] e −t /

From the boundary conditions

_________________________________________________________________________________________

The method described in the last section can be readily extended to higher order systems. For example, we wish to minimize

∫

J=

T

∫f 0

10.3 Extension to Higher Order Systems

0 (x, u, t ) dt

2

x& i = f i (x, u, t )

(i = 1,2, L , n)

(10.28)

x 2 L x n ]T and u = [u1 u 2 L u n ]T

(10.13)

0

(10.27)

subject to the constraints

where, x = [x1

( x + α u ) dt 2

(10.24)

The basic problem is to find the optimal controls u which yield extremal values of

cosh(1 − t ) sinh(1 − t ) and u = − cosh 1 cosh 1

∞

(10.23)

−t / 2

u = [− a / 2 − ( 2 − 1)b − (b + a / 2 )(1 − 1 / 2 ) t / 2 ] e −t / 2 (10.25) General Problem

x = a sinh t + b cosh t

1 J= 2

2

x 2 = [b − (b + a / 2 ) t / 2 ] e

Hence we obtain &x& − x = 0 with the solution is

x=

1 ⎛ 2 1⎞ = ⎜m − ⎟ = 0 2 ⎝ 2⎠

We introduce the Lagrange multipliers, pi (i = 1,2, L , n) usually called the adjoint variables and form the augmented functional

where, &x& = − x& + u , α is constant

(10.14)

and initially x = a , x& = b , and as t → ∞ , both x and x& → 0 . Following the usual technique, we define x1 ≡ x, x 2 ≡ x& so that we have two constraints x 2 − x&1 = 0 − x 2 + u − x& 2 = 0

J* =

T

∫

0

⎡ ⎢ f0 + ⎢ ⎣

∑p (f i

i

i =1

⎤ − x& i )⎥ dt ⎥ ⎦

We also define the Hamintonian, H, is

(10.15) (10.16)

n

n

H = f0 +

∑p f

(10.29)

i i

i =1

The augmented function is J* =

∫

so that

∞ ⎡1

0

⎤ 2 2 ⎢ 2 ( x1 + α u ) + λ1 ( x 2 − x&1 ) + λ 2 (− x 2 + u − x& 2 )⎥ dt ⎣ ⎦

J* =

T

∫

0

⎡ ⎢H − ⎢ ⎣

⎤

n

∑ p x& ⎥⎥ dt i i

i =1

⎦

We now have three Euler equations, namely,

n

The integrand, F = H −

∂F d ⎛ ∂F ⎞ ⎟=0 − ⎜ ∂x1 dt ⎜⎝ ∂x&1 ⎟⎠ ∂F d ⎛ ∂F − ⎜ ∂x 2 dt ⎜⎝ ∂x& 2

⇒

⎞ ⎟=0 ⇒ ⎟ ⎠

∂F d ⎛ ∂F ⎞ − ⎜ ⎟=0 ∂u dt ⎝ ∂u& ⎠

⇒

x1 + λ&1 = 0

(10.17)

λ1 − λ 2 + λ&2 = 0

(10.18)

α u + λ2 = 0

(10.19)

dt

4

−

d 2 x1 dt 2

+

x1

α

i i

depends on x, u, t, and we

i =1

form (n+m) Euler equations, namely

∂F d ⎛ ∂F − ⎜ ∂xi dt ⎜⎝ ∂x& i

⎞ ∂H ⎟ = 0 (i = 1,2, L , n) that is, p& i = − ⎟ ∂xi ⎠ (10.30)

which are known as the adjoint equations; and

Eliminating λi , we obtain d 4 x1

∑ p x&

=0

d ⎛ ∂F ∂F − ⎜ ∂u j dt ⎜⎝ ∂u& j

(10.20)

⎞ ⎟ = 0 ( j = 1,2, L , m) that is, ∂H = 0 ⎟ ∂u j ⎠ (10.31)

___________________________________________________________________________________________________________ 51 Chapter 10 Optimal Control with Unbounded Continuous Controls

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

_________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

The optimal solutions are determined from (10.28,30,31). Using boundary conditions given in section 8.4, we have xi (0) given (i = 1,2, L , n) and xl (T ), (l = 1,2, L , q ) are given. The remaining values x q +1 (T ), L , x n (T ) are free, and so we

Example 10.3______________________________________

Find the optimal control u which gives an extremum value of the functional 1

must apply the free end condition (9.23):

∫

J = u 2 dt

∂F = 0 (k = q + 1, L , n) at t = T ∂x& i

where x&1 = x 2 , x& 2 = u

0

with x1 (0) = 1 , x1 (1) = 0 , x 2 (0) = 1 and x 2 (1) is not specified.

This gives

The Hamintonian is

p k (T ) = 0 (k = q + 1, L , n)

(10.32)

H = u 2 + p1 x 2 + p 2 u

which is known as the transversality condition. In general pi (T ) = 0 for any xi not specified at t = T .

Example 10.2______________________________________

(10.31): (10.30):

∂H = 0 = 2u + p 2 ⇒ ∂u p& 1 = 0 ⇒ p& 2 = − p1 = − A

1 p2 2 p1 = A

u=−

p 2 = − At + B

Find the control u which minimizes

Hence

1

∫

J = u dt 2

where x& = u + ax , a is constant, and

0

(i) x(0) = 1 (ii) x(0) = 1 , x(1) = 0 The Hamintonian is H = u 2 + p (u + ax) where p is the adjoint variable. (10.31): (10.30):

2u + p = 0 p& = − ap

u=

A B t− 2 2

A B x& 2 = t − 2 2

⇒

A 2 ⎧ ⎪⎪ x 2 = 4 t − ⎨ ⎪x = A t 3 − ⎪⎩ 1 12

B t +C 2 B 2 t + Ct + D 4

Boundary conditions give A = 12, B = 12, C = 1, D = 1 . The optimal control is u = 6(t − 1), 0 ≤ t ≤ 1 _________________________________________________________________________________________

⎧ p = Ae − at A − at ⎪ ⎨ A −at ⇒ x& − ax = − 2 e u e = − ⎪ 2 ⎩ A −at at and x = Be + , the constants A and B will be e 4a calculated from boundary conditions.

⇒

(i). x(0) = 1 and x(1) is not specified Here B + A / 4a = 1 and since x(1) is not specified, we apply the transversality condition (10.32) at t = 1 , which is p (1) = 0 , that is A = 0 . The optimal control u = 0, 0 ≤ t ≤ 1 . This gives J = 0 , which is clearly a minimum as J ≥ 0 . (ii) x(0) = 1 , x(1) = 0 Now, we have B + A / 4a = 1 Be a + Ae − a / 4a = 0 Hence, A = 4a /(1 − e −2a ) and the optimal control is u = −2ae − at /(1 − e −2a ), 0 ≤ t ≤ 1 _________________________________________________________________________________________

___________________________________________________________________________________________________________ 52 Chapter 10 Optimal Control with Unbounded Continuous Controls

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.11 Bang-bang Control 11.1 Introduction This chapter deals with the control with restrictions: is bounded and might well be possible to have discontinuities. To illustrate some of the basic concepts involved when controls are bounded and allowed to have discontinuities we start with a simple physical problem: Derive a controller such that a car move a distance a with minimum time.

dt 2

⎞ ⎟=0 ⇒ ⎟ ⎠

∂F d ⎛ ∂F ⎞ − ⎜ ⎟=0 ∂u dt ⎝ ∂u& ⎠

⇒

∂F d ⎛ ∂F ⎞ − ⎜ ⎟=0 ∂v dt ⎝ ∂v& ⎠

p& 2 = − p1

(11.10)

p 2 = µ ( β − α − 2u )

(11.11)

2vµ = 0

(11.12)

(11.12) ⇒ v = 0 or µ = 0 . We will consider these two cases.

The motion equation of the car

d 2x

∂F d ⎛ ∂F − ⎜⎜ ∂x 2 dt ⎝ ∂x& 2

=u

(11.1)

(i)

µ =0

⎧p = 0 ⇒⎨ 1 ⇒ be impossible. ⎩ p2 = 0 (ii) v = 0

where u = u (t ), − α ≤ u ≤ β

(11.2)

represents the applied acceleration or deceleration (braking) and x the distance traveled. The problem can be stated as minimize

(11.5): v 2 = (u + α )( β − u ) ⇒ u = −α or u = β Hence ⎧β x& 2 = ⎨ ⎩− α

0 ≤ t ≤τ τ
T

∫

T = 1 dt

(11.3)

o

the switch taking place at time τ . Integrating using boundary conditions on x 2

subject to (10.11) and (10.12) and boundary conditions x(0) = 0, x& (0) = 0, x(T ) = a, x& (T ) = 0

(11.4)

The methods we developed in the last chapter would be appropriate for this problem except that they cannot cope with inequality constraints of the form (11.2). We can change this constraint into an equality constraint by introducing another control variable, v, where

v = (u + α )( β − u ) 2

(11.5)

⎧β t x 2 = x&1 = ⎨ ⎩ − α (t − T )

0 ≤ t ≤τ τ
(11.13)

Integrating using boundary conditions on x1 ⎧1 2 ⎪⎪ β t x1 = ⎨ 2 ⎪− 1 α (t − T ) 2 + a ⎪⎩ 2

0 ≤ t ≤τ (11.14)

τ
Since v is real, u must satisfy (11.2). We introduce th usual state variable notation x1 = x so that

Both distance, x1 , and velocity, x 2 , are continuous at t = τ , we must have

x&1 = x2 x1 (0) = 0, x2 (T ) = a x& 2 = u x 2 (0) = 0, x 2 (T ) = 0

(11.14) ⇒ βτ = −α (τ − T )

(11.6) (11.7)

(11.15) ⇒

We now form the augmented functional T* =

∫

T

0

1 1 βτ 2 = a − α (τ − T ) 2 2 2

Eliminating T gives the switching time as {1 + p1 ( x 2 − x&1 ) + p 2 (u − x& 2 ) + µ [v 2

−(u + α )( β − u )]}dt

(11.8)

where p1 , p 2 ,η are Lagrange multipliers associated with the constraints (11.6), (11.7) and (11.5) respectively. The Euler equations for the state variables x1 , x 2 and control variables u and v are ∂F d ⎛ ∂F ⎞ ⎟=0 − ⎜ ∂x1 dt ⎜⎝ ∂x&1 ⎟⎠

⇒

p& 1 = 0

Chapter 11 Bang-bang Control

(11.9)

τ=

2aα β (α + β )

(11.15)

and the final time is T=

2a (α + β )

αβ

(11.16)

The problem now is completely solved and the optimal control is specified by

53

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

⎧β u=⎨ ⎩− α

0 ≤ t ≤τ τ
(11.17)

∂H ∂xi

(i = 1,2, L , n)

(11.23)

The Euler equations for the control variables, u i , do not follow as in section 10.4 as it is possible that there are discontinuities in u i , and so we cannot assume that the partial

this is illustrated in Fig. 11.1 control u

β

p& i = −

derivaties ∂H / ∂u i exist. On the other hand, we can apply the free end point condition (9.23):

0

τ

T

time t

p k (T ) = 0

∂F = 0 to obtain ∂x& i

k = q + 1, L , n

(11.24)

that is, the adjoint variable is zero at every end point where the corresponding state variable is not specified. As before, we refer to (11.24) as tranversality conditions.

−α

Fig. 11.1 Optimal Control Our difficulty now lies in obtaining the analogous equation to ∂H / ∂u i = 0 for continuous controls. For the moment, let us assume that we can differentiate H with respect to u, and consider a small variation δu in the control u such that u + δu still belong to U , the admissible control region. Corresponding to the small change in u , there will be small change in x , say δx , and in p , say δp . The change in the value of J * will be δJ * , where

This figure shows that the control: - has a switch (discontinuity) at time t = τ - only take its maximum and minimum values This type of control is called bang-bang control. 11.2 Pontryagin’s Principle (early 1960s)

Problem: We are seeking extremum values of the functional J =

∫

0

f 0 (x, u, t )dt

(i = 1,2, L , n)

∫

n

(11.19)

The small change operator, δ , obeys the same sort of properties as the differential operator d / dx .

initial conditions x = x 0 and final conditions on x1 , x 2 , L , x q (q ≤ n) and subject to u ∈ U , the admissible control region. For example, in the previous problem, the admissible control region is defined by

{H −

o

Assuming we can interchange the small change operator, δ , and integral sign, we obtain T

∫

δJ * =

U = {u : −α ≤ u ≤ β }

n

[δ ( H −

T

∫

0

⎡ ⎢ f0 + ⎢ ⎣

n

∑ i =1

⎤ pi ( f i − x& i )⎥ dt ⎥ ⎦

T

∫

=

δH =

∫

0

⎞

n

∑ p x& ⎟⎟ dt i i

i =1

(11.22)

⎠

and evaluating the Euler equations for xi , we obtain as in section 10.4 the adjoint equations

Chapter 11 Bang-bang Control

∑ p δ x& ]dt i

i

i =1

∂H δu j + ∂u j

n

∑ i =1

∂H δxi + ∂xi

n

∂H

∑ ∂p δp i =1

i

i

so that

δJ * =

⎡ ⎢ ⎢ ⎣

m

∂H

∫ ∑ ∂u T

0

⎜H − ⎜ ⎝

∑

n

x& i δ pi −

i =1

∑ j =1

For simplicity, we consider that the Hamintonian is a function of the state vector x, control vector u, and adjoint vector p, that is, H = H (x, u, p) . We can express J * as T⎛

[δH −

m

(11.21)

pi f i

i =1

J* =

i i

i =1 n

0

n

∑

T

Using chain rule for partial differentiation

and define the Hamintonian H = f0 +

i i

i =1 n

0

(11.20)

∑ p x& )]dt

∫ [δH − ∑ δ ( p x& )]dt

=

As in section 10.4, we form the augmented functional

i i

i =1

0

J* =

∑ p x& }dt

δJ * = δ

subject to state equations x& i = f i(x, u, t )

T

(11.18)

T

j =1

n

+

δu j

j

⎛ ∂H

∑ ⎜⎜⎝ ∂p δ p i =1

i

i

⎞⎤ − p& i δ xi − xi δ pi − pi δ x& i ⎟⎟⎥ dt ⎠⎥⎦

since p& i = −∂H / ∂xi . Also, from (11.21)

∂H = f i = x& i ∂pi

54

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

using (11.19): x& i = f i(x, u, t ), (i = 1,2, L , n) . Thus ⎡ ⎢ 0 ⎢ ⎣ T⎡ ⎢ = 0 ⎢ ⎣

δJ * =

m

∫ ∑ T

j =1 m

∫ ∑ j =1

∂H δu j − ∂u j ∂H δu j − ∂u j

We first illustrate its use by examining a simple problem. We required to minimize

⎤

n

∑ ( p& iδ xi + piδ x&i )⎥⎥ dt ∑ i =1

0

⎦

i =1 n

T

∫

J = 1dt

⎤ d ( piδ xi )⎥⎥ dt dt ⎦

subject to x&1 = x 2 , x& 2 = u where −α ≤ u ≤ β and x1 (T ) = a , x2 (T ) = 0 , x1 (0) = x 2 (0) = 0 .

We can now integrate the second part of the integrand to yield T

n

δJ * = −

T ⎛⎜ m

∑ ( p δ x ) + ∫ ⎜⎜ ∑ i

i

i =1

0

0

⎝

j =1

⎞ ∂H ⎟ δu j ⎟ dt ∂u j ⎟ ⎠

At t = 0 :

xi (i = 1,L, n) are specified ⇒ δxi (0) = 0

At t = T :

xi (i = 1,L , q ) are fixed

H = 1 + p1 x 2 + p 2 u

(11.25)

⇒ δxi (T ) = 0

For i = q + 1, L, n, from the transversality conditions, (11,24)

for i = 1,2, L , q, q + 1, L , n . We now have m

∫ ∑ j =1

We must minimize H with respect to u, and where u ∈U = [−α , β ] , the admissible control region.. Since H is linear in u, it clearly attains its minimum on the boundary of the control region, that is, either at u = −α or u = β . This illustrated in Fig.11.3. In fact we can write the optimal control as ⎧−α u=⎨ ⎩ β

pi (T ) δxi (T ) = 0

T ⎛⎜ δJ * = ⎜ 0 ⎜ ⎝

Introducing adjoint variables p1 and p 2 , the Hamiltonian is given by

H

control vector u . Since all these variations are independent, and we require δJ * = 0 for a turning point when the controls are continuous, we conclude that ∂H = 0 ( j = 1,2, L , m) ∂u j

(11.26)

But this is only valid when the controls are continuous and not constrained. In our present case when u ∈ U , the admissible control region and discontinuities in u are allowed. The arguments presented above follow through in the same way, except that (∂H / ∂u j )du j must be replaced by

control u

We thus obtain

∫∑ 0

[ H (x; u1 ,L, u j + δu j ,L, u m ; p) − H (x, u, p)] dt

j =1

In order for u to be a minimizing control, we must have δJ * ≥ 0 for all admissible controls u + δu . This implies that H (x; u1 ,L, u j + δu j ,L, u m ; p) ≥ H (x, u, p)

(11.27)

for all admissible δu j and for j = 1,L, m . So we have established that on the optimal control H is minimized with respect to the control variables, u1 , u 2 ,L, u m . This is known as Pontryagin’s minimum principle. Chapter 11 Bang-bang Control

−α

β Fig. 11.3 The case p 2 > 0

But p 2 will vary in time, and satisfies the adjoint equations, p& 1 = −

∂H =0 ∂x1

p& 2 = −

∂H = − p1 ∂x 2

Thus p1 = A , a constant, and p 2 = − At + B , where B is constant. Since p 2 is a linear function of t, there will at most be one switch in the control, since p 2 has at most one zero, and from the physical situation there must be at least one switch. So we conclude that

H (x; u1 , u 2 ,L, u j + δu j ,L, u m ; p) − H (x, u, p)

T m

if p 2 < 0

⎞ ∂H ⎟ δu j ⎟ dt ∂u j ⎟ ⎠

where δu j is the small variation in the j th component of the

δJ * =

if p 2 > 0

(i) the control u = −α or u = β , that is, bang-bang control; (ii) there is one and only one switch in the control. Again, it is clear from the basic problem that initially u = β , followed by u = −α at the appropriate time.

11.3 Switching Curves In the last section we met the idea of a switch in a control. The time (and position) of switching from one extremum value of the control to another does of course depend on the initial starting point in the phase plane. Byconsidering a specific example we shall show how these switching positions define a switching curve.

55

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

Suppose a system is described by the state variables x1 , x 2 where x&1 = − x1 + u x& 2 = u

(11.28) (11.29)

Here u is the control variable which is subject to the constraints −1 ≤ u ≤ 1 . Given that at t = 0, x1 = a, x 2 = b , we wish to find the optimal control which takes the system to x1 = 0, x 2 = 0 in minimum time; that is, we wish to minimize T

∫

J = 1dt

(11.30)

x 2 − B = − k log | ( x1 − k ) / A |

Now if u = 1 , that is, k = 1 , then the trajectories are of the form x 2 − B = − log | ( x1 − 1) / A |

that is x 2 = − log | x1 − 1 | +C

(11.32)

where C is constant. The curves for different values of C are illustrated in Fig.11.4

0

u =1

while moving from (a, b) to (0,0) in the x1 − x 2 phase plane and subject to (11.28), (11.29) and −1 ≤ u ≤ 1

x2

(11.31)

Following the procedure outline in section 11.2, we introduce the adjoint variables p1 and p 2 and the Hamiltonian

x1

0

H = 1 + p1 (− x1 + u ) + p 2 u

1

= u ( p1 + p 2 ) + 1 − p1 x1 Since H is linear in u, and | u |≤ 1 , H is minimized with respect to u by taking ⎧+1 if u=⎨ ⎩− 1 if

p1 + p 2 < 0

A

p1 + p 2 > 0

So the control is bang-bang and the number of switches will depend on the sign changes in p1 + p 2 . As the adjoint equations, (11.23), are p& 1 = −

∂H = p1 ∂x1

∂H p& 2 = − =0 ∂x 2

Fig. 11.4 Trajectories for u = 1 Follow the same procedure for u = −1 , giving x2 = log | x1 − 1 | +C

⇒

p1 = Ae p2 = B

(11.33)

t

, A and B are constant

and the curves are illustrated in Fig. 11.5.

B

x2

and p1 + p 2 = Ae t + B , and this function has at most one sign change. So we know that from any initial point (a, b) , the optimal control will be bang-bang, that is, u = ±1 , with at most one switch in the control. Now suppose u = k , when k = ±1 , then the state equations for the system are

x1 −1

x&1 = − x1 + k x& 2 = k

0

u = −1

We can integrate each equation to give x1 = Ae −t + k x2 = k t + B

, A and B are constants

The x1 − x 2 plane trajectories are found by eliminating t , giving

Chapter 11 Bang-bang Control

Fig. 11.4 Trajectories for u = −1

56

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

The basic problem is to reach the origin from an arbitrary initial point. All the possible trajectories are illustrated in Figs. 11.4 and 11.5, and we can see that these trajectories are only two possible paths which reach the origin, namely AO in Fig. 11.4 and BO in Fig. 11.5. x2

Following the usual procedure, we form the Hamiltonian H = 1 + p1 x2 + p 2 u

(11.39)

We minimize H with respect to u, where −1 ≤ u ≤ 1 , which gives

B

⎧+1 if u=⎨ ⎩− 1 if

u = −1

p2 < 0 p2 > 0

The adjoint variables satisfy x1

0

−1

1

∂H ⎧& ⎪ p1 = − ∂x = 0 ⎧p = A ⎪ 1 ⇒ ⎨ 1 ⎨ H ∂ ⎩ p 2 = − At + B ⎪ p& = − = − p1 ⎪⎩ 2 ∂x2 Since x2 (T ) is not specified, the transversality condition

u =1

A

becomes p 2 (T ) = 0 . Hence 0 = − At + B and p 2 = A(T − t ) . For 0 < t < T , there is no change in the sign of p 2 , and hence no switch in u. Thus either u = +1 or u = −1 , but with no switches. We have

Fig. 11.4 Trajectories for u = −1 Combining the two diagrams we develop the Fig. 11.6. The curve AOB is called switching curve. For initial points below AOB, we take u = +1 until the switching curve is reached, followed by u = −1 until the origin is reached. Similarly for the points above AOB, u = −1 until the switching curve is reached, followed by u = +1 until the origin is reached. So we have solved the problem of finding the optimal trajectory from an arbitrary starting point.

x22 = 2( x1 − B) x22

= −2( x1 − B)

when

u =1

(11.40)

when

u = −1

(11.41)

These trajectories are illustrated in Fig. 11.7, the direction of the arrows being determined from x& 2 = u u = +1

x2

x2

u = −1

Thus the switching curve has equation ⎧ log(1 + x1 ) x2 = ⎨ ⎩− log(1 + x1 )

for

x1 > 0

for

x1 < 0

(11.34)

x1

0

x1

0

11.4 Transversarlity conditions To illustrate how the transversality conditions ( pi (T ) = 0 if xi is not specified) are used, we consider the problem of

finding the optimum control u when | u | ≤ 1 for the system described by x&1 = x 2 x& 2 = u

Fig. 11.7 Possible Optimal Trajectories x2 x1

(11.35) 0

(11.36)

( a, b )

which takes the system from an arbitrary initial point x1 (0) = a, x 2 (0) = b to any point on the x2 axis, that is, x1 (T ) = 0 but x2 (T ) is not given, and minimize

∫

J = 1.dt

{

T+

T

(11.37)

T−{

( a, b )

0

subject to (11.35), (11.36), the above boundary conditions on x1 and x 2 , and such that −1 ≤ u ≤ 1

Chapter 11 Bang-bang Control

(11.38)

A

Fig. 11.8 Initial point below OA We first consider initial points (a, b) for which a > 0 . For points above the curve OA, there is only one trajectory,

57

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

u = −1 which reaches the x 2 -axis, and this must be the optimal curve. For points below the curve OA, there are two possible curves, as shown in Fig. 11.8. From (11.36): x& 2 = u , that is x& 2 = ±1 , and the integrating between 0 and T gives x 2 (T ) − x 2 (0) = ±T

that is T = | x 2 (T ) − x 2 (0) |

(11.42)

Hence the modulus of the difference in final and initial values of x 2 given the time taken. This is shown in the diagram as T+ for u = +1 and T− for u = −1 . The complete set of optimal trajectories is illustrated in Fig. 11.9. x2 u = +1

u = −1

0

x1

Fig. 11.9 Optimal Trajectories to reach x 2 -axis in minimum time

11.5 Extension to the Boltza problem

Chapter 11 Bang-bang Control

58

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

C.12 Applications of Optimal Control Problem Minimizing the fuel (cost function)

12.6 Rocket Trajectories The governing equation of rocket motion is

J=

dv m = m F + Fext dt

(8.22)

T

∫ β dt

(12.41)

o

subject to the differential constraints where m v F Fext

dv1 cβ = cos φ dt m dv 2 cβ = sin φ − g dt m

: rocket’s mass : rocket’s velocity : thrust produced by the rocket motor : external force

It can be seen that

(12.43)

and also

c β m

F =

(12.42)

(8.23)

where c : relative exhaust speed β = −dm / dt : burning rate

dx = v1 dt dy = v2 dt

(12.44) (12.45)

The boundary conditions

Let φ be the thrust attitude angle, that is the angle between the rocket axis and the horizontal, as in the Fig. 8.4, then the equations of motion are

t = 0 : x = 0, y = 0, v1 = 0, v 2 = 0 t = T : x not specified, y = y , v1 = v , v 2 = 0 Thus we have four state variables, namely x, y , v1 , v 2 and two

dv1 cβ = cos φ dt m dv 2 cβ = sin φ − g dt m

(8.24)

Here v = (v1 , v 2 ) and the external force is the gravity force only Fext = (0,− mg ) .

controls β (the rate at which the exhaust gases are emitted) and φ (the thrust attitude angle). In practice we must have bounds on β , that is, 0≤β ≤β

(12.46)

so that β = 0 corresponds to the rocket motor being shut

y

down and β = β corresponds to the motor at full power.

path of rocket

The Hamilton for the system is r v

cβ ⎛ cβ ⎞ cos φ + p 4 ⎜ sin φ − g ⎟ m ⎝ m ⎠ (12.47) where p1 , p 2 , p 3 , p 4 are the adjoint variables associated H = β + p1v1 + p 2 v 2 + p 3

φ x

Fig. 8.4 Rocket Flight Path

(12.47) ⇒

The minimum fuel problem is to choose the controls, β and φ , so as to take the rocket from initial position to a prescribed height, say, y , in such a way as to minimize the fuel used. The fuel consumed is T

∫ β dt o

where T is the time at which y is reached.

Chapter 12 Applications of Optimal Control

with x, y , v1 , v 2 respectively.

(8.25)

c c ⎛ ⎞ H = p1v1 + p 2 v 2 − p 4 g + β ⎜1 + p 3 cos φ + p 4 sin φ ⎟ m m ⎝ ⎠ c c ⎛ ⎞ If we assuming that ⎜1 + p 3 cos φ + p 4 sin φ ⎟ ≠ 0 , we m m ⎝ ⎠ see that H is linear in the control β , so that β is bang-bang. That is β = 0 or β = β . We must clearly start with β = β so that, with one switch in β , we have

59

Introduction to Control Theory Including Optimal Control

Nguyen Tan Tien - 2002.5

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

⎧⎪ β ⎪⎩0

β =⎨

0 ≤ t < t1

for

(12.48)

for t1 < t ≤ T

Now H must also be maximized with respect to the second control, φ , that is, ∂H / ∂φ = 0 giving − p3

cβ cβ sin φ + p 4 cos φ = 0 m m

p& 1 = −

∂H = − p1 = − A ∂v1

∂H p& 1 = − = − p2 ∂v 2

(12.50)

subject to d 2θ 2

+α

dθ + ω 2θ = u dt

(12.51)

Boundary conditions

t = 0 : θ = θ 0 , θ& = θ 0′ t = T : θ = 0, θ& = 0

The adjoint variables satisfy the equations ∂H =0 ∂x ∂H p& 1 = − =0 ∂y

T

0

dt

which yields tan φ = p 4 / p 3 .

p& 1 = −

∫

T = 1dt

⇒

p1 = A

Constraint on control: | u | ≤ 1 .

⇒

p2 = B

Introduce the state variables: x1 = θ , x 2 = θ& . Then (12.51) becomes

⇒

p 3 = C − At

⇒

p 4 = D − Bt

x&1 = x 2 x& 2 = a x 2 − ω 2 x1 + u

(12.52)

As usual, we form the Hamiltonian

where, A, B, C, D are constant. Thus H = 1 + p1 y 2 + p 2 (−a x 2 − ω 2 x1 + u )

D − Bt tan φ = C − At

where p1 , p 2 are adjoint variables satisfying

Since x is not specified at t = T , the transversality condition is p1 (T ) = 0 , that is, A = 0 , and so tan φ = a − bt

(12.49)

where a = D / C , b = B / C .

r v

y=y

β =0 β =β

∂H = ω 2 p2 ∂x ∂H p& 1 = − = − p1 + ap 2 ∂y p& 1 = −

(12.54) (12.55)

Since H is linear in u, and | u | ≤ 1 , we again immediately see

The problem, in principle, I now solved. To complete it requires just integration and algebraic manipulation. With φ given by (12.49) and β given by (12.48), we can integrate (12.42) to (12.45). This will bring four further constants of integration; together with a, b and the switchover time t1 we have seven unknown constants. These are determined from the seven end-point conditions at t = 0 and t = T . A typical trajectory is shown in Fig. 12.5. y

(12.53)

⎧+1 if u=⎨ ⎩− 1 if

p2 < 0 p2 > 0

To ascertain the number of switches, we must solve for p 2 . From (12.54) and (12.55), we see that &p& 2 = − p& 1 + a p& 2 = −ω 2 p 2 + a p& 2

that is,

null thrust arc

t = t1

&p&2 − a p& 2 + ω 2 p 2 = 0

(12.56)

The solution of (12.56) gives us the switching time.

maximum thrust arc 0

that the control u is bang-bang, that is, u = ±1 , in fact

x

Fig. 12.5 Typical rocket trajectory

12.7 Servo Problem The problem here is to minimize

Chapter 12 Applications of Optimal Control

59