Spring 2008 EEE 8007

Chapter #1

EEE 8007

Digital Control

ƒ State Space Representation of Discrete Time Control

Systems
ƒ State Feedback
ƒ Discrete Linear Quadratic Regulator (LQR)
ƒ Closed Loop Estimators

Chapter 1 1/14
Spring 2008 EEE 8007

State Space Representation of Discrete Time Control Systems

Most controllers are digital, hence we need to transform the continuous

systems to discrete.

Many methods for that

1) First Order Hold (exists in Matlab)

2) Tustin’s Bilinear Method (exists in Matlab)
3) Bilinear with pre-warping (exists in Matlab)
4) Zero pole matched (exists in Matlab)
5) Impulse (exists in Matlab)
6) Forward where Φ=I+AT, Γ=BT
7) Backward where Φ=(I-AT)-1BT, Γ=C(I-AT)-1, C=C(I-AT)-1, D=D+(I-
AT)-1 BT

Forward:

• X((k + 1)T ) − X(kT ) X((k + 1)T ) − X(kT )

X(kT ) = =
(k + 1)T − kT T

X((k + 1)T ) − X(kT )

Hence = f (X(kT ), U(kT ), kT )
T

or X((k + 1)T ) = Tf (X(kT ), U(kT ), kT ) + X(kT ) or

X((k + 1)T ) = g (X(kT ), U(kT ), kT )

Hence, for time-varying (linear or not) discrete time systems the state space
representation is:

Chapter 1 2/14
Spring 2008 EEE 8007

X(k + 1) = f (X(k ), U(k ), k )

.
Y(k ) = g (X(k ), U(k ), k )

If the system is linear and there is no direct coupling between the input and
the output then:

X(k + 1) = Φ(k ) X(k ) + Γ(k )U(k )

Y(k ) = C(k ) X(k )

And for LTI systems:

X(k + 1) = ΦX(k ) + ΓU(k )

Y(k ) = CX(k )

The relation of the matrices Φ and A, Γ and B can be found in any text book
and is:

A 2T 2 A 3T 3
Φ = e AT = I + AT + + + ...
2! 3!
T
Γ = ∫ e An dn B
0

In Matlab use: c2d

All the other concepts are the same:

TF of digital system: C[zI − Φ ]−1 Γ and hence poles of system: eig (Φ)

A=[-1 0;0 -2]; B=[1;0]; C=eye(2); D=0;

sys_c=ss(A,B,C,D);
sys_d=c2d(sys_c, 0.1);
initial(sys_c,[1,1]);
figure

Chapter 1 3/14
Spring 2008 EEE 8007

initial(sys_d,[1,1]);

Response to Initial Conditions

1

0.8
To: Out(1)
0.6

0.4

0.2
Amplitude

0
1

0.8
To: Out(2)

0.6

0.4

0.2

0
0 1 2 3 4 5 6
Time (sec)

Response to Initial Conditions

1

0.8
To: Out(1)

0.6

0.4

0.2
Amplitude

0
1

0.8
To: Out(2)

0.6

0.4

0.2

0
0 1 2 3 4 5 6
Time (sec)

Chapter 1 4/14
Spring 2008 EEE 8007

State Feedback

By following the same steps, as in continuous time systems:

•
X(k + 1) = Φ cl X(k ) + Γ cl U(k ) where Φcl=Φ-ΓK, Γcl=ΓF, Ccl=C-DK,
Dcl=DF. Thus we have described the closed loop system into the same form
as the open loop system.

Since the poles of the OL can be found by the |zI-Φ|=0, the close loop poles
are: |zI- (Φ-ΓK)|=0.

If we want the close loop poles at P=[P1 P2 … Pn]T, (n is the order of the
system) then:

Chapter 1 5/14
Spring 2008 EEE 8007

|zI-Φ + ΓK|=(z-P1)(z-P2)…(z-P3)Ù K=[K1 K2 …]T. Pole placement

Check CONTROLLABILITY

[
S = Γ ΦΓ Φ 2 Γ ... Φ n −1 Γ ]
Example:

clc
close all
clear all
A=[0 1 0;0 -1 1;0 0 -5];B=[0 0 5]';C=eye(3);D=0
sys_c=ss(A,B,C,D)
initial(sys_c,[1 1 1])

K=place(A,B,[-10 -5+2*j -5-2*j])

Acl=A-B*K;
sys_c_cl=ss(Acl,[],C,[])
figure
initial(sys_c_cl,[1 1 1])

Ts=0.01;
sys_d=c2d(sys_c,Ts)

figure
initial(sys_d,[1 1 1])

