Chapter 2 State Space Models General format : (valid for any - - PDF document

✬ ✫ ✩ ✪

Chapter 2 State Space Models

General format :

(valid for any nonlinear causal system)

CT : ˙ x = f(x, u), DT : xk+1 = f(xk, uk), y = h(x, u). yk = h(xk, uk). where x = the state of the system, n × 1-vector u = the input of the system, m × 1-vector y = the output of the system, p × 1-vector f = state equation vector function h = output equation vector function n = number of states ⇒ n-th order system m = number of inputs p = number of outputs Single-Input Single-Output system (SISO) : m = p = 1, Multi-Input Multi-Output system (MIMO) : m, p > 1, Multi-Input Single-Output system (MISO) : m > 1, p = 1, Single-Input Multi-Output system (SIMO) : m = 1, p > 1.

ESAT–SCD–SISTA

CACSD

• pag. 34

✬ ✫ ✩ ✪

Linear Time Invariant (LTI) systems

(only valid for linear causal systems)

Continuous-time system: ˙ x = Ax + Bu y = Cx + Du

u B ˙ x x C D A

• y

Discrete-time system: xk+1 = Axk + Buk yk = Cxk + Duk

uk B xk+1 xk C D A z−1 yk

A: n × n system matrix B: n × m input matrix C: p × n output matrix D: p × m direct transmission matrix

ESAT–SCD–SISTA

CACSD

• pag. 35

✬ ✫ ✩ ✪

An example of linear state space modeling:

Example: Tape drive control - state space modeling Process description:

p1 p3 p2 r r θ1 θ2 K D K D Kt Kt β β J J i1 i2 + + − − e1 e2

The drive motor on each end of the tape is independently controllable by voltage resources e1 and e2. The tape is modeled as a linear spring with a small amount of viscous damping near to the static equilibrium with tape tension

• 6N. The variables are defined as deviations from this equi-

librium point.

ESAT–SCD–SISTA

CACSD

• pag. 36

✬ ✫ ✩ ✪

The equations of motion of the system are: capstan 1: J dω1 dt = Tr − βω1 + Kti1, speed of x1: ˙ p1 = rω1, motor 1: Ldi1 dt = −Ri1 − Keω1 + e1, capstan 2: J dω2 dt = −Tr − βω + Kti2, speed of x1: ˙ p2 = rω2, motor 2: Ldi2 dt = −Ri2 − Keω2 + e2, Tension of tape: T = K 2 (p2 − p1) + D 2 ( ˙ p2 − ˙ p1), Position of the head: p3 = p1 + p2 2 .

ESAT–SCD–SISTA

CACSD

• pag. 37

✬ ✫ ✩ ✪

Description of variables and constants: D = damping in the tape-stretch motion = 20 N/m·sec, e1,2 = applied voltage, V, i1,2 = current into the capstan motor, J = inertia of the wheel and the motor = 4×10−5kg·m2, β = capstan rotational friction, 400kg·m2/sec, K = spring constant of the tape, 4×104N/m, Ke = electrical constant of the motors = 0.03V·sec, Kt = torque constant of the motors = 0.03N·m/A, L = armature inductance = 10−3 H, R = armature resistance = 1Ω, r = radius of the take-up wheels, 0.02m, T = tape tension at the read/write head, N, p1,2,3 = tape position at capstan 1,2 and the head, ˙ p1,2,3 = tape velocity at capstan 1,2 and the head, θ1,2 = angular displacement of capstan 1,2, ω1,2 = speed of drive wheels = ˙ θ1,2.

ESAT–SCD–SISTA

CACSD

• pag. 38

✬ ✫ ✩ ✪

With a time scaling factor of 103 and a position scaling fac- tor 10−5 for numerical reasons, the state equations become: ˙ x = Ax + Bu, y = Cx + Du, where x =            p1 ω1 p2 ω2 i1 i2            , A =            2 −0.1 −0.35 0.1 0.1 0.75 2 0.1 0.1 −0.1 −0.35 0 0.75 0 −0.03 −1 0 −0.03 −1            B =            0 0 0 0 0 0 0 0 1 0 0 1            , C =

• 0.5

0 0.5 0 0 0 −0.2 −0.2 0.2 0.2 0 0

• ,

D =

• 0 0

0 0

• , y =
• p3

T

• , u =
• e1

e2

• .

ESAT–SCD–SISTA

CACSD

• pag. 39

✬ ✫ ✩ ✪

State-space model, transfer matrix and impulse response

Continuous-time system:

state-space equations

• ˙

x = Ax + Bu y = Cx + Du Laplace

x(0)=0

− → G(s) = Y (s) U(s) = D + C(sI − A)−1B

• transfer matrix

in practice : D = 0

⇓ G(t) = CeAtB

• impulse response matrix

Laplace ← → G(s) =

∞

• i=1

CAi−1Bs−i

• transfer matrix

Discrete-time system:

state-space equations

• xk+1 = Axk + Buk

yk = Cxk + Duk z − transf

x0=0

← → G(z) = Y (z) U(z) = D + C(zI − A)−1B

• transfer matrix

⇓ G(k) =

• D

: k = 0 CAk−1B : k ≥ 1

• impulse response matrix

z − transf ← → G(z) =

∞

• i=1

CAi−1Bz−i + D

• transfer matrix

In case of a SISO system, G(t) or G(k) is the impulse response. For MIMO G(t), G(k) are matrices containing the m × p impulse responses (one for every input-output pair).

ESAT–SCD–SISTA

CACSD

• pag. 40

✬ ✫ ✩ ✪

Linearization of a nonlinear system about an equilibrium point

Consider a general nonlinear system in continuous time : dx dt = f(x, u) y = h(x, u) For small deviations about an equilibrium point (xe, ue, ye) for which f(xe, ue) = 0 ye = h(xe, ue) we define x = xe + ∆x, u = ue + ∆u, y = ye + ∆y, and obtain dx dt = d∆x dt = f(x, u) = f(xe + ∆x, ue + ∆u) and ye + ∆y = h(x, u) = h(xe + ∆x, ue + ∆u).

ESAT–SCD–SISTA

CACSD

• pag. 41

✬ ✫ ✩ ✪

By first order approximation we obtain a linear state space model from ∆u to ∆y : d∆x dt = f(xe + ∆x, ue + ∆u) ⇓ f(xe, ue) = 0 d∆x dt

(1)

= ∂f ∂x

• xe,ue

n × n ⇓ A ∆x + ∂f ∂u

• xe,ue

n × m ⇓ B ∆u and ye + ∆y = h(xe + ∆x, ue + ∆u) ⇓ ye = h(xe, ue) ∆y

(1)

= ∂h ∂x

• xe,ue

p × n ⇓ C ∆x + ∂h ∂u

• xe,ue

p × m ⇓ D ∆u Conclusion : use a linear approximation about equilibrium.

ESAT–SCD–SISTA

CACSD

• pag. 42

✬ ✫ ✩ ✪

Example : Consider a decalcification plant which is used to reduce the concentration of calcium hydroxide in water by forming a calcium carbonate precipitate. The following equations hold (simplified model) :

• chemical reaction : Ca(OH)2 + CO2 → CaCO3 + H2O
• reaction speed : r = c[Ca(OH)2][CO2]
• rate of change of concentration :

d[Ca(OH)2] dt = k V − r V d[CO2] dt = u V − r V k and u are the inflow rates in moles/second of calcium hydroxide and carbon dioxide respectively. V is the tank volume in liters. Let the inflow rate of carbon dioxide be the input and the concentration of calcium hydroxide be the output.

ESAT–SCD–SISTA

CACSD

• pag. 43

✬ ✫ ✩ ✪

A nonlinear state space model for this reactor is : d[Ca(OH)2] dt = k V − c V [Ca(OH)2][CO2] d[CO2] dt = u V − c V [Ca(OH)2][CO2] y = [Ca(OH)2] The state variables are x1 = [Ca(OH)2] and x2 = [CO2]. In equilibrium we have : k V − c V [Ca(OH)2]eq[CO2]eq = k V − c V X1X2 = 0 ueq V − c V [Ca(OH)2]eq[CO2]eq = U V − c V X1X2 = 0 Y = [Ca(OH)2]eq = X1 For small deviations about the equilibrium : d∆x1 dt = − c V X2∆x1 − c V X1∆x2 d∆x2 dt = − c V X2∆x1 − c V X1∆x2 + 1 V ∆u ∆y = ∆x1 so, A =

• −cX2

V

−cX1

V

−cX2

V

−cX1

V

• , B =
• 1

V

• , C =
• 1 0
• , D = 0

ESAT–SCD–SISTA

CACSD

• pag. 44

✬ ✫ ✩ ✪

If k = 0.1 mole sec , c = 0.5

l2

sec· mole, U = 0.1 mole sec , X1 = 0.25 mole l , X2 = 0.8 mole l , V = 5 l then

• A B

C D

• =

   −0.08 −0.025 −0.08 −0.025 0.2 1    The corresponding transfer function is ∆y(s) ∆u(s) = −0.005 s2 + 0.105s and its bode plot

Frequency (rad/sec) Phase (deg); Magnitude (dB) Bode Diagrams −50 50 10

−3

10

−2

10

−1

10 −400 −350 −300 −250

ESAT–SCD–SISTA

CACSD

• pag. 45

✬ ✫ ✩ ✪

One can also obtain a linear state space model for this chemical plant from linear system identification. For small deviations about the equilibrium point the dy- namics can be described fairly well by a linear (A, B, C, D)-

• model. Hence, a small white noise disturbance ∆u was

generated and was added to the equilibrium value U. We applied this signal to the input of a nonlinear model of the chemical reactor (Simulink model for instance), i.e. u(t) = U + ∆u. The following input-output set was obtained :

100 200 300 400 500 600 700 800 900 1000 0.096 0.098 0.1 0.102 0.104 u(t) time (s) 100 200 300 400 500 600 700 800 900 1000 0.249 0.2495 0.25 0.2505 time (s) y(t)

ESAT–SCD–SISTA

CACSD

• pag. 46

✬ ✫ ✩ ✪

By applying a linear system identification algorithm (N4SID), the following 2nd-order model was obtained :

• A B

C D

• =

   0.0015 −0.0589 −0.0867 0.0027 −0.1066 0.1713 0.2546 0.129    The corresponding transfer function is y(s) u(s) = −1.299.10−7s − 0.005 s2 + 0.1051s + 1.346.10−8 It has 2 poles at −1.281.10−7 and −1.051. As −1.281.10−7 lies close to 0, and taking in account the properties of the manually derived linear model of page 45, we conclude that the plant has one integrator pole. Hence, it might be better to fix one pole at s = 0. In this way it is guaranteed that the linear model obtained by system identification is stable. The transfer function which was obtained using this mod- ified identification procedure is y(s) u(s) = −1.68.10−8s − 0.005 s2 + 0.1051s

ESAT–SCD–SISTA

CACSD

• pag. 47

✬ ✫ ✩ ✪

The different models are validated by comparing their re- sponse to the following input : ∆u(t) = 0.003 if 100 < t < 200 = 0 elsewhere Four responses are shown :

• the nonlinear system (—)
• the linearised model obtained by hand (page 45) (- -)
• the 2nd-order linear model obtained from N4SID (· · ·)
• the linear model obtained from N4SID having a fixed

integrator pole by construction (-.)

100 200 300 400 500 600 700 800 900 1000 0.235 0.24 0.245 0.25 y(t) mole/l 900 905 910 915 920 925 930 935 940 945 950 0.2356 0.2358 0.236 0.2362 0.2364 y(t) : detail 1 mole/l 900 905 910 915 920 925 930 935 940 945 950 0.23572 0.23573 0.23574 mole/l y(t) : detail 2 time(s)

ESAT–SCD–SISTA

CACSD

• pag. 48

✬ ✫ ✩ ✪

Digitization

in this course we want to control physical plants using a digital computer :

Plant A/D controller digital sensors clock actuators + D/A computer

• analog-to-digital converter :

clock A/D x(t) x[k] x(t) filter lowpass sampler quantizer

Q

x[k]

– the analog anti-aliasing filter filters out high frequent components (> Nyquist frequency) – the filtered analog signal is sampled and quantized in order to obtain a digital signal. For quantization 10 to 12 bits are common.

ESAT–SCD–SISTA

CACSD

• pag. 49

✬ ✫ ✩ ✪

• digital-to-analog converter :

– zero-order hold

10 20 30 40 50 60 70 80 90 100 −4 −3 −2 −1 1 2 3 time ZOH approximation

∗ introduces a delay ∗ frequency spectrum deformation – first-order hold – n-th order polynomial (more general case) : fit a n- th order polynomial through the n+1 most recent samples and extrapolate to the next time instance

ESAT–SCD–SISTA

CACSD

• pag. 50

✬ ✫ ✩ ✪

Digitization of state space models

• 1. discretization by applying numerical integration rules :

a continuous-time integrator ˙ x(t) = e(t) ⇔ sX(s) = E(s) can be approximated using

• the forward rectangular rule or Euler’s method

xk+1 − xk Ts = ek ⇔ z − 1 Ts X(z) = E(z) with xk = x(kTs)

t e(t)

k-1 k k+1

• the backward rectangular rule

xk+1 − xk Ts = ek+1 ⇔ z − 1 zTs X(z) = E(z)

t e(t)

k-1 k k+1

ESAT–SCD–SISTA

CACSD

• pag. 51

✬ ✫ ✩ ✪

• the trapezoid rule or bilinear transformation

xk+1 − xk Ts = ek+1 + ek 2 ⇔ 2 Ts z − 1 z + 1X(z) = E(z)

t e(t)

k-1 k k+1

Change a continuous model G(s) into a discrete model Gd(z) by replacing all integrators with their discrete equiv- alents : = projection of the left half-plane Euler backward rect. bilinear transf. s → z−1

Ts

s → z−1

zTs

s → 2

Ts z−1 z+1

Re Im Re Im Re Im

Except for the forward rectangular rule stable continuous poles (grey zone) are guaranteed to be placed in stable discrete areas, i.e. within the unit circle. The bilinear transformation maps stable poles → stable poles and unstable poles → unstable poles.

ESAT–SCD–SISTA

CACSD

• pag. 52

✬ ✫ ✩ ✪

The following continuous-time model is given : ˙ x = Ax + Bu y = Cx + Du ⇔ sX = AX + BU Y = CX + DU following Euler’s method s is replaced by z−1

Ts , so

z − 1 Ts X = AX + BU

• r

zX = (I + ATs)X + BTsU the output equation Y = CX + DU remains. A similar calculation can be done for the backward rectangular rule and the bilinear transformation resulting in the following table : Euler backward rect. bilinear transf.

Ad I + ATs (I − ATs)−1 (I − ATs

2 )−1(I + ATs 2 )

Bd BTs (I − ATs)−1BTs (I − ATs

2 )−1BTs

Cd C C(I − ATs)−1 C(I − ATs

2 )−1

Dd D D + C(I − ATs)−1BTs D + C(I − ATs

2 )−1BTs 2

ESAT–SCD–SISTA

CACSD

• pag. 53

✬ ✫ ✩ ✪

• 2. discretization by assuming zero-order hold ...

˙ x = Ax + Bu ⇒ x(t) = eAtx(0) + t eA(t−τ)Bu(τ)dτ

• convolution integral

Let the sampling time be Ts, then x(t + Ts) = eA(t+Ts)x(0) + t+Ts eA(t+Ts−τ)Bu(τ)dτ = eATsx(t) + eATs t+Ts

t

eA(t−τ)Bu(τ)dτ

• approximated

ZOH

⇒ xk+1

t=kTs

= eATsxk + eATsA−1(I − e−ATs)Buk = eATs

• Ad

xk + A−1(eATs − I)B

• Bd

uk One can prove that Gd(z) can be expressed as Gd(z)

ZOH

= (1 − z−1)Z

• L−1

G(s) s

• For this reason this discretization method is sometimes

called step invariance mapping.

ESAT–SCD–SISTA

CACSD

• pag. 54

✬ ✫ ✩ ✪

• 3. discretization by zero-pole mapping :

following the previous method the poles of Gd(z) are related to the poles of G(s) according to z = esTs. If we assume by simplicity that also the zeros undergo this transformation, the following heuristic may be applied: (a) map poles and zeros according to z = esTs (b) if the numerator is of lower order than the denomi- nator, add discrete zeros at -1 until the order of the numerator is one less than the order of the denomi-

• nator. A lower numerator order corresponds to zeros

at ∞ in continuous time. By discretization they are put at -1. (c) adjust the DC gain such that lim

s→0 G(s) = lim z→1 Gd(z) ESAT–SCD–SISTA

CACSD

• pag. 55

✬ ✫ ✩ ✪

Example : given the following SISO system

• A B

C D

• =

   −0.2 −0.5 1 1 1 0.4    ⇒ G(s) = C(sI − A)−1B + D = s + 0.4 s2 + 0.2s + 0.5. Now compare the following discretization methods (Ts = 1 sec.) :

• 1. bilinear transformation :

(a) Ad, Bd, Cd and Dd are calculated using the conver- sion table (b) Gbilinear(z) = Cd(zI − Ad)−1Bd + Dd. = 0.4898 + 0.1633z−1 − 0.3265z−2 1 − 1.4286z−1 + 0.8367z−2

• 2. ZOH :

(a) Ad = eATs, Bd = A−1(eATs − I)B, Cd = C and Dd = D (b) GZOH(z) = Cd(zI − Ad)−1Bd + Dd. = 1.0125z−1 − 0.6648z−2 1 − 1.3841z−1 + 0.8187z−2

ESAT–SCD–SISTA

CACSD

• pag. 56

✬ ✫ ✩ ✪

• 3. pole-zero mapping :

(a) the poles of G(s) are −0.1 + j0.7 and −0.1 − j0.7 (b) there is one zero at −0.4 (c) Gpz = K z − e−0.4 (z − e−0.1+j0.7)(z − e−0.1−j0.7) = K z−1 − 0.6703z−2 1 − 1.3841z−1 + 0.8187z−2 hence K = lims→0 G(s) limz→1 Gpz(z) = 1.0546 Compare the bode plots :

10

−2

10

−1

10

−1

10 10

1

frequency (Hz) frequency amplitude spectrum __ : continuous system −− : bilinear approximation −. : zero−order hold assumption .. : pole−zero mapping 10

−2

10

−1

−150 −100 −50 50 frequency (Hz) phase spectrum (degr.) __ : continuous system −− : bilinear approximation −. : zero−order hold assumption .. : pole−zero mapping

ESAT–SCD–SISTA

CACSD

• pag. 57

✬ ✫ ✩ ✪

Sampling rate selection :

• the lower the sampling rate, the lower the implementa-

tion cost (cheap microcontroller & A/D converters), the rougher the response and the larger the delay between command changes and system response.

• an absolute lower bound is set by the specification to

track a command input with a certain frequency : fs ≥ 2BWcl fs is the sampling rate and BWcl is the closed-loop sys- tem bandwidth.

• if the controller is designed for disturbance rejection,

fs > 20BWcl

• when the sampling rate is too high, finite-word effects

show up in small word-size microcontrollers (< 10 bits).

• for systems where the controller adds damping to a

lightly damped mode with resonant frequency fr, fs > 2fr In practice, the sampling rate is a factor 20 to 40 higher than BWcl.

ESAT–SCD–SISTA

CACSD

• pag. 58

✬ ✫ ✩ ✪

Advantages of state space models

• More general models: LTI and Nonlinear Time Varying

(NTV).

• Geometric concepts: more mathematical tools (linear

algebra).

• Internal and external descriptions: “divide and con-

quer” strategy.

• Unified framework: the same for SISO and MIMO.

ESAT–SCD–SISTA

CACSD

• pag. 59

✬ ✫ ✩ ✪

Geometric properties of linear state space models Canonical forms

Control canonical form (SISO): ˙ x = Acx + Bcu, y = Ccx. Ac =         −a1 −a2 · · · · · · −an 1 · · · · · · ... ... . . . . . . ... ... ... . . . · · · 1         , Bc =       1 . . .       , Cc =

• b1 b2 · · · bn
• Corresponding transfer function:

G(s) = b(s) a(s) b(s) = b1sn−1 + b2sn−2 + · · · + bn a(s) = sn + a1sn−1 + a2sn−2 + · · · + an

ESAT–SCD–SISTA

CACSD

• pag. 60

✬ ✫ ✩ ✪

Block diagram for the control canonical form

x1 x2 xn−1 b1 b2 bn−1 bn −a1 −a2 −an−1 −an u y

• xn

ESAT–SCD–SISTA

CACSD

• pag. 61

✬ ✫ ✩ ✪

Modal canonical form: ˙ x = Amx + Bmu, y = Cmx. Am =       −s1 −s2 ... −sn       , Bm =       1 1 . . . 1       Cm =

• r1 r2 · · · rn
• Corresponding transfer function:

G(s) =

n

• i=1

ri s + si Block diagram for the modal canonical form:

y u

• r1

rn −s1

• −s2

−sn r2 x1 x2 xn

ESAT–SCD–SISTA

CACSD

• pag. 62

✬ ✫ ✩ ✪

Similarity transformation

A state space model (A, B, C, D) is not unique for a phys- ical system. Let x = T ¯ x, det(T) = 0, ¯ A = T −1AT, ¯ B = T −1B, ¯ C = CT, ¯ D = D ⇒ ˙ ¯ x = T −1AT ¯ x + T −1Bu = ¯ A¯ x + ¯ Bu, y = CT ¯ x + Du = ¯ C¯ x + Du T is chosen to give the most convenient state-space descrip- tion for a given problem (e.g.control or modal canonical forms). Example: Eigenvalue decomposition of A: A = XΛX−1, Λ =         α1 β1 −β1 α1 ... λ1 ...         ˙ x = Ax + Bu y = Cx + Du

X−1x = z

= ⇒ ˙ z = Λz + X−1Bu y = CXz + Du

ESAT–SCD–SISTA

CACSD

• pag. 63

✬ ✫ ✩ ✪

Controllability

An n-th order system is called controllable if one can reach any state x from any given initial state x0 in a finite time. Controllability matrix: C =

• B AB · · · An−1B
• rank(C) = n ⇔ (A, B) controllable

Explanation: Consider a linear (discrete) system of the form xk+1 = Axk + Buk Then xk+1 = Ak+1x0 + AkBu0 + · · · + Buk such that xk+1 − Ak+1x0 =

• B AB · · · AkB
• 

     uk uk−1 . . . u0      

ESAT–SCD–SISTA

CACSD

• pag. 64

✬ ✫ ✩ ✪

Note that rank

• B AB · · · AkB
• = rank
• B AB · · · An−1B
• for k ≥ n − 1 (Cayley-Hamilton Theorem).

There exists always a vector       uk uk−1 . . . u0       if rank(C) = n. Remarks: Controllability matrices for a continuous time system and a discrete time system are of the same form.

ESAT–SCD–SISTA

CACSD

• pag. 65

✬ ✫ ✩ ✪

Observability

An n-th order system is called observable if knowledge of the input u and the output y over a finite time interval is sufficient to determine the state x. Observability matrix: O =       C CA . . . CAn−1       rank(O) = n ⇔ (A, C) observable Explanation: Consider an autonomous system xk+1 = Axk yk = Cxk Then we find easily yk = CAkx0

ESAT–SCD–SISTA

CACSD

• pag. 66

✬ ✫ ✩ ✪

and       y0 y1 . . . yk       =       C CA . . . CAk       x0 Note that rank       C CA · · · CAk       = rank       C CA · · · CAn−1       for k ≥ n − 1. One can always determine x0 if rank(O) = n. Remarks: Observability matrices for a continuous time sys- tem and a discrete time system are of the same form.

ESAT–SCD–SISTA

CACSD

• pag. 67

✬ ✫ ✩ ✪

The Popov-Belevitch-Hautus tests (PBH)

PBH controllability test: (A, B) is not controllable if and only if there exists a left eigenvector q = 0 of A such that ATq = qλ BTq = 0 In other words, (A, B) is controllable if and only if there is no left eigenvector q of A that is orthogonal to B. If there is a λ and q satisfying the PBH test, then we say that the mode corresponding to eigenvalue λ is uncon- trollable (uncontrollable mode), otherwise it is controllable (controllable mode). How to understand (for SISO) ? Put A in the modal canon- ical form: A =    ∗ λi ∗    , B =    ∗ ∗    , ⇒ q =    1    Uncontrollable mode: ˙ xi = λixi.

ESAT–SCD–SISTA

CACSD

• pag. 68

✬ ✫ ✩ ✪

PBH observability test: (A, C) is not observable if and only if there exists a right eigenvector p = 0 of A such that Ap = pλ Cp = 0 If there is a λ and p satisfy the PBH test, then we say that the mode corresponding to eigenvalue λ is unobservable (unobservable mode), otherwise it is observable (observable mode). How to understand for SISO ? Put A in a modal canonical form.

ESAT–SCD–SISTA

CACSD

• pag. 69

✬ ✫ ✩ ✪

Stability/Stabilizability/Detectability

Stability: A system with system matrix A is unstable if A has an eigenvalue λ with real(λ) ≥ 0 for a continuous time system and |λ| ≥ 1 for a discrete time system. Stabilizability: (A, B) is stabilizable if all unstable modes are controllable. Detectability: (A, C) is detectable if all unstable modes are observable.

ESAT–SCD–SISTA

CACSD

• pag. 70

✬ ✫ ✩ ✪

Example xk+1 = Axk + Buk yk = Cxk + Duk

• A B

C D

• =

        1.1 0 0 1 0 −0.5 0.6 0 2 0 −0.6 −0.5 0 −1 0 2 2 1 2 3 1         Modes Stable? Stabilizable? Detectable? 1.1 no yes no −0.5 ± 0.6i yes yes yes 2 no yes yes

ESAT–SCD–SISTA

CACSD

• pag. 71

✬ ✫ ✩ ✪

Kalman decomposition and minimal realization

Kalman decomposition Given a state space system [A, B, C, D]. Then we can always find an invertible similarity transformation T such that the transformed matrices have the structure TAT −1 =       r1 r2 r3 r4 r1 Aco A13 r2 A21 Aco A23 A24 r3 Aco r4 A43 Aco       TB =       m r1 Bco r2 Bco r3 r4       , CT −1 = r1 r2 r3 r4 p Cco Cco

• where

r1 = rank (OC), r2 = rank (C) − r1, r3 = rank (O) − r1, r4 = n − r1 − r2 − r3.

ESAT–SCD–SISTA

CACSD

• pag. 72

✬ ✫ ✩ ✪

The subsystem [ Aco, Bco, Cco, D ] is controllable and observable. The subsystem [

• Aco

A21 Aco

• ,
• Bco

Bco

• ,
• Cco 0
• , D ]

is controllable. The subsystem [

• Aco A13

Aco

• ,
• Bco
• , ( Cco Cco ), D ]

is observable. The subsystem [ Aco , 0, 0, D ] is neither controllable nor observable.

ESAT–SCD–SISTA

CACSD

• pag. 73

✬ ✫ ✩ ✪

Kalman decomposition

u y CO C ¯ O ¯ CO ¯ C ¯ O ESAT–SCD–SISTA

CACSD

• pag. 74

✬ ✫ ✩ ✪

Minimal realization: A minimal realization is one that has the smallest-size A matrix for all triplets [A, B, C] satisfying G(s) = D + C(sI − A)−1B where G(s) is a given transfer matrix.

• A B C D
• is minimal ⇔ controllable and observable.

Let

• Ai Bi Ci D
• , i = 1, 2, be two minimal realizations
• f a transfer matrix, then
• A1 B1 C1 D
• T ⇓ ⇑ T −1
• A2 B2 C2 D
• with

T = C1CT

2 (C2CT 2 )−1 ESAT–SCD–SISTA

CACSD

• pag. 75

✬ ✫ ✩ ✪

Input/output properties of state-space models Transfer matrix

General transfer matrix for a state space model : G(ξ) = D + C(ξI − A)−1B

ξ can be s (CT) or z (DT).

Poles

Characteristic polynomial of matrix A: α(ξ) = det(ξI − A) = αnξn + αn−1ξn−1 + · · · + α1ξ + α0 Characteristic equation: α(ξ) = 0 The eigenvalues λi, i = 1, · · · , n of the system matrix A are called the poles of the system. The pole polynomial is defined as π(ξ) = Πn

i=1(ξ − λi)

(= characteristic equation up to within a scalar)

ESAT–SCD–SISTA

CACSD

• pag. 76

✬ ✫ ✩ ✪

Physical interpretation of a pole :

Consider the following second order system :

• A B

C D

• =

   α β b1 −β α b2 c1 c2 0    The transfer matrix is given by (try to verify it) : G(ξ) = C(ξI − A)−1B + D = ξ(b1c1 + b2c2) + β(b2c1 − b1c2) − α(b1c1 + b2c2) ξ2 − 2αξ + α2 + β2 There are 2 poles at α ± jβ. There is a zero at α(b1c1 + b2c2) − β(b2c1 − b1c2) b1c1 + b2c2

ESAT–SCD–SISTA

CACSD

• pag. 77

✬ ✫ ✩ ✪

In continuous time the impulse response takes on the form: g(t) = L−1 {G(s)} = L−1 s(b1c1 + b2c2) + β(b2c1 − b1c2) − α(b1c1 + b2c2) s2 − 2αs + α2 + β2

• = (Am cos(βt) + Bm sin(βt))eαt = Cmeαt cos(βt + γ)

with Am = b1c1 + b2c2 and Bm = b2c1 − b1c2 ⇒ Cm =

• (b2

1 + b2 2)(c2 1 + c2 2)

⇒ γ = atan −Bm Am

• Define the damping ratio ζ and natural frequency ωn:

α = −ζωn, β = ωn

• 1 − ζ2

jω θ α + jβ α − jβ β α ζ = sin θ ωn stability degree r ESAT–SCD–SISTA

CACSD

• pag. 78

✬ ✫ ✩ ✪

For a second order system G(s) = ω2

n

s2 + 2ζωns + ω2

n

the time response looks like : rise time tr ≃ 1.8

ωn

settling time ts = 4.6

ζωn

peak time tp = π

ωd, ωd = ωn

• 1 − ζ2
• vershoot Mp = e−πζ/√

1−ζ2, 0 ≤ ζ < 1

In general : if a continuous-time system (A, B, C, 0) has poles at α ± jβ ⇒ the impulse response will have time modes of the form Ameαt cos(βt + γ) Am: amplitude γ: phase

• determined by B and C.

ESAT–SCD–SISTA

CACSD

• pag. 79

✬ ✫ ✩ ✪

In discrete time the impulse response matrix Gd(k) takes

• n the form :

Gd(k) = CdAk−1

d

Bd, k ≥ 1 which can be proven to be a sum of terms of the form : Cmbk−1 cos(ω(k − 1) + γ) each of which satisfies a second order linear system

• A B

C D

• =

   α β b1 −β α b2 c1 c2 0    Parameterize A as A =

• b cos ω

b sin ω −b sin ω b cos ω

• then

Ak = bk

• cos ωk

sin ωk − sin ωk cos ωk

• ESAT–SCD–SISTA

CACSD

• pag. 80

✬ ✫ ✩ ✪

Now, CAkB = (Am cos ωk + Bm sin ωk)bk = Cmbk cos(ωk + γ) with Am = Cm cos γ = b1c1 + b2c2 and Bm = −Cm sin γ = b2c1 − b1c2 ⇒ Cm =

• (b2

1 + b2 2)(c2 1 + c2 2)

⇒ b =

• α2 + β2

⇒ ω = atan β α

• ⇒ γ = atan

−Bm Am

• ESAT–SCD–SISTA

CACSD

• pag. 81

✬ ✫ ✩ ✪

Transmission zeros

Definition: The zeros of a LTI system are defined as those values ζ ∈ C for which the rank of the transfer matrix G(ξ) is lower than its normal rank (=rank of G(ξ) for almost all values of ξ): rank(G(ζ)) < normal rank Property: Let ζ be a zero of G(ξ) (p × m), then rank (G(ζ)) < normal rank, ⇓ if ζ is not a pole of G(ξ) ∃v = 0 s.t.

• D + C(ζI − A)−1B
• v = 0,

⇓ define w

∆

= (ζI − A)−1Bv Dv + Cw = 0, (ζI − A)w − Bv = 0, ⇓

• ζI − A −B

C D w v

• = 0

ESAT–SCD–SISTA

CACSD

• pag. 82

✬ ✫ ✩ ✪

How to find zeros for square MIMO systems (p = m) with invertible D ? (ζI − A)w − Bv = 0, Cw + Dv = 0 ⇓ Aw + Bv = wζ, v = −D−1Cw ⇓ (A − BD−1C)w = wζ ζ = eigenvalue of A − BD−1C, w = corresponding eigenvector ! For other cases :

• Generalized eigenvalue problem.
• Use Matlab function “tzero”.

Minimum and non-minimum phase systems: If a system has an unstable zero (in the right half–plane (RHP) or outside the unit circle), then it is a non-minimum phase (NMP) system, otherwise it is a minimum phase (MP) system.

ESAT–SCD–SISTA

CACSD

• pag. 83

✬ ✫ ✩ ✪

Physical explanation (continuous time) : Let ζ be a real zero, then there will exist vectors x0 and u0 such that:

• ζI − A −B

C D x0 u0

• = 0

This means if we have an input u0eζt, there exists an initial state x0 such that the response is y(t) = 0. For complex zeros : if ζ is a complex zero, then also its complex conjugate ζ∗ is a zero. Try to prove that if

• ζI − A −B

C D x0 u0

• = 0

the output y(t) will be exactly zero if the input is u(t) = |u0|eRe{ζ}t · cos(Im{ζ}t + φu0) and the initial state is chosen to be Re{x0}. “·” is the elementwise multiplication, both “cos()” and “e()” are as- sumed to be elementwise operators, Re{x0} is the real part

• f x0, Im{ζ} is the imaginary part of ζ and u0 = |u0|·ejφu0.

Try to derive equivalent formulas for the discrete-time case.

ESAT–SCD–SISTA

CACSD

• pag. 84

✬ ✫ ✩ ✪

Why zeros are important :

• Limited control system performance

k C G

If G is NMP, k can not go to ∞, since unstable zeros (which become unstable poles) put a limit on high gains.

• Stability of a stabilizing controller, Parity Interlacing

Property (PIP) : Let G be unstable. Then, G can be stabilized by a controller C which itself is stable ⇔ between every pair

• f real RHP zeros of G (including at ∞), there is an

even number of poles. Example: G(s) = s s2 − 1

×

• ×

−1 1 pole zero pole ∞

• zero
• ∞

G(s) can not be stabilized by a stable controller.

ESAT–SCD–SISTA

CACSD

• pag. 85

✬ ✫ ✩ ✪

RHP zeros can cause undershoot ! Let G(s) = Πn

i=1(1 − s ζi)

Πn+r

i=1 (1 − s λi),

where r = relative degree of G(s). Let y(t) be the step response. Fact :      y(0) and its first r − 1 derivative are 0; y(r) is the first non-zero derivative; y(∞) = G(0). The system has undershoot ⇔ the steady state value y(∞) has a sign opposite to y(r)(0), i.e. y(r)(0)y(∞) < 0. It can be proven that the system has undershoot ⇔ there is an odd number of RHP zeros.

ESAT–SCD–SISTA

CACSD

• pag. 86

✬ ✫ ✩ ✪

An example of state-space analysis

Example: Tape drive control - state space analysis Controllability/observability Direct rank test using SVD: The singular values of C are 4.5762, 4.3861, 1.2940, 1.2758, 0.5737, 0.4787, so the system is controllable (rank=6). The singular values of O are 2.5695, 1.4504, 0.7145, 0.7071, 0.5171, 0.2356, so the system is observable (rank=6).

ESAT–SCD–SISTA

CACSD

• pag. 87

✬ ✫ ✩ ✪

PBH test: There are 6 (independent) left eigenvectors associated with the 6 eigenvalues of A :            −0.1480 ± 0.0826i 0.0703 ± 0.5380i 0.1480 ∓ 0.0826i −0.0703 ∓ 0.5380i 0.2994 ± 0.2954i −0.2994 ∓ 0.2954i            ,            0.0766 0.5624 0.0766 0.5624 0.4218 0.4218                       0.4892 0.0000 0.4892 0.5105 0.5105            ,            −0.0046 −0.0227 0.0046 0.0227 −0.7067 0.7067            ,            0.0000 −0.0295 0.0000 −0.0295 −0.7065 −0.7065            It can be easily checked that none of them is orthogonal to

• B. The result is the same as the rank test above.

The same can be done to investigate the observability.

ESAT–SCD–SISTA

CACSD

• pag. 88

✬ ✫ ✩ ✪

Poles and zeros:

Poles: the eigenvalues of A are −0.2370 ± 0.5947i, 0, −0.2813, −0.9760, −0.9687 The system is not stable since there is a pole at zero. The first two oscillatory poles are from the spring-mass system consisting of the tape and the motor/capstan in-

• ertias. The zero pole reflects the pure integrative effect of

the tape on the capstans.

ESAT–SCD–SISTA

CACSD

• pag. 89

✬ ✫ ✩ ✪

Zeros: it can be checked in Matlab using function “tzero” that there is only one zero at -2. The eigenvector of the matrix

• −2I − A −B

C

• associated with the eigenvalue 0 is
• x0

u0

• with x0 =

           −0.1959 0.1959 0.1959 −0.1959 −0.4571 0.4571            , u0 =

• 0.4630

−0.4630

• .

Thus, let x(0) = x0, u(t) = u0e−2t, then the output y(t) will be always zero, ∀t ≥ 0.

ESAT–SCD–SISTA

CACSD

• pag. 90

✬ ✫ ✩ ✪

Physical explanation of the zero:

Since the initial values i1(0) = −i2(0) and the inputs e1(t) = −e2(t), and the motor drive systems in both sides are symmetric, the two capstans will rotate at the same speed, but in opposite directions. So the tape position over the read/write head will remain at 0 if the initial value is at 0. The tension in the tape is the sum of the tensions in the tape spring and the tape damping. The input voltages vary such a way that the tensions in the tape spring and the tape damping cancel each other out.

ESAT–SCD–SISTA

CACSD

• pag. 91

✬ ✫ ✩ ✪

Matlab Functions

c2d c2dm canon expm tf2ss ss2tf (d)impulse (d)step minreal tzero ctrb

• bsv

ctrbf

• bsvf

(d)lsim ss ssdata

ESAT–SCD–SISTA

CACSD

• pag. 92