Computer Implementation of Control Systems
Karl-Erik rzen, Anton Cervin
Session outline
Sampled-data control
Discretization of continuous-time controllers
Implementation of PID Controllers
Sampled-data control systems
y (t )
u (t )
t
u (t )
Process
y (t )
Hold
Sampler
yk
uk
uk
D-A
Computer
A-D
yk
Mix of continuous-time and discrete-time signals
Networked control systems
.
.
. ...
u(t)
y(t)
u(t)
y(t)
Process
t
DA
and
Hold
.
.
. ...
Sampler
and
AD
Communication network
uk
yk
uk
Computer
Extra delay, possibly lost packets
...
.
.
.
yk
Sampling
Computer
A/D
u
Algorithm
D/A
y
Process
AD-converter acts as sampler
A/D
DA-converter acts as a hold device
Normally, zero-order-hold is used piecewise constant control
signals
Aliasing
1
0
1
0
s =
N =
2
h
s
2
5
Time
10
= sampling frequency
= Nyquist frequency
Frequencies above the Nyquist frequency are folded and appear as
low-frequency signals.
The fundamental alias for a frequency f1 is given by
f = ( f1 + f N ) mod ( fs ) f N
Above: f1 = 0.9, fs = 1, f N = 0.5, f = 0.1
Anti-aliasing filter
Analog low-pass filter that eliminates all frequencies above the
Nyquist frequency
Analog filter
2-6th order Bessel or Butterworth filter
Difficulties with changing h (sampling interval)
Analog + digital filter
Fixed, fast sampling with fixed analog filter
Downsampling using digital LP-filter
Control algorithm at the lower rate
Easy to change sampling interval
The filter may have to be included in the control design
Example Prefiltering
(a)
(b)
1
0
1
0
(c)
10
20
30
(d)
1
10
20
30
Time
d = 0.9, N = 0.5, alias = 0.1
6th order Bessel with B = 0.25
10
20
30
10
20
30
Time
Design approaches
Digital controllers can be designed in two different ways:
Discrete-time design sampled control theory
Sample the continuous system
Design a digital controller for the sampled system
Z-transform domain
state-space domain
Continuous time design + discretization
Design a continuous controller for the continuous system
Approximate the continuous design
Use fast sampling
Disk drive example
Control of the arm of a disk drive
k
G (s) =
Js2
Continuous time controller
s+b
bK
Uc (s) K
Y (s)
U (s) =
a
s+a
Discrete time controller (continuous time design + discretization)
u(t k) = K ( ab uc (t k) y(t k) + x(t k))
x(t k + h) = x(t k) + h ((a b) y(t k) ax(t k))
Disk drive example
y: = adin(in2)
u:=K*(b/a*uc-y+x)
dout(u)
x:=x+h*((a-b)*y-a*x)
Clock
Algorithm
Output
Sampling period h = 0.2/ 0
1
10
5
Time (0t)
10
Input
0.5
0
0.5
Increased sampling period
a) h = 0.5/ 0 b) h = 1.08/ 0
(b)
Output
Output
(a)
1
10
10
5
10
Time (0t)
15
0.5
Input
0.5
Input
0
0.5
0
0.5
5
10
Time (0t)
15
Better performance?
Dead-beat control h = 1.4/ 0
Position
u(t k) = t0 uc (t k) + t1 uc (t k1 ) s0 y(t k) s1 y(t k1 ) r1 u(t k1 )
1
Velocity
10
10
5
Time (0t)
10
0.5
Input
0.5
0
0.5
However, long sampling periods also have problems
open loop between the samples
disturbance and reference changes that occur between samples will remain undetected until the next sample
Sampled control theory
Computer
Clock
A-D
{y(t k )}
{u(t k )}
Algorithm
u(t)
D-A
y(t )
Process
Basic idea: look at the sampling instances only
System theory analogous to continuous-time systems
Better performance can be achieved
Potential problem with intersample behaviour
Sampling of systems
Look at the system from the point of view of the computer
Clock
{u(t k )}
y(t)
u(t)
D-A
System
A-D
{ y (tk )}
Zero-order-hold sampling of a system
Let the inputs be piecewise constant
Look at the sampling points only
Solve the system equation
Sampling a continuous-time system
Process:
dx
= Ax(t) + Bu(t)
dt
y(t) = Cx(t) + Du(t)
Solve the system equation:
t
A(ttk )
A(ts )
x ( t) = e
x(t k) +
e
Bu(s ) ds
=e
A(ttk )
x(t k) +
= eA(ttk) x(tk) +
tk
t
A(ts )
tk
ttk
ds Bu(t k) (u const.)
eAs ds Bu(t k) (variable change)
= (t, tk) x(tk) + (t, tk)u(tk)
Periodic sampling
Assume periodic sampling, i.e. tk = kh. Then
x( kh + h) = x( kh) + u( kh)
y( kh) = Cx( kh) + Du( kh)
where
= eAh
h
=
eAs ds B
0
Time-invariant linear system!
Example: Sampling of inverted pendulum
0 1
0
dx
x+
u
dt
1 0
1
y = 1 0x
We get
cosh h sinh h
Ah
=e =
sinh h cosh h
h
sinh s
cosh h 1
ds =
cosh s
sinh h
0
Several ways to calculate and . Matlab
Sampling a system with a time delay
Sampling the system
dx(t)
dt
= Ax(t) + Bu(t ),
we get the discrete-time system
x( kh + h) = x( kh) + 0 u( kh) + 1 u( kh h)
where
= eAh
h
0 =
eAs ds B
0
1 = eA(h )
eAs ds B
We get one extra state u( kh h) in the sampled system
Stability region
In continuous time the stability region is the complex left half
plane, i.e., the system is stable if all the poles are in the left
half plane.
In discrete time the stability region is the unit circle.
Digital control design
Similar to analog control design, but
Z-transform instead of Laplace transform
zX ( z) x(t k+1 )
z1 X ( z) x(t k1 )
Poles are placed within the unit circle
The frequency response is more difficult to compute
The sampling interval h is a new design parameter
Choice of sampling interval
Nyquists sampling theorem:
We must sample at least twice as fast as the highest
frequency we are interested in
What frequencies are we interested in?
Typical loop transfer function L(i ) = P(i ) C(i ):
1
Frstrkning
10
10
10
10
10
10
50
Fas
100
150
200
250
1
10
10
Frekvens [rad/s]
c = cross-over frequency, m = phase margin
We should have s 2 c
Sampling interval rule of thumb
A sample-and-hold (S&H) circuit can be approximated by a delay of
h/2.
G S&H (s) esh/2
This will decrease the phase margin by
arg G S&H (i c ) = arg ei c h/2 = c h/2
Assume we can accept a phase loss between 5 and 15 . Then
0.15 < c h < 0.5
This corresponds to a Nyquist frequency about 6 to 20 times larger
than the crossover frequency
Example: control of inverted pendulum
h = 0.1,
c h = 0.28
h = 0.3,
c h = 0.78
h = 0.5,
c h = 1.12
10
10
10
10
10
10
20
5
Time
20
5
Time
20
Time
Large c h may seem OK, but beware!
Digital design assuming perfect model
Controller perfectly synchronized with initial disturbance
Pendulum with non-synchronized disturbance
h = 0.1,
c h = 0.28
h = 0.3,
c h = 0.78
h = 0.5,
c h = 1.12
10
10
10
10
10
10
20
5
Time
20
5
Time
20
Time
Accounting for the anti-aliasing filter
Assume we also have a second-order Butterworth anti-aliasing
filter with a gain of 0.1 at the Nyquist frequency. The filter gives an
additional phase margin loss of 1.4 c h.
Again assume we can accept a phase loss of 5 to 15 . Then
0.05 < c h < 0.14
This corresponds to a Nyquist frequency about 23 to 70 times
larger than the crossover frequency
Session outline
Sampled-data control
Discretization of continuous-time controllers
Implementation of PID Controllers
Discretization of continuous-time controllers
Basic idea: Reuse the analog design
H ( z) G (s)
u(t)
A-D
{u ( kh )}
Algorithm
{ y ( kh )}
y (t)
D-A
Clock
Want to get:
A/D + Algorithm + D/A G (s)
Methods:
Approximate s, i.e., H ( z) = G (s )
Other discretization methods (Matlab)
Approximation methods
Forward Difference (Eulers method):
dx(t)
x(t k+1 ) x(t k)
dt
h
s =
z1
h
Backward Difference:
dx(t)
x(t k) x(t k1 )
dt
h
s =
Tustin:
dx(t)
dt
z1
zh
tk+1 )
+ dx(dt
x(t k+1 ) x(t k)
2
h
1
s = h2 zz
+1
Stability of approximations
How is the continuous-time stability region (left half plane) mapped?
Forward differences
Backward differences
Tustin
Discretization example
Controller designed in continuous-time:
b
E(s)
U (s) =
s+a
Discretization using Forward Euler (s =
z1
):
h
b
e( k)
u( k) =
( z 1)/h + a
( z 1 + ha)u( k) = bhe( k)
u( k + 1) = (1 ha)u( k) + bhe( k)
u( k) = (1 ha)u( k 1) + bhe( k 1)
Controller stable if 1 < (1 ha) < 1, i.e., 0 < h < 2/a (does not 3
imply that the closed loop system is stable, though)
Controller Synthesis
Process Model
G(s)
x = Ax + Bu
y = Cx + Du
Control Design in Continuous-Time
Loop shaping
Pole placement
PID
.
Discretize the Controller
Euler
Tustin
.
Discretize the process
e.g. ZOH Sampling
x(k + 1) = x(k) + u(k)
y(k) = Cx(k) + Du(k)
Control Design in Discrete-Time
Pole placement
LQG
.
Difference Equation
Software algorithm
Session outline
Sampled-data control
Discretization of continuous-time controllers
Implementation of PID Controllers
PID Algorithm
Textbook Algorithm:
u(t)
K ( e(t)
U (s) = K ( E(s) +
1
TI
1
sTI
e( )d +
TD dedt(t) )
E(s)
+ TD sE(s))
Algorithm Modifications
Modifications are needed to make the controller practically useful
Limitations of derivative gain
Derivative weighting
Setpoint weighting
Limitations of derivative gain
We do not want to apply derivation to high frequency measurement
noise, therefore the following modification is used:
sTD
sTD
1 + sTD / N
N = maximum derivative gain, often 10 20
Derivative weighting
The setpoint is often constant for long periods of time
Setpoint often changed in steps D-part becomes very large.
Derivative part applied on part of the setpoint or only on the measurement signal.
sTD
D (s) =
( Ysp(s) Y (s))
1 + sTD / N
Often, = 0 in process control, = 1 in servo control
Setpoint weighting
An advantage to also use weighting on the setpoint.
u = K ( ysp y)
replaced by
u = K ( ysp y)
0 1
A way of introducing feedforward from the reference signal (position
a closed loop zero)
Improved set-point responses.
A better algorithm
1
TD s
Y (s))
U (s) = K ( yr y +
E(s)
sTI
1 + sTD / N
Modifications:
Setpoint weighting ( ) in the proportional term improves setpoint response
Limitation of the derivative gain (low-pass filter) to avoid derivation of measurement noise
Derivative action only on y to avoid bumps for step changes in
the reference signal
Control Signal Limitations
All actuators saturate.
Problems for controllers with integration.
When the control signal saturates the integral part will continue to grow
integrator (reset) windup.
When the control signal saturates the integral part will integrate up to a
very large value. This may cause large overshoots.
2 Output y and yref
1.5
1
0.5
0
0
10
20
Control variable u
0.2
0
0.2
4
0
10
20
Anti-Reset Windup
Several solutions exist:
controllers on velocity form (not discussed here))
limit the setpoint variations (saturation never reached)
conditional integration (integration is switched off when the
control is far from the steady-state)
tracking (back-calculation)
Tracking
when the control signal saturates, the integral is recomputed so
that its new value gives a control signal at the saturation limit
to avoid resetting the integral due to, e.g., measurement noise,
the re-computation is done dynamically, i.e., through a LP-filter
with a time constant Tr .
Tracking
y
KTds
Actuator
e = r y
K
Ti
1
s
es
1
Tt
y
KT d s
e = r y
Actuator
model
K
K
Ti
Actuator
1
s
1
Tt
+
es
Tracking
1
0.5
0
10
20
30
10
20
30
10
20
30
0.15
0.05
0.05
0
0.4
0.8
Discretization
P-part:
u P ( k) = K ( ysp( k) y( k))
Discretization
I-part:
K
I ( t) =
TI
t
0
K
dI
e( )d ,
=
e
dt
TI
Forward difference
I (t k+1 ) I (t k)
K
=
e(t k)
h
TI
I(k+1) := I(k) + (K*h/Ti)*e(k)
The I-part can be precalculated in UpdateStates
Backward difference
The I-part cannot be precalculated, i(k) = f(e(k))
Others
Discretization
D-part (assume = 0):
sTD
D=K
( Y (s))
1 + sTD / N
TD dD
dy
+ D = K TD
N dt
dt
Forward difference (unstable for small TD )
Discretization, cont.
D-part:
Backward difference
TD D (t k) D (t k1 )
+ D (t k)
N
h
y(t k) y(t k1 )
= K TD
h
TD
D (t k1 )
D (t k) =
TD + Nh
K TD N
( y(tk) y(tk1 ))
TD + Nh
Discretization
Tracking:
v := P + I + D;
u := sat(v,umax,umin);
I := I + (K*h/Ti)*e + (h/Tr)*(u - v);
PID code
PID-controller with anti-reset windup
y = yIn.get(); // A-D conversion
e = yref - y;
D = ad * D - bd * (y - yold);
v = K*(beta*yref - y) + I + D;
u = sat(v,umax,umin)}
uOut.put(u);
// D-A conversion
I = I + (K*h/Ti)*e + (h/Tr)*(u - v);
yold = y
ad and bd are precalculated parameters given by the backward
difference approximation of the D-term.
Execution time for CalculateOutput can be minimized even further.
Alternative PID Implementation
The PID controller described so far has constant gain, K (1 + N ), a
high frequencies, i.e., no roll-off at high frequencies.
An alternative is to instead of having a low pass filter only on the
derivative part use a second-order low-pass filter on the measured
signal before it enters a PID controller with a pure derivative part.
1
Y (s)
Y f (s) = 2 2
T f s + 1.4T f s + 1
1
U (s) = K ( Yre f (s) Y f (s) +
( Yre f (s) Y f (s)) TD sY f (s))
TI s
Class SimplePID
public class SimplePID {
private double u,e,v,y;
private double K,Ti,Td,Beta,Tr,N,h;
private double ad,bd;
private double D,I,yOld;
public SimplePID(double nK, double nTi, double NTd,
double nBeta, double nTr, double nN, double nh) {
updateParameters(nK,nTi,nTd,nBeta,nTr,nN,nh);
}
public void updateParameters(double nK, double nTi, double NTd,
double nBeta, double nTr, double nN, double nh) {
K = nK; Ti = nTi; Td = nTd; Beta = nBeta;
Tr = nTr
N = nN;
h = nh;
ad = Td / (Td + N*h);
bd = K*ad*N;
}
public double calculateOutput(double yref, double newY) {
y = newY;
e = yref - y;
D = ad*D - bd*(y - yOld);
v = K*(Beta*yref - y) + I + D;
return v;
}
public void updateState(double u) {
I = I + (K*h/Ti)*e + (h/Tr)*(u - v);
yOld = y;
}
}
Extract from Regul
public class Regul extends Thread {
private SimplePID pid;
public Regul() {
pid = new SimplePID(1,10,0,1,10,5,0.1);
}
public void run() {
// Other stuff
while (true) {
y = getY();
yref = getYref():
u = pid.calculateOutput(yref,y);
u = limit(u);
setU(u);
pid.updateState(u);
// Timing Code
}
}
}
Task Models
Further reading
B. Wittenmark, K. J. strm, K.-E. rzn: Computer Control:
An Overview. IFAC Professional Brief, 2002.
(Available at http://www.control.lth.se)
K. J. strm, B. Wittenmark: Computer-Controlled Systems,
Third Ed. Prentice Hall, 1997.
K. J. strm, Tore Hgglund: Advanced PID Control. The
Instrumentation, Systems, and Automation Society, 2005.