[F,G,C,D,Ts1]=ssdata(sys_d)
rank(ctrb(F,G))
K_d=place(F,G,[-0.5 -0.1+0.2*j -0.1-0.2*j])

Fcl=F-G*K_d;
figure
sys_d_cl=ss(Fcl,[],C,[])
initial(sys_d_cl,[1 1 1])

Chapter 1 6/14
Spring 2008 EEE 8007

Response to Initial Conditions

2.5

To: Out(1)
2

1.5

1
1
To: Out(2)
Amplitude

0.5

0
1
To: Out(3)

0.5

0
0 1 2 3 4 5 6 7 8 9 10
Time (sec)

Response to Initial Conditions

1.5
To: Out(1)

0.5

0
1
To: Out(2)

0
Amplitude

-1

-2
10
To: Out(3)

-10

-20
0 0.2 0.4 0.6 0.8 1 1.2
Time (sec)

Chapter 1 7/14
Spring 2008 EEE 8007

Response to Initial Conditions

2

To: Out(1) 1.5

1
1
To: Out(2)

0.9
Amplitude

0.8

0.7
1
To: Out(3)

0.5

0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Time (sec)

Response to Initial Conditions

2
To: Out(1)

-1
500
To: Out(2)
Amplitude

-500 4
x 10
10
To: Out(3)

-5
0 10 20 30 40 50 60
Time (sec)

Chapter 1 8/14
Spring 2008 EEE 8007

LQR Control

To use an LQR controller use the command:

>>dlqr()

Chapter 1 9/14
Spring 2008 EEE 8007

Discrete Time Estimators

The same method can be applied at the DTE:

ΓU X(k ) Y
Γ z −1I
X(k + 1)

Φ
ΦX

z −1I
Γ ~
X( k + 1) ~
ΓU ~
Y
~ X( k )
FX ( k )

End the error dynamics are: e(k + 1) = (Φ − G p C )e(k ) .

10% error of ICs:

F=[1 0.01 0; 0 0.99 0.0097; 0 0 0.9512];G=[0 0.002

0.0488]'; C=eye(1,3);D=0;
rank(obsv(F,C))
Gp=place(F',C',[-0.1 -0.1+0.1*j -0.1-0.1*j])' ;

Chapter 1 10/14
Spring 2008 EEE 8007

K_d=place(F,G,[-0.2 -0.2+0.2*j -0.2-0.2*j])

As it can be seen we estimate the X(k) from the estimation of X(k+1) after a
delay. To calculate the estimation of X(k+1) we use Y(k) and not Y(k+1).
So this estimator predicts X(k+1) from Y(k). For this reason is called
prediction Estimator (or priory estimator).

Chapter 1 11/14
Spring 2008 EEE 8007

Current Estimators:

Until now the estimated state vector X(k) was calculated by the one sample
delay of X(k+1). The X(k+1) was estimated by the measurements of Y(k).
What about if we use the sample Y(k+1) to calculate the X(k+1) and hence
the X(k). Then this estimator is called current estimator – or “posteriori”
estimator (used in Kalman Filter).

ˆ ( k ) and its relation with the

The current estimator will be symbolised by: X
~ ~
ˆ (k ) = X
X(k ) is: X
~
( k ) + Gc Y (k ) − CX( k ) ( )
Hence it is the same as the prediction plus some extra information. But it
must be noticed that for the calculation of the predictor the current is used:
~ ˆ ( k − 1) + ΓU ( k − 1) (OL estimation of current value)
X ( k ) = ΦX

The estimator is predicting the states and then it is using the current
information of Y(k) to find the correct the estimating value

Hence:

U ΓU + X(k ) Y
Γ z −1I C
+ X(k + 1)

Φ
ΦX Y

Plant

ΓU
+ +
z −1 I -
Γ ~ C
X(k + 1) ~
ΓU + ~ Y
~ X (k )
ΦX ( k ) +
Xˆ (k ) Gc
Φ
+
Estimator

Chapter 1 12/14
Spring 2008 EEE 8007

And the error dynamics are: e(k + 1) = (Φ − G c CΦ )e(k )

The previous example with a current estimator:

>> Gc=place(F',F'*C',[-0.1 -0.1+0.1*j -0.1-0.1*j])'

2.5

1.5

0.5

-0.5

-1

-1.5

-2
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

400

300

200

100

-100

-200

-300

-400
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

4
x 10
1.5

0.5

-0.5

-1

-1.5

-2
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

Chapter 1 13/14
Spring 2008 EEE 8007

1
K*u
z
Unit Delay1 Gain6 Scope1

Gain3

K*uve

K*u

Gain7

1
K*u
z
Unit Delay2 Gain5 Scope2

Gain1

K*uve

States

Chapter 1 14/14