Advanced Control using

M ATLAB

or Stabilising the
unstabilisable
David I. Wilson
Auckland University of Technology
New Zealand
February 12, 2014
Copyright © 2014 David I. Wilson
Auckland University of Technology
New Zealand

Creation date: February, 2014.

All rights reserved. No part of this work may be reproduced, stored in a retrieval system, or
transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or
otherwise, without prior permission.
C ONTENTS

1 Introduction 1
1.1 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Matlab for computer aided control design . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Alternative computer design aids . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Economics of control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Laboratory equipment for control tests . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.1 Plants with one input and one output . . . . . . . . . . . . . . . . . . . . . . 6
1.4.2 Multi-input and multi-output plants . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Slowing down Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 From differential to difference equations 13

2.1 Computer in the loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.1 Sampling an analogue signal . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1.2 Selecting a sample rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.3 The sampling theorem and aliases . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.4 Discrete frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Finite difference models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.1 Difference equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 The z transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1 z-transforms of common functions . . . . . . . . . . . . . . . . . . . . . . . . 23
2.4 Inversion of z-transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4.1 Inverting z-transforms with symbolically . . . . . . . . . . . . . . . . . . . . 24
2.4.2 The partial fraction method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4.3 Long division . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.4 Computational approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4.5 Numerically inverting the Laplace transform . . . . . . . . . . . . . . . . . . 31
2.5 Discretising with a sample and hold . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.5.1 Converting Laplace transforms to z-transforms . . . . . . . . . . . . . . . . 36
2.5.2 The bilinear transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.6 Discrete root locus diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.7 Multivariable control and state space analysis . . . . . . . . . . . . . . . . . . . . . . 44
2.7.1 States and state space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.7.2 Converting differential equations to state-space form . . . . . . . . . . . . . 47
2.7.3 Interconverting between state space and transfer functions . . . . . . . . . . 49
2.7.4 Similarity transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.7.5 Interconverting between transfer functions forms . . . . . . . . . . . . . . . 54
2.7.6 The steady state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.8 Solving the vector differential equation . . . . . . . . . . . . . . . . . . . . . . . . . 57
2.8.1 Numerically computing the discrete transformation . . . . . . . . . . . . . . 61
2.8.2 Using M ATLAB to discretise systems . . . . . . . . . . . . . . . . . . . . . . . 63
2.8.3 Time delay in state space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
2.9 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

iii
iv CONTENTS

2.9.1 Stability in the continuous domain . . . . . . . . . . . . . . . . . . . . . . . . 69

2.9.2 Stability of the closed loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
2.9.3 Stability of discrete time systems . . . . . . . . . . . . . . . . . . . . . . . . . 73
2.9.4 Stability of nonlinear differential equations . . . . . . . . . . . . . . . . . . . 74
2.9.5 Expressing matrix equations succinctly using Kronecker products . . . . . . 80
2.9.6 Summary of stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
2.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

3 Modelling dynamic systems with differential equations 85

3.1 Dynamic system models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
3.1.1 Steady state and dynamic models . . . . . . . . . . . . . . . . . . . . . . . . 86
3.2 A collection of illustrative models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.2.1 Simple models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.3 Chemical process models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
3.3.1 A continuously-stirred tank reactor . . . . . . . . . . . . . . . . . . . . . . . 92
3.3.2 A forced circulation evaporator . . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.3.3 A binary distillation column model . . . . . . . . . . . . . . . . . . . . . . . 95
3.3.4 Interaction and the Relative Gain Array . . . . . . . . . . . . . . . . . . . . . 103
3.4 Curve fitting experimental data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
3.4.1 Polynomial regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
3.4.2 Least squares improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
3.4.3 Nonlinear least-squares model identification . . . . . . . . . . . . . . . . . . 112
3.4.4 Parameter confidence intervals . . . . . . . . . . . . . . . . . . . . . . . . . . 116
3.5 Numerical tools for modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
3.5.1 The solution of nonlinear ordinary differential equations . . . . . . . . . . . 120
3.5.2 Differential/Algebraic equation systems and algebraic loops . . . . . . . . . 123
3.6 Linearisation of nonlinear dynamic equations . . . . . . . . . . . . . . . . . . . . . . 125
3.6.1 Linearising a nonlinear tank model . . . . . . . . . . . . . . . . . . . . . . . 127
3.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

4 The PID controller 131

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
4.1.1 P, PI or PID control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
4.2 The industrial PID algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
4.2.1 Implementing the derivative component . . . . . . . . . . . . . . . . . . . . 133
4.2.2 Variations of the PID algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 134
4.2.3 Integral only control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.3 Simulating a PID process in S IMULINK . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.4 Extensions to the PID algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.4.1 Avoiding derivative kick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.4.2 Input saturation and integral windup . . . . . . . . . . . . . . . . . . . . . . 140
4.5 Discrete PID controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
4.5.1 Discretising continuous PID controllers . . . . . . . . . . . . . . . . . . . . . 144
4.5.2 Simulating a PID controlled response in Matlab . . . . . . . . . . . . . . . . 146
4.5.3 Controller performance as a function of sample time . . . . . . . . . . . . . 148
4.6 PID tuning methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
4.6.1 Open loop tuning methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
4.6.2 Closed loop tuning methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
4.6.3 Closed loop single-test tuning methods . . . . . . . . . . . . . . . . . . . . . 159
4.6.4 Summary on closed loop tuning schemes . . . . . . . . . . . . . . . . . . . . 166
4.7 Automated tuning by relay feedback . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
4.7.1 Describing functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
4.7.2 An example of relay tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
4.7.3 Self-tuning with noise disturbances . . . . . . . . . . . . . . . . . . . . . . . 173
CONTENTS v

4.7.4 Modifications to the relay feedback estimation algorithm . . . . . . . . . . . 176

4.8 Drawbacks with PID controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
4.8.1 Inverse response processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
4.8.2 Approximating inverse-response systems with additional deadtime . . . . 184
4.9 Dead time compensation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
4.10 Tuning and sensitivity of control loops . . . . . . . . . . . . . . . . . . . . . . . . . . 188
4.11 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191

5 Digital filtering and smoothing 193

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
5.1.1 The nature of industrial noise . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
5.1.2 Differentiating without smoothing . . . . . . . . . . . . . . . . . . . . . . . . 196
5.2 Smoothing measured data using analogue filters . . . . . . . . . . . . . . . . . . . . 197
5.2.1 A smoothing application to find the peaks and troughs . . . . . . . . . . . . 197
5.2.2 Filter types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197
5.2.3 Classical analogue filter families . . . . . . . . . . . . . . . . . . . . . . . . . 201
5.3 Discrete filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209
5.3.1 A low-pass filtering application . . . . . . . . . . . . . . . . . . . . . . . . . . 209
5.3.2 Digital filter approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
5.3.3 Efficient hardware implementation of discrete filters . . . . . . . . . . . . . 214
5.3.4 Numerical and quantisation effects for high-order filters . . . . . . . . . . . 217
5.4 The Fourier transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219
5.4.1 Fourier transform definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
5.4.2 Orthogonality and frequency spotting . . . . . . . . . . . . . . . . . . . . . . 224
5.4.3 Using M ATLAB’s FFT function . . . . . . . . . . . . . . . . . . . . . . . . . . 225
5.4.4 Periodogram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 227
5.4.5 Fourier smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228
5.5 Numerically differentiating industrial data . . . . . . . . . . . . . . . . . . . . . . . 230
5.5.1 Establishing feedrates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 231
5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233

6 Identification of process models 235

6.1 The importance of system identification . . . . . . . . . . . . . . . . . . . . . . . . . 235
6.1.1 Basic definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236
6.1.2 Black, white and grey box models . . . . . . . . . . . . . . . . . . . . . . . . 237
6.1.3 Techniques for identification . . . . . . . . . . . . . . . . . . . . . . . . . . . 238
6.2 Graphical and non-parametric model identification . . . . . . . . . . . . . . . . . . 239
6.2.1 Time domain identification using graphical techniques . . . . . . . . . . . . 239
6.2.2 Experimental frequency response analysis . . . . . . . . . . . . . . . . . . . 245
6.2.3 An alternative empirical transfer function estimate . . . . . . . . . . . . . . 252
6.3 Continuous model identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254
6.3.1 Fitting transfer functions using nonlinear least-squares . . . . . . . . . . . . 254
6.3.2 Identification of continuous models using successive differentiation . . . . 256
6.3.3 Practical continuous model identification . . . . . . . . . . . . . . . . . . . . 259
6.4 Popular discrete-time linear models . . . . . . . . . . . . . . . . . . . . . . . . . . . 260
6.4.1 Extending the linear model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262
6.4.2 Output error model structures . . . . . . . . . . . . . . . . . . . . . . . . . . 264
6.4.3 General input/output models . . . . . . . . . . . . . . . . . . . . . . . . . . . 265
6.5 Regressing discrete model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 266
6.5.1 Simple offline system identification routines . . . . . . . . . . . . . . . . . . 269
6.5.2 Bias in the parameter estimates . . . . . . . . . . . . . . . . . . . . . . . . . . 270
6.5.3 Using the System Identification toolbox . . . . . . . . . . . . . . . . . . . . . 271
6.5.4 Fitting parameters to state space models . . . . . . . . . . . . . . . . . . . . 275
6.6 Model structure determination and validation . . . . . . . . . . . . . . . . . . . . . 277
vi CONTENTS

6.6.1 Estimating model order and deadtime . . . . . . . . . . . . . . . . . . . . . . 280

6.6.2 Robust model fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282
6.6.3 Common nonlinear model structures . . . . . . . . . . . . . . . . . . . . . . 283
6.7 Online model identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283
6.7.1 Recursive least squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284
6.7.2 Recursive least-squares in M ATLAB . . . . . . . . . . . . . . . . . . . . . . . 289
6.7.3 Tracking the precision of the estimates . . . . . . . . . . . . . . . . . . . . . . 293
6.8 The forgetting factor and covariance windup . . . . . . . . . . . . . . . . . . . . . . 295
6.8.1 The influence of the forgetting factor . . . . . . . . . . . . . . . . . . . . . . . 297
6.8.2 Covariance wind-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 298
6.9 Identification by parameter optimisation . . . . . . . . . . . . . . . . . . . . . . . . . 300
6.10 Online estimating of noise models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304
6.10.1 A recursive extended least-squares example . . . . . . . . . . . . . . . . . . 306
6.10.2 Recursive identification using the SI toolbox . . . . . . . . . . . . . . . . . . 309
6.10.3 Simplified RLS algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 310
6.11 Closed loop identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312
6.11.1 Closed loop RLS in Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . 313
6.12 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313

7 Adaptive Control 321

7.1 Why adapt? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 321
7.1.1 The adaption scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322
7.1.2 Classification of adaptive controllers . . . . . . . . . . . . . . . . . . . . . . . 323
7.2 Gain scheduling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323
7.3 The importance of identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325
7.3.1 Polynomial manipulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325
7.4 Self tuning regulators (STRs) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 326
7.4.1 Simple minimum variance control . . . . . . . . . . . . . . . . . . . . . . . . 327
7.5 Adaptive pole-placement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 329
7.5.1 The Diophantine equation and the closed loop . . . . . . . . . . . . . . . . . 330
7.5.2 Solving the Diophantine equation in Matlab . . . . . . . . . . . . . . . . . . 331
7.5.3 Adaptive pole-placement with identification . . . . . . . . . . . . . . . . . . 334
7.5.4 Adaptive pole-placement in S IMULINK . . . . . . . . . . . . . . . . . . . . . 339
7.6 Practical adaptive pole-placement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 340
7.6.1 Dealing with non-minimum phase systems . . . . . . . . . . . . . . . . . . . 340
7.6.2 Separating stable and unstable factors . . . . . . . . . . . . . . . . . . . . . . 343
7.6.3 Experimental adaptive pole-placement . . . . . . . . . . . . . . . . . . . . . 346
7.6.4 Minimum variance control with dead time . . . . . . . . . . . . . . . . . . . 346
7.7 Summary of adaptive control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352

8 Multivariable controller design 357

8.1 Controllability and observability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 358
8.1.1 Controllability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 358
8.1.2 Observability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 360
8.1.3 Computing controllability and observability . . . . . . . . . . . . . . . . . . 361
8.1.4 State reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362
8.2 State space pole-placement controller design . . . . . . . . . . . . . . . . . . . . . . 365
8.2.1 Poles and where to place them . . . . . . . . . . . . . . . . . . . . . . . . . . 369
8.2.2 Deadbeat control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 371
8.3 Estimating the unmeasured states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373
8.4 Combining estimation and state feedback . . . . . . . . . . . . . . . . . . . . . . . . 374
8.5 Generic model control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 376
8.5.1 The tuning parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 378
8.5.2 GMC control of a linear model . . . . . . . . . . . . . . . . . . . . . . . . . . 380
CONTENTS vii

8.5.3 GMC applied to a nonlinear plant . . . . . . . . . . . . . . . . . . . . . . . . 382

8.6 Exact feedback linearisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 386
8.6.1 The nonlinear system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 387
8.6.2 The input/output feedback linearisation control law . . . . . . . . . . . . . 387
8.6.3 Exact feedback example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 389
8.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393

9 Classical optimal control 395

9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 395
9.2 Parametric optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 396
9.2.1 Choosing a performance indicator . . . . . . . . . . . . . . . . . . . . . . . . 396
9.2.2 Optimal tuning of a PID regulator . . . . . . . . . . . . . . . . . . . . . . . . 397
9.2.3 Using S IMULINK inside an optimiser . . . . . . . . . . . . . . . . . . . . . . . 403
9.2.4 An optimal batch reactor temperature policy . . . . . . . . . . . . . . . . . . 404
9.3 The general optimal control problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 405
9.3.1 The optimal control formulation . . . . . . . . . . . . . . . . . . . . . . . . . 407
9.3.2 The two-point boundary problem . . . . . . . . . . . . . . . . . . . . . . . . 409
9.3.3 Optimal control examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 410
9.3.4 Problems with a specified target set . . . . . . . . . . . . . . . . . . . . . . . 414
9.4 Linear quadratic control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 417
9.4.1 Continuous linear quadratic regulators . . . . . . . . . . . . . . . . . . . . . 417
9.4.2 Analytical solution to the LQR problem . . . . . . . . . . . . . . . . . . . . . 419
9.4.3 The steady-state solution to the matrix Riccati equation . . . . . . . . . . . . 423
9.4.4 The discrete LQR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 426
9.4.5 A numerical validation of the optimality of LQR . . . . . . . . . . . . . . . . 431
9.4.6 An LQR with integral states . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435
9.5 Estimation of state variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 439
9.5.1 Random processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 439
9.5.2 Combining deterministic and stochastic processes . . . . . . . . . . . . . . . 445
9.5.3 The Kalman filter estimation scheme . . . . . . . . . . . . . . . . . . . . . . . 446
9.5.4 The steady-state form of the Kalman filter . . . . . . . . . . . . . . . . . . . . 450
9.5.5 Current and future prediction forms . . . . . . . . . . . . . . . . . . . . . . . 452
9.5.6 An application of the Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . 455
9.5.7 The role of the Q and R noise covariance matrices in the state estimator . . 456
9.5.8 Extensions to the basic Kalman filter algorithm . . . . . . . . . . . . . . . . . 460
9.5.9 The Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463
9.5.10 Combining state estimation and state feedback . . . . . . . . . . . . . . . . . 465
9.5.11 Optimal control using only measured outputs . . . . . . . . . . . . . . . . . 466
9.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 467

10 Predictive control 469

10.1 Model predictive control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 469
10.1.1 Constrained predictive control . . . . . . . . . . . . . . . . . . . . . . . . . . 472
10.1.2 Dynamic matrix control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 476
10.2 A Model Predictive Control Toolbox . . . . . . . . . . . . . . . . . . . . . . . . . . . 482
10.2.1 A model predictive control GUI . . . . . . . . . . . . . . . . . . . . . . . . . 482
10.2.2 MPC toolbox in M ATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 482
10.2.3 Using the MPC toolbox in S IMULINK . . . . . . . . . . . . . . . . . . . . . . 484
10.2.4 MPC control of an unstable 3 degree of freedom helicopter . . . . . . . . . . 485
10.2.5 Further readings on MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 489
10.3 Optimal control using linear programming . . . . . . . . . . . . . . . . . . . . . . . 489
10.3.1 Development of the LP problem . . . . . . . . . . . . . . . . . . . . . . . . . 490
viii CONTENTS

11 Expert systems and neural networks 499

11.1 Expert systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 499
11.1.1 Where are they used? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 500
11.1.2 Features of an expert system . . . . . . . . . . . . . . . . . . . . . . . . . . . 501
11.1.3 The user interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 502
11.1.4 Expert systems used in process control . . . . . . . . . . . . . . . . . . . . . 503
11.2 Neural networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 505
11.2.1 The architecture of the neural network . . . . . . . . . . . . . . . . . . . . . . 507
11.2.2 Curve fitting using neural networks . . . . . . . . . . . . . . . . . . . . . . . 511
11.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 515

A List of symbols 517

B Useful utility functions in Matlab 519

C Transform pairs 521

D A comparison of Maple and MuPad 523

D.1 Partial fractions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523
D.2 Integral transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523
D.3 Differential equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 524
D.4 Vectors and matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 525

E Useful test models 527

E.1 A forced circulation evaporator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 527
E.2 Aircraft model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 528
L IST OF F IGURES

1.1 Traditional vs. Advanced control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2 Economic improvements of better control . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Blackbox configuration. The manual switch marked will toggle between either 7
or 9 low-pass filters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 The “Black-box” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.5 Balance arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.6 The flapper wiring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.7 Flapper . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.8 Helicopter plant with 2 degrees of freedom. See also Fig. 1.9(a). . . . . . . . . . . . 10
1.9 Helicopter control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.10 Helicopter flying . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.11 Real-time Simulink simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.1 The computer in the control loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.2 3 bit sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 A time series with unknown frequency components . . . . . . . . . . . . . . . . . . 18
2.4 The frequency component of a sampled signal . . . . . . . . . . . . . . . . . . . . . 18
2.5 Frequency aliases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.6 The Scarlet Letter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.7 Hénon’s attractor in Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.8 Hénon’s attractor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.9 Inverting z-transforms using dimpulse . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.10 Numerically inverting the Laplace transform using the Bromwich integral . . . . . 32
2.11 Numerically inverting the Laplace transform . . . . . . . . . . . . . . . . . . . . . . 33
2.12 Numerically inverting Laplace transforms . . . . . . . . . . . . . . . . . . . . . . . . 34
2.13 Ideal sampler and zeroth-order hold . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.14 Zeroth-order hold effects on the discrete Bode diagram . . . . . . . . . . . . . . . . 41
2.15 Bode plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.16 The discrete root locus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.17 Various discrete closed loop responses . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.18 A binary distillation column with multiple inputs and multiple outputs . . . . . . 45
2.19 A block diagram of a state-space dynamic system, (a) continuous system: ẋ =
Ax + Bu, and (b) discrete system: xk+1 = Φxk + ∆uk . (See also Fig. 2.20.) . . . . 46
2.20 A complete block diagram of a continuous state-space dynamic system with out-
puts and direct measurement feed-through, Eqn. 2.41. . . . . . . . . . . . . . . . . . 47
2.21 Unsteady and steady states for level systems . . . . . . . . . . . . . . . . . . . . . . 56
2.22 Submarine step response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
2.23 Issues in assessing system stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
2.24 Nyquist diagram of Eqn. 2.94 in (a) three dimensions and (b) as typically presented
in two dimensions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
2.25 Liapunov (1857–1918) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
2.26 Regions of stability for the poles of continuous (left) and discrete (right) systems . 83

ix
x LIST OF FIGURES

3.1 A stable and unstable pendulum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

3.2 Simple buffer tank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
3.3 The UK growth based on the GDP . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.4 A continually stirred tank reactor (CSTR) reactor . . . . . . . . . . . . . . . . . . . . 93
3.5 A forced circulation evaporator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3.6 Schematic of a distillation column . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
3.7 Wood-Berry step response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
3.8 Wood-Berry column in S IMULINK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.9 Distillation tower . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
3.10 Sparsity of the distillation column model . . . . . . . . . . . . . . . . . . . . . . . . 101
3.11 Open loop distillation column control . . . . . . . . . . . . . . . . . . . . . . . . . . 102
3.12 Distillation column control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
3.13 Distillation column control (in detail) . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
3.14 Distillation interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
3.15 Dynamic RGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
3.16 Dynamic RGA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
3.17 Density of Air . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3.18 Fitting a high-order polynomial to some physical data . . . . . . . . . . . . . . . . . 111
3.19 A bio-chemical reaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
3.20 Model of compressed water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
3.21 Experimental pressure-rate data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
3.22 Parameter confidence regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
3.23 Linear and nonlinear trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
3.24 Linearising a nonlinear tank model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

4.1 Comparing PI and integral-only control for the real-time control of a noisy flapper
plant with sampling time T = 0.08. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.2 PID simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.3 PID internals in Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.4 Block diagram of PID controllers as implemented in S IMULINK (left) and classical
text books (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.5 Realisable PID controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
4.6 PID controller in Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
4.7 PID controller with anti-derivative kick. . . . . . . . . . . . . . . . . . . . . . . . . . 139
4.8 Avoiding derivative kick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
4.9 Illustrating the improvement of anti-derivative kick schemes for PID controllers
when applied to the experimental electromagnetic balance. . . . . . . . . . . . . . . 140
4.10 Derivative control and noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.11 Anti-windup comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
4.12 Discrete PID controller in S IMULINK . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
4.13 Headbox control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
4.14 Headbox controlled response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
4.15 A PID controlled process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
4.16 Sample time and discrete PID control . . . . . . . . . . . . . . . . . . . . . . . . . . 149
4.17 The parameters T and L to be graphically estimated for the openloop tuning method
relations given in Table 4.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
4.18 Cohen-Coon model fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
4.19 Cohen-Coon tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
4.20 PID tuning using a GUI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
4.21 Solving for the ultimate frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
4.22 Ziegler-Nichols tuned responses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
4.23 Typical response of a stable system to a P-controller. . . . . . . . . . . . . . . . . . . 160
4.24 A Yuwana-Seborg closed loop step test . . . . . . . . . . . . . . . . . . . . . . . . . . 163
4.25 Closed loop responses using the YS scheme . . . . . . . . . . . . . . . . . . . . . . . 164
LIST OF FIGURES xi

4.26 A self-tuning PID controlled process . . . . . . . . . . . . . . . . . . . . . . . . . . . 167

4.27 A process under relay tuning with the PID regulator disabled. . . . . . . . . . . . . 167
4.28 An unknown plant under relay feedback exhibits an oscillation . . . . . . . . . . . 169
4.29 Nyquist & Bode diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
4.30 PID Relay tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172
4.31 Relay tuning with noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
4.32 Relay tuning of the blackbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
4.33 Relay tuning results of the blackbox . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
4.34 A relay with hysteresis width h and output amplitude d. . . . . . . . . . . . . . . . 177
4.35 Relay feedback with hysteresis width h. . . . . . . . . . . . . . . . . . . . . . . . . . 178
4.36 Relay feedback with hysteresis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
4.37 Relay feedback with an integrator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
4.38 2-point Relay identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
4.39 Monotonicity index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
4.40 The J curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
4.41 An inverse response process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
4.42 A NMP plant controlled with a PI controller . . . . . . . . . . . . . . . . . . . . . . 183
4.43 Approximating inverse-response systems with additional deadtime . . . . . . . . . 185
4.44 The Smith predictor structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186
4.45 The Smith predictor structure from Fig. 4.44 assuming no model/plant mis-match. 186
4.46 Smith predictor in Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
4.47 Dead time compensation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
4.48 Deadtime compensation applied to the blackbox . . . . . . . . . . . . . . . . . . . . 189
4.49 Closed loop with plant G(s) and controller C(s) subjected to disturbances and
measurement noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
4.50 Sensitivity transfer functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
4.51 Sensitivity robustness measure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191

5.1 A filter as a transfer function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193

5.2 A noisy measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
5.3 Noise added to a true, but unknown, signal . . . . . . . . . . . . . . . . . . . . . . . 195
5.4 Derivative action given noisy data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196
5.5 Smoothing industrial data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198
5.6 Low-pass filter specification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
5.7 Three single low-pass filters cascaded together to make a third-order filter. . . . . . 199
5.8 Amplitude response for ideal, low-pass, high pass and band-pass filters. . . . . . . 200
5.9 Analogue Butterworth filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
5.10 Analogue Chebyshev filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207
5.11 Butterworth and Chebyshev filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
5.12 Using a Butterworth filter to smooth noisy data . . . . . . . . . . . . . . . . . . . . . 211
5.13 The frequency response for Butterworth filters . . . . . . . . . . . . . . . . . . . . . 212
5.14 Various Butterworth filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213
5.15 Advantages of frequency pre-warping . . . . . . . . . . . . . . . . . . . . . . . . . . 214
5.16 Hardware difference equation in Direct Form I . . . . . . . . . . . . . . . . . . . . . 215
5.17 An IIR filter with a minimal number of delays, Direct Form II . . . . . . . . . . . . 216
5.18 Cascaded second-order sections to realise a high-order filter. See also Fig. 5.19. . . 217
5.19 A second-order section (SOS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 218
5.20 Comparing single precision second-order sections with filters in direct form II
transposed form. Note that the direct form II filter is actually unstable when run
in single precision. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 219
5.21 Approximating square waves with sine waves . . . . . . . . . . . . . . . . . . . . . 223
5.22 The Fourier approximation to a square wave . . . . . . . . . . . . . . . . . . . . . . 223
5.23 Two signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224
5.24 Critical radio frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228
xii LIST OF FIGURES

5.25 Power spectrum for a signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 229

5.26 Smoothing by Fourier transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230
5.27 Differentiating and smoothing noisy measurement . . . . . . . . . . . . . . . . . . . 232
5.28 Filtering industrial data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233

6.1 The prediction problem: Given a model and the input, u, can we predict the out-
put, y? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235
6.2 A good model, M, duplicates the behaviour of the true plant, S. . . . . . . . . . . . 237
6.3 An experimental setup for input/output identification. We log both the input and
the response data to a computer for further processing. . . . . . . . . . . . . . . . . 238
6.4 Typical open loop step tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 240
6.5 Areas method for model identification . . . . . . . . . . . . . . . . . . . . . . . . . . 241
6.6 Examples of the Areas method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242
6.7 Identification of the Blackbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242
6.8 Balance arm step test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243
6.9 Random signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244
6.10 A 5-element binary shift register to generate a pseudo-random binary sequence. . 245
6.11 Pseudo-random binary sequence generator in S IMULINK . . . . . . . . . . . . . . . 245
6.12 Black box experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 246
6.13 Black box response analysis using a series of sine waves . . . . . . . . . . . . . . . . 247
6.14 Black box response using an input chirp signal. . . . . . . . . . . . . . . . . . . . . . 248
6.15 Black box frequency response analysis using a chirp signal. . . . . . . . . . . . . . . 249
6.16 Flapper response to a chirp signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 249
6.17 Experimental setup to subject a random input into an unknown plant. The in-
put/output data was collected, processed through Listing 6.2 to give the frequency
response shown in Fig. 6.18. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250
6.18 The experimental frequency response compared to the true analytical Bode dia-
gram. See the routine in Listing 6.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250
6.19 Black box response given a pseudo-random input sequence. . . . . . . . . . . . . . 251
6.20 Black box frequency response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252
6.21 Empirical transfer function estimate . . . . . . . . . . . . . . . . . . . . . . . . . . . 253
6.22 Experimental data from a continuous plant . . . . . . . . . . . . . . . . . . . . . . . 254
6.23 A continuous-time model fitted to input/output data . . . . . . . . . . . . . . . . . 255
6.24 Continuous model identification strategy . . . . . . . . . . . . . . . . . . . . . . . . 257
6.25 Continuous model identification simulation . . . . . . . . . . . . . . . . . . . . . . . 258
6.26 Continuous model identification of the blackbox . . . . . . . . . . . . . . . . . . . . 259
6.27 Identification using Laguerre functions . . . . . . . . . . . . . . . . . . . . . . . . . 261
6.28 A signal flow diagram of an auto-regressive model with exogenous input or ARX
model. Compare this structure with the similar output-error model in Fig. 6.31. . . 262
6.29 A signal flow diagram of a ARMAX model. Note that the only difference between
this, and the ARX model in Fig. 6.28, is the inclusion of the C polynomial filtering
the noise term. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263
6.30 Discrete time models implemented in S IMULINK . . . . . . . . . . . . . . . . . . . . 264
6.31 A signal flow diagram of an output-error model. Compare this structure with the
similar ARX model in Fig. 6.28. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 265
6.32 A general input/output model structure . . . . . . . . . . . . . . . . . . . . . . . . . 265
6.33 Kaos, Control and regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 268
6.34 ARX estimation exhibiting bias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 271
6.35 Offline system identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272
6.36 Fitting a polynomial to some (x, y) data . . . . . . . . . . . . . . . . . . . . . . . . . 278
6.37 Using the AIC to establish model structure . . . . . . . . . . . . . . . . . . . . . . . 278
6.38 Using the AIC to establish model structure using separate validation data . . . . . 279
6.39 Identification of deadtime from the step response . . . . . . . . . . . . . . . . . . . 281
6.40 Deadtime estimation at fast sampling . . . . . . . . . . . . . . . . . . . . . . . . . . 281
LIST OF FIGURES xiii

6.41 Deadtime estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 281

6.42 Blackbox step model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282
6.43 Hammerstein and Wiener model structures . . . . . . . . . . . . . . . . . . . . . . . 284
6.44 Ideal RLS parameter estimation. (See also Fig. 6.46(a).) . . . . . . . . . . . . . . . . 289
6.45 Recursive least squares estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 290
6.46 RLS under S IMULINK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 291
6.47 RLS under S IMULINK (version 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292
6.48 Estimation of a two parameter plant . . . . . . . . . . . . . . . . . . . . . . . . . . . 294
6.49 Confidence limits for estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294
6.50 RLS and an abrupt plant change . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 295
6.51 The ‘memory’ of the forgetting factor . . . . . . . . . . . . . . . . . . . . . . . . . . . 297
6.52 Identification using various forgetting factors . . . . . . . . . . . . . . . . . . . . . . 298
6.53 Covariance wind-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 300
6.54 Variable forgetting factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 301
6.55 The MIT rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304
6.56 Optimising the adaptation gain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305
6.57 Addition of coloured noise to a dynamic process. See also Fig. 6.32. . . . . . . . . . 305
6.58 RLS with coloured noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 307
6.59 Recursive extended least-squares estimation . . . . . . . . . . . . . . . . . . . . . . 308
6.60 Recursive extended least-squares estimation . . . . . . . . . . . . . . . . . . . . . . 309
6.61 A simplified recursive least squares algorithm . . . . . . . . . . . . . . . . . . . . . 311
6.62 Furnace input/output data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312
6.63 Closed loop estimation using Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . 313
6.64 RLS in Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 319

7.1 The structure of an indirect adaptive controller . . . . . . . . . . . . . . . . . . . . . 323

7.2 Varying process gain of a spherical tank . . . . . . . . . . . . . . . . . . . . . . . . . 324
7.3 Simple minimum variance control . . . . . . . . . . . . . . . . . . . . . . . . . . . . 328
7.4 Simple minimum variance control (zoomed) . . . . . . . . . . . . . . . . . . . . . . 329
7.5 Adaptive pole-placement control structure . . . . . . . . . . . . . . . . . . . . . . . 329
7.6 Adaptive pole-placement control structure with RLS identification . . . . . . . . . 335
7.7 Control of multiple plants with an adapting controller. We desire the same closed
loop response irrespective of the choice of plant. . . . . . . . . . . . . . . . . . . . . 336
7.8 Three open loop plants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336
7.9 Desired closed loop response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336
7.10 Adaptive pole-placement with identification . . . . . . . . . . . . . . . . . . . . . . 338
7.11 Comparing the adaptive pole-placement with the reference trajectory . . . . . . . . 338
7.12 Adaptive pole-placement with identification in S IMULINK . . . . . . . . . . . . . . 339
7.13 Internals of the Adaptive pole-placement design routine . . . . . . . . . . . . . . . 339
7.14 Bursting in adaptive control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 341
7.15 A plant with poorly damped zeros . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342
7.16 Adaptive pole-placement with an unstable B . . . . . . . . . . . . . . . . . . . . . . 344
7.17 Migration of poles and zeros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 345
7.18 Areas of well-damped poles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 346
7.19 Adaptive pole-placement of the black-box . . . . . . . . . . . . . . . . . . . . . . . . 347
7.20 Adaptive pole-placement of the black-box . . . . . . . . . . . . . . . . . . . . . . . . 348
7.21 A non-minimum phase plant with an unstable zero which causes an inverse step
response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 350
7.22 Moving average control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 351

8.1 Reconstructing states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 366

8.2 Pole-placement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 370
8.3 Pole-placement degradation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 370
8.4 Deadbeat control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 372
xiv LIST OF FIGURES

8.5 Simultaneous control and state estimation . . . . . . . . . . . . . . . . . . . . . . . . 375

8.6 Control and estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 377
8.7 GMC tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 379
8.8 GMC tuning characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 379
8.9 Linear GMC response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 381
8.10 Linear GMC comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 381
8.11 GMC CSTR control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 385
8.12 A CSTR phase plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 386
8.13 The configuration of an input/output feedback linearisation control law . . . . . . 388
8.14 Exact feedback linearisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392

9.1 IAE areas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 397

9.2 ITAE breakdown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 398
9.3 Optimal responses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 398
9.4 Optimal PID tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 401
9.5 Optimum PI tuning of the blackbox plant . . . . . . . . . . . . . . . . . . . . . . . . 402
9.6 S IMULINK model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403
9.7 Production of a valuable chemical in a batch reactor. . . . . . . . . . . . . . . . . . . 404
9.8 Temperature profile optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405
9.9 Temperature profile optimisation using 3 temperatures . . . . . . . . . . . . . . . . 406
9.10 Optimum temperature profile comparison for different number of temperatures . 406
9.11 Optimal control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 412
9.12 Optimal control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414
9.13 Optimal control with targets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 416
9.14 Steady-state and time-varying LQR control . . . . . . . . . . . . . . . . . . . . . . . 422
9.15 Steady-state continuous LQR controller . . . . . . . . . . . . . . . . . . . . . . . . . 425
9.16 Comparing discrete and continuous LQR controllers . . . . . . . . . . . . . . . . . . 428
9.17 LQR control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 431
9.18 Pole-placement and LQR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 432
9.19 Pole-placement and LQR showing the input . . . . . . . . . . . . . . . . . . . . . . 433
9.20 Trial pole locations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433
9.21 Trial pole-placement performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434
9.22 State feedback control system with an integral output state . . . . . . . . . . . . . . 436
9.23 State feedback with integral states . . . . . . . . . . . . . . . . . . . . . . . . . . . . 437
9.24 Black box servo control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 438
9.25 A state-based estimation and control scheme . . . . . . . . . . . . . . . . . . . . . . 439
9.26 PDF of a random variable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 440
9.27 Correlated noisy x, y data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443
9.28 2D ellipse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444
9.29 A block diagram of a stochastic state-space system . . . . . . . . . . . . . . . . . . . 445
9.30 Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 451
9.31 A block diagram of a steady-state prediction-type Kalman filter applied to a linear
discrete plant. Compare this with the alternative form in Fig. 9.32. . . . . . . . . . . 453
9.32 A block diagram of a steady-state current estimator-type Kalman filter applied to
a linear discrete plant. Compare with the alternative prediction form in Fig. 9.31. . 454
9.33 Kalman filter demonstration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 457
9.34 The performance of a Kalman filter for different q/r ratios . . . . . . . . . . . . . . 460
9.35 A random walk process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464
9.36 LQG in S IMULINK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 466

10.1 Horizons used in model predictive control . . . . . . . . . . . . . . . . . . . . . . . 470

10.2 Predictions of the Reserve Bank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 470
10.3 Inverse plant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473
10.4 Predictive control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474
LIST OF FIGURES xv

10.5 Acausal response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 475

10.6 Varying the horizons of predictive control . . . . . . . . . . . . . . . . . . . . . . . . 476
10.7 MPC on the blackbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 476
10.8 Step response coefficients, gi , for a stable system. . . . . . . . . . . . . . . . . . . . . 477
10.9 DMC control details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 479
10.10DMC control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 480
10.11Adaptive DMC of the blackbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 481
10.12An MPC graphical user interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 482
10.13Multivariable MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484
10.14S IMULINK and MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484
10.15MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485
10.16A 3 degree of freedom helicopter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 486
10.17MPC control of a helicopter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 487
10.18MPC control of a helicopter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 489
10.19LP constraint matrix dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 492
10.20LP optimal control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 496
10.21LP optimal control with active constraints . . . . . . . . . . . . . . . . . . . . . . . . 496
10.22Non-square LP optimal control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 497
10.23LP optimal control showing acausal behaviour . . . . . . . . . . . . . . . . . . . . . 497

11.1 Possible neuron activation functions. . . . . . . . . . . . . . . . . . . . . . . . . . . 507

11.2 A single neural processing unit with multiple inputs . . . . . . . . . . . . . . . . . . 508
11.3 Single layer feedforward neural network . . . . . . . . . . . . . . . . . . . . . . . . 509
11.4 A 3 layer fully interconnected feedforward neural network . . . . . . . . . . . . . . 509
11.5 Single layer network with feedback . . . . . . . . . . . . . . . . . . . . . . . . . . . . 510
11.6 An unknown input/output function . . . . . . . . . . . . . . . . . . . . . . . . . . . 511
11.7 Fitting a Neural-Network to an unknown input/output function . . . . . . . . . . 513
11.8 Tide predictions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514
xvi LIST OF FIGURES
L IST OF TABLES

1.1 Computer aids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.1 Final and initial value theorems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.2 Inverting a z transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3 Laplace transform pairs used for testing . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.1 Standard nomenclature used in modelling dynamic systems . . . . . . . . . . . . . 90

3.2 Parameters of the CSTR model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.3 The important variables in the forced circulation evaporator . . . . . . . . . . . . . 94
3.4 Compressed water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
3.5 The parameter values for the CSTR model . . . . . . . . . . . . . . . . . . . . . . . . 122
3.6 The initial state and manipulated variables for the CSTR simulation . . . . . . . . . 123

4.1 Alternative PID tuning parameter conventions . . . . . . . . . . . . . . . . . . . . . 132

4.2 Ziegler-Nichols open-loop PID tuning rules . . . . . . . . . . . . . . . . . . . . . . . 152
4.3 PID controller settings based on IMC for a small selection of common plants where
the control engineer gets to chose a desired closed loop time constant, τc . . . . . . . 154
4.4 Various alternative ‘Ziegler-Nichols’ type PID tuning rules as a function of the
ultimate gain, Ku , and ultimate period, Pu . . . . . . . . . . . . . . . . . . . . . . . . 155
4.5 Closed-loop single-test PID design rules . . . . . . . . . . . . . . . . . . . . . . . . . 161
4.6 Relay based PID tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172

5.1 Filter transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201

6.1 Experimentally determined frequency response of the blackbox . . . . . . . . . . . 247

6.2 Identification in state-space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 277

8.1 The relationship between regulation and estimation . . . . . . . . . . . . . . . . . . 374

8.2 Litchfield nonlinear CSTR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382

9.1 Common integral performance indices . . . . . . . . . . . . . . . . . . . . . . . . . . 399

11.1 Comparing expert systems and neural networks . . . . . . . . . . . . . . . . . . . . 506

xvii
xviii LIST OF TABLES
L ISTINGS

2.1 Symbolic Laplace to z-transform conversion . . . . . . . . . . . . . . . . . . . . . . 37

2.2 Symbolic Laplace to z-transform conversion with ZOH . . . . . . . . . . . . . . . . 37
2.3 Extracting the gain, time constants and numerator time constants from an arbitrary
transfer function format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
2.4 Submarine simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
2.5 Example of the Routh array using the symbolic toolbox . . . . . . . . . . . . . . . . 71
2.6 Solve the continuous matrix Lyapunov equation using Kronecker products . . . . 77
2.7 Solve the matrix Lyapunov equation using the lyap routine . . . . . . . . . . . . . 78
2.8 Using Lyapunov to establish stability of a linear system . . . . . . . . . . . . . . . . 79
2.9 Solve the discrete matrix Lyapunov equation using Kronecker products . . . . . . 79
3.1 Computing the dynamic relative gain array analytically . . . . . . . . . . . . . . . . 105
3.2 Computing the dynamic relative gain array numerically as a function of ω. See
also Listing 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
3.3 Curve fitting using polynomial least-squares . . . . . . . . . . . . . . . . . . . . . . 109
3.4 Polynomial least-squares using singular value decomposition. This routine fol-
lows from, and provides an alternative to Listing 3.3. . . . . . . . . . . . . . . . . . 111
3.5 Curve fitting using a generic nonlinear optimiser . . . . . . . . . . . . . . . . . . . . 113
3.6 Curve fitting using the O PTI optimisation toolbox. (Compare with Listing 3.5.) . . 113
3.7 Fitting water density as a function of temperature and pressure . . . . . . . . . . . 115
3.8 Parameter confidence limits for a nonlinear reaction rate model . . . . . . . . . . . 117
3.9 Comparing the dynamic response of a pendulum to the linear approximation . . . 120
3.10 Using linmod to linearise an arbitrary S IMULINK module. . . . . . . . . . . . . . . 127
4.1 Constructing a transfer function of a PID controller . . . . . . . . . . . . . . . . . . 133
4.2 Constructing a discrete (filtered) PID controller . . . . . . . . . . . . . . . . . . . . . 145
4.3 A simple PID controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
¿ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
4.4 Ziegler-Nichols PID tuning rules for an arbitrary transfer function . . . . . . . . . . 159
4.5 Identifies the characteristic points for the Yuwana-Seborg PID tuner from a trial
closed loop response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
4.6 Compute the closed loop model from peak and trough data . . . . . . . . . . . . . 165
4.7 Compute the ultimate gain and frequency from the closed loop model parameters. 165
4.8 Compute the open loop model, Gm , Eqn. 4.31. . . . . . . . . . . . . . . . . . . . . . 165
4.9 Compute appropriate PI or PID tuning constants based on a plant model, Gm ,
using the IMC schemes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
4.10 Calculates the period and amplitude of a sinusoidal time series using least-squares. 174
5.1 Designing Butterworth Filters using Eqn. 5.4. . . . . . . . . . . . . . . . . . . . . . . 203
5.2 Designing a low-pass Butterworth filter with a cut-off frequency of fc = 800 Hz. . 203
5.3 Designing a high-pass Butterworth filter with a cut-off frequency of fc = 800 Hz. . 203
5.4 Designing Chebyshev Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
5.5 Computing a Chebyshev Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207
5.6 Converting a 7th-order Butterworth filter to 4 second-order sections . . . . . . . . . 218
5.7 Comparing DFII and SOS digital filters in single precision. . . . . . . . . . . . . . . 219
5.8 Routine to compute the power spectral density plot of a time series . . . . . . . . . 226

xix
xx LISTINGS

5.9 Smoothing and differentiating a noisy signal . . . . . . . . . . . . . . . . . . . . . . 232

6.1 Identification of a first-order plant with deadtime from an openloop step response
using the Areas method from Algorithm 6.1. . . . . . . . . . . . . . . . . . . . . . . 241
6.2 Frequency response identification of an unknown plant directly from input/out-
put data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250
6.3 Non-parametric frequency response identification using etfe. . . . . . . . . . . . . 252
6.4 Function to generate output predictions given a trial model and input data. . . . . 255
6.5 Optimising the model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255
6.6 Validating the fitted model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255
6.7 Continuous model identification of a non-minimum phase system . . . . . . . . . . 257
6.8 Generate some input/output data for model identification . . . . . . . . . . . . . . 269
6.9 Estimate an ARX model from an input/output data series using least-squares . . . 270
6.10 An alternative way to construct the data matrix for ARX estimation using Toeplitz
matrices. See also Listing 6.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 270
6.11 Offline system identification using arx from the System Identification Toolbox . . 272
6.12 Offline system identification with no model/plant mismatch . . . . . . . . . . . . . 272
6.13 Demonstrate the fitting of an AR model. . . . . . . . . . . . . . . . . . . . . . . . . . 273
6.14 Create an input/output sequence from an output-error plant. . . . . . . . . . . . . 274
6.15 Parameter identification of an output error process using oe and arx. . . . . . . . 274
6.16 A basic recursive least-squares (RLS) update (without forgetting factor) . . . . . . . 287
6.17 Tests the RLS identification scheme using Listing 6.16. . . . . . . . . . . . . . . . . . 290
6.18 A recursive least-squares (RLS) update with a forgetting factor. (See also List-
ing 6.16.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 296
6.19 Adaption of the plant gain using steepest descent . . . . . . . . . . . . . . . . . . . 303
6.20 Create an ARMAX process and generate some input/output data suitable for sub-
sequent identification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 307
6.21 Identify an ARMAX process from the data generated in Listing 6.20. . . . . . . . . 308
6.22 Recursively identify an ARMAX process. . . . . . . . . . . . . . . . . . . . . . . . . 308
6.23 Kaczmarz’s algorithm for identification . . . . . . . . . . . . . . . . . . . . . . . . . 310
7.1 Simple minimum variance control where the plant has no time delay . . . . . . . . 327
7.2 A Diophantine routine to solve F A + BG = T for the polynomials F and G. . . . . 332
7.3 Alternative Diophantine routine to solve F A + BG = T for the polynomials F and
G. Compare with Listing 7.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332
7.4 Constructing polynomials for the Diophantine equation example . . . . . . . . . . 333
7.5 Solving the Diophantine equation using polynomials generated from Listing 7.4. . 334
7.6 Adaptive pole-placement control with 3 different plants . . . . . . . . . . . . . . . . 337
7.7 The pole-placement control law when H = 1/B . . . . . . . . . . . . . . . . . . . . 341
7.8 Factorising an arbitrary polynomial B(q) into stable, B + (q), and unstable and
poorly damped, B − (q), factors such that B = B + B − and B + is defined as monic. . 345
7.9 Minimum variance control design . . . . . . . . . . . . . . . . . . . . . . . . . . . . 349
8.1 A simple state reconstructor following Algorithm 8.1. . . . . . . . . . . . . . . . . . 365
8.2 Pole-placement control of a well-behaved system . . . . . . . . . . . . . . . . . . . . 369
8.3 A deadbeat controller simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 371
8.4 Pole placement for controllers and estimators . . . . . . . . . . . . . . . . . . . . . . 375
8.5 GMC on a Linear Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 380
8.6 GMC for a batch reactor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383
8.7 The dynamic equations of a batch reactor . . . . . . . . . . . . . . . . . . . . . . . . 384
8.8 Find the Lie derivative for a symbolic system . . . . . . . . . . . . . . . . . . . . . . 390
8.9 Establish relative degree, r (ignore degree 0 possibility) . . . . . . . . . . . . . . . . 390
8.10 Design Butterworth filter of order r. . . . . . . . . . . . . . . . . . . . . . . . . . . . 390
8.11 Symbolically create the closed loop expression . . . . . . . . . . . . . . . . . . . . . 391
9.1 Returns the IAE performance for a given tuning. . . . . . . . . . . . . . . . . . . . . 400
9.2 Optimal tuning of a PID controller for a non-minimum phase plant. This script file
uses the objective function given in Listing 9.1. . . . . . . . . . . . . . . . . . . . . . 400
LISTINGS xxi

9.3 Returns the ITSE using a S IMULINK model. . . . . . . . . . . . . . . . . . . . . . . . 403

9.4 Analytically computing the co-state dynamics and optimum input trajectory as a
function of states and co-states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413
9.5 Solving the reaction profile boundary value problem using the boundary value
problem solver, bvp4c.m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413
9.6 Computes the full time-evolving LQR solution . . . . . . . . . . . . . . . . . . . . . 421
9.7 The continuous time differential Riccati equation. This routine is called from List-
ing 9.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 421
9.8 Solves the continuous time differential Riccati equation using a numerical ODE
integrator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 422
9.9 Calculate the continuous optimal steady-state controller gain. . . . . . . . . . . . . 424
9.10 Closed loop simulation using an optimal steady-state controller gain. . . . . . . . . 424
9.11 Solving the algebraic Riccati equation for P∞ using Kronecker products and vec-
torisation given matrices A, B, Q and R. . . . . . . . . . . . . . . . . . . . . . . . . 425
9.12 Calculate the discrete optimal steady-state gain by ‘iterating until exhaustion’.
Note it is preferable for numerical reasons to use lqr for this computation. . . . . 427
9.13 Comparing the continuous and discrete LQR controllers. . . . . . . . . . . . . . . . 427
9.14 An LQR controller for the blackbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . 431
9.15 Comparing an LQR controller from Listing 9.14 with a pole-placement controller . 432
9.16 Computing the closed loop poles from the optimal LQR controller from Listing 9.14.434
9.17 Comparing the actual normally distributed random numbers with the theoretical
probability density function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 440
9.18 Probability and inverse probability distributions for the F -distribution. . . . . . . . 442
9.19 Generate some correlated random data. . . . . . . . . . . . . . . . . . . . . . . . . . 442
9.20 Plot a 3D histogram of the random data from Listing 9.19. . . . . . . . . . . . . . . 443
9.21 Compute the uncertainty regions from the random data from Listing 9.20. . . . . . 444
9.22 Validating the uncertainty regions computed theoretically from Listing 9.21. . . . . 444
9.23 Solving the discrete time Riccati equation using exhaustive iteration around Eqn. 9.98
or alternatively using the dare routine. . . . . . . . . . . . . . . . . . . . . . . . . . 451
9.24 Alternative ways to compute the Kalman gain . . . . . . . . . . . . . . . . . . . . . 453
9.25 State estimation of a randomly generated discrete model using a Kalman filter. . . 456
9.26 Computing the Kalman gain using dlqe. . . . . . . . . . . . . . . . . . . . . . . . . 457
9.27 Demonstrating the optimality of the Kalman filter. . . . . . . . . . . . . . . . . . . . 458
9.28 Potter’s algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 462
10.1 Predictive control with input saturation constraints using a generic nonlinear op-
timiser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473
10.2 Objective function to be minimised for the predictive control algorithm with input
saturation constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474
10.3 Dynamic Matric Control (DMC) control . . . . . . . . . . . . . . . . . . . . . . . . . 479
10.4 Setting up an MPC controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483
10.5 Nonlinear dynamics of a helicopter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 486
10.6 Optimal control using linear programming . . . . . . . . . . . . . . . . . . . . . . . 494
11.1 Generate some arbitrary data to be used for subsequent fitting . . . . . . . . . . . . 512
B.1 Polynomial addition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 519
B.2 Multiple convolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 519
B.3 Strip leading zeros from a polynomial. . . . . . . . . . . . . . . . . . . . . . . . . . . 520
xxii LISTINGS
C HAPTER 1

I NTRODUCTION

Mathematicians may flatter themselves that they posses new ideas which mere human language is as yet
unable to express. Let them make the effort to express those ideas in appropriate words without the aid of
symbols, and if they succeed, they will not only lay us laymen under a lasting obligation but, we venture to
say, they will find themselves very much enlightened during the process, and will even be doubtful
whether the ideas expressed as symbols had ever quite found their way out of the equations into their
minds.
James Clerk Maxwell, 1890

Control, in an engineering sense, is where actions are taken to ensure that a particular physi-
cal process responds in some desired manner. Automatic control is where we have relieved the
human operator from the tedium of consistently monitoring the process and supplying the nec-
essary corrections. Control as a technical discipline is therefore important not only in the fields
of engineering, but also in economics, sociology and indeed in most aspects of our life. When
studying control, we naturally assume that we do conceivably have some chance of influencing
things. For example, it is worthwhile to study the operation of a coal fired power plant in order
to minimise possibly polluting emissions, but it is not worth our time to save the world from the
next ice age, or as the results of a special study group who investigated methods designed to pro-
tect the world from a stray comet (such as the one postulated to have wiped out the dinosaurs 80
million years ago) concluded, there was nothing feasible we could do, such as change the earth’s
orbit, or blast the asteroid, to avoid the collision. In these latter examples, the problem exists, but
our influence is negligible.

The teaching of control has changed in emphasis over the last decade from that of linear single-
input/single-output systems elegantly described in the Laplace domain, to general nonlinear
multiple-input/multiple-output systems best analysed in the state space domain. This change
has been motivated by the increasing demands by industry and public to produce more, faster or
cleaner and is now much more attractive due to the impressive improvements in computer aided
tools, such as M ATLAB used in these notes. This new emphasis is called advanced (or modern)
control as opposed to the traditional or classical control. This set of notes is intended for students
who have previously attended a first course in automatic control covering the usual continuous
control concepts like Laplace transforms, Bode diagrams, stability of linear differential equations,
PID controllers, and perhaps some exposure to discrete time control topics like z transforms and
state space.

This book attempts to describe what advanced control is, and how it is applied in engineering
applications with emphasis on the construction of controllers using computer aided design tools
such as the numerical programming environment M ATLAB from the MathWorks, [136]. With

1
2 CHAPTER 1. INTRODUCTION

this tool, we can concentrate on the intentions behind the design procedures, rather than the
mechanics to follow them.

1.1 S COPE

Part one contains some revision material in z-transforms, modelling and PID controller tuning.
The discrete domain, z-transforms, and stability concepts with a brief discussion of appropriate
numerical methods are introduced in chapter 2. A brief potpourri of modelling is summarised
in chapter 3. Chapter 4 is devoted to the most common industrial controller, the three term PID
controller with emphasis on tuning, implementation and limitations. Some basic concepts from
signal processing such as filtering and smoothing are introduced in chapter 5. Identification and
the closely related adaptive control are together in chapters 6 and 7. State space analysis and
optimal control design are given in the larger chapters 8 and 9.

N OTATION CONVENTIONS

Throughout these notes I have used some typographical conventions. In mathematical expres-
sions, scalar variables are written in italic such as a, b, c, or if Greek, γ, ϕ while vectors x, y are
upright bold lower case and matrices, A, ∆ are bold upper case. More notation is introduced as
required.

Computer commands and output and listings are given in a fixed-width font as A=chol(B*B’).
In some cases where you are to type an interactive command, the M ATLAB prompt >> is given,
and the computed solution returned. If no ambiguity exists such as in the case for functions, the
prompt is omitted.

1.2 M ATLAB FOR COMPUTER AIDED CONTROL DESIGN

Modern control design has heavy computing requirements. In particular one needs to:

1. manipulate symbolic algebraic expressions, and

2. perform intensive numerical calculations and simulations for proto-typing and testing quickly
and reliably, and finally

3. to implement the controller at high speed in special hardware such as an embedded con-
troller or a digital signal processing (DSP) chip perhaps using assembler.

To use this new theory, it is essential to use computer aided design (CAD) tools efficiently as
real-world problems can rarely be solved manually. But as [172] point out, ‘the use of computers
in the design of control systems has a long and fairly distinguished history’. This book uses
M ATLAB for the design, simulation and prototyping of controllers.

M ATLAB, (which is short for MATrix LABoratory), is a programming environment that grew
out of an effort to create an easy user-interface to the very popular and well regarded public
domain F ORTRAN linear algebra collection of programmes, L INPACK and E ISPACK. With this
direct interpretive interface, one can write quite sophisticated algorithms in a very high level
language, that are consequently easy to read and maintain. M ATLAB is supported with a variety
1.2. MATLAB FOR COMPUTER AIDED CONTROL DESIGN 3

of toolboxes comprising of collections of source code subroutines organised in areas of specific

interest. The toolboxes we are most interested in, and used in this book are:

Control toolbox containing functions for controller design, frequency domain analysis, conver-
sions between various models forms, pole placement, optimal control etc. (Used through-
out)
Symbolic toolbox which contains a gateway to the symbolic capabilities of M APLE.
Signal processing toolbox containing filters, wave form generation and spectral analysis. (Used
principally in chapter 5.)
System identification toolbox for identifying the parameters of various dynamic model types.
(Used in chapter 6.) You may also find the following free statistics toolbox useful available
at: www.maths.lth.se/matstat/stixbox/.
Optimisation which contains routines to maximise or minimise algebraic or dynamic functions.
One suitable free toolbox is available from www.i2c2.aut.ac.nz/Wiki/OPTI/. This
toolbox is used primarily in chapters 9 and 10.
Real-time toolbox can be used to interface M ATLAB to various analogue to digital converters. A
free version suitable with Measurement Computing’s MCC DAQs is available from www.i2c2.aut.ac.nz/R

Additional documentation to that supplied with M ATLAB is the concise and free summary notes
[186] or the more recent [69]. Recently there has been exponential growth of other texts that heav-
ily use M ATLAB (such as this one), and a current list is available from the Mathwork’s anonymous
ftp server at www.mathworks.com. This server also contains many user contributed codes, as
well as updates, bug fixes etc.

If M ATLAB, or even programming a high level language is new to you, then [205] is a cheap
recommended compendium, similar in form to this, covering topics in numerical analysis, again
with many M ATLAB examples.

1.2.1 A LTERNATIVE COMPUTER DESIGN AIDS

Table 1.1 lists a number of alternatives computer-aided design and modelling environments sim-
ilar and complimentary to M ATLAB.

Product WWW site comment

S CI L AB www.scilab.org Free Matlab/Simulink clone
O CTAVE www.octave.org Free Matlab clone, inactive
RL AB rlabplus.sourceforge.net Matlab clone, Linux
V ISUAL M ODEL Q www.qxdesign.com shareware Simulink clone
M ATH V IEWS www.mathwizards.com shareware
M U PAD www.mupad.de Interfaces with S CI L AB
M APLE www.maplesoft.com commercial CAS
M ATHEMATICA www.mathematica.com commercial CAS

Table 1.1: Shareware or freeware Matlab lookalikes and computer algebra systems

Unlike M ATLAB, symbolic manipulators are computer programs that by manipulating symbols
can perform algebra. Such programs are alternatively known as computer algebra systems or
4 CHAPTER 1. INTRODUCTION

CAS. The most well known examples are M ATHEMATICA, M APLE, M U PAD, and M ACSYMA,
(see Table 1.1). These programs can find analytical solutions to many mathematical problems
involving integrals, limits, special functions and so forth. They are particularly useful in the
controller design stage.

The Numerics in Control1 group in Europe have collect together a freeware F ORTRAN subroutine
library S LICOT for routines relevant in systems and control.

Problem 1.1 1. Familiarise yourself with the fundamentals of M ATLAB. Run the M ATLAB
demo by typing demo once inside M ATLAB.

2. Try the M ATLAB tutorial (part 1).

3. Read through the M ATLAB primer, [186] or [69], and you should get acquainted with the
M ATLAB users manual.

1.3 E CONOMICS OF CONTROL

Most people would agree that Engineers apply technology, but what do these two words really
mean? Technology is derived from two Greek words, techne which means skill or art, and logia
which means science or study. The interesting point here is that the art component is included.
The English language unfortunately confuses the word engine with engineering so that many peo-
ple have the mistaken view that engineers drive engines (mostly). Actually engineer is derived
from the Latin ingeniatorium which means one who is ingenious at devising. A far cry from the
relatively simple act of piloting jumbos. An interesting American perspective of the professional
engineer and modern technology is given as light reading in [2] and Florman’s The Existential
Pleasures of Engineering, [67].

Chemical engineering is, succinctly put, chemistry constrained by cost. The chemist wants the
reaction to proceed, the chemical engineer takes that for granted, but is interested in increasing
the rate, or pushing the equilibrium, or in most cases both at the same time. As the incentive to
produce better products increases, accompanied by an awareness of potent global competition
driving one to reduce costs, process control becomes an important aspect of engineering.

Obviously modern computerised process control systems are expensive. They are especially ex-
pensive compared with other computers such as office or financial computers because the market
is smaller, the environment harsher, the graphical requirements more critical, the duties more var-
ied, and the potential payoffs larger. Process control has at least two main duties to perform; first
to ensure that the plant is operated safely (that is protect the plant, environment and people), and
second that the product quality is consistent with some customer or regulatory body demanded
specifications. There is always a trade off between how much control you apply and the bene-
fits that result. An automobile manufacturer could produce an almost totally indestructible car
(ie, tank), but the expense of raw materials required, and the high running costs would certainly
deem the project an economic failure.

On the other hand, in 1965 Ralph Nader complained in the aptly named Unsafe at Any Speed
about the poor quality of American automotive engineering, the lack of controls and the unsafe
result. This influential book challenged the balance from between commercial profits and more
quality control. A product with less variation in quality may be worth more than a product that
has a higher average quality, but more variation. A potentially dangerous production facility that
regularly destroys equipment, people or surroundings is not usually tolerated by the licensing
authorities.
1 The home page is located at http://www.win.tue.nl/niconet/niconet.html
1.3. ECONOMICS OF CONTROL 5

Safety concerns motivate

better control.

Fig. 1.1, adapted from [6], gives an industrial perspective of the status of process control in 1994.
The techniques are divided into those considered “classical” or traditional, which demand only
modest digital online computing power, (if any), little in the way of explicit process models or
understanding, and those termed loosely “advanced”.

Advanced control

Optimisation Direct search Rule definition Online simulation

Constraint Steady state dynamic methods single variable multivariable

Regulatory Signal condition PID algorithm Deadtime feedforward control

compensation

Basic analogue control PLCs DCS Process computer

Field Valves Transmitters smart transmitters Onstream analysers

Traditional control

Figure 1.1: A comparison of traditional vs. advanced process control techniques. Adapted from
[6].

One of the major concerns for the process control engineer is to reduce the output variance. If
the variation about the setpoint is small, then the setpoint can be shifted closer to the operating
constraint, without increasing the frequency of alarms. Fig. 1.2 demonstrates the ideal case that
while popular in the advertising literature, is harder to achieve unambiguously in practice.

Many text books in the control field are very vague about the actual configuration of real pro-
cess control systems used today. Other books that take a more trade oriented approach, are
vague about the academic side of the control performance. There are a number of reasons for
this. First many texts try hard to describe only the theoretical aspects of digital control, and any-
thing remotely applied is not considered worthy of their attention. Secondly, the control systems
are rapidly changing as the cost of micro-processors drop in price, and different programming
methods come into flavour. Thirdly many industries are deliberately vague about publishing
the details of their control system since they perceive that this information could help their com-
petitors. However some information of this type is given in [64, pp131-149] and [17]. One good
6 CHAPTER 1. INTRODUCTION

upper quality constraint

violations

spt #3
loss ($) spt #2

setpoint #1

process output advanced control

regulatory

manual control

time

Figure 1.2: Economic improvements owing to better control. If the control scheme can reduce
the variance, the setpoint can be shifted closer to the operating or quality constraint, thereby
decreasing operating costs.

balance for the practitioner is [145].

1.4 L ABORATORY EQUIPMENT FOR CONTROL TESTS

Obviously if we are to study automatic control with the aim to control eventually chemical plants,
manufacturing processes, robots, undertake filtering to do active noise cancellation and so forth,
we should practice, preferably on simpler, more well understood, and potentially less hazardous
equipment.

In the Automatic Control Laboratory in the Department of Electrical Engineering at the Karlstad
University, Sweden we have a number of simple bench-scale plants to test identification and
control algorithms on.

1.4.1 P LANTS WITH ONE INPUT AND ONE OUTPUT

T HE BLACKBOX

Fig. 1.3 and Fig. 1.4(a) shows what we perhaps unimaginatively refer to as a “black-box’. It is a
box, and it is coloured black. Subjecting the box to an input voltage from 0 to 5 volts delivers an
output voltage also spanning from around 0 to 5 volts, but lagging behind the input voltage since
the internals of the blackbox are simply either 7 or 9 (depending on the switch position) low-pass
passive filters cascaded together.

The blackbox is a relatively well behaved underdamped stable system with dominant time con-
stants of around 5 to 10 seconds. Fig. 1.4(b) shows the response of the blackbox to two input
steps. The chief disadvantage of using this device for control studies is that the output response
1.4. LABORATORY EQUIPMENT FOR CONTROL TESTS 7

Input Output
From computer
To computer
D/A ✲ B LACK -B OX Sluggish ✲ A/D
(connector #10) (connector #1)

✲
Fast
GND ✲ input indicator blue ✲ Earth
blue
(connector #11)

Figure 1.3: Blackbox configuration. The manual switch marked will toggle between either 7 or 9
low-pass filters.

Blackbox step response

0.45

0.4

0.35

0.3
ipnut/output

0.25

0.2

0.15 input
output
0.1

0.05

−0.05
15 20 25 30 35 40 45 50
time (s)

(a) Black-box wiring to the National Instruments (b) The response of the blackbox to 2 step inputs
LabPC terminator.

Figure 1.4: The “Black-box”

is not ‘visible’ to the naked eye, and that we cannot manually introduce disturbances. One com-
plication you can do is to cascade two blackboxes together to modify the dynamics.

E LECTRO - MAGNETIC BALANCE ARM

The electromagnetic balance arm shown in Fig. 1.5(a) is a fast-acting, highly oscillatory, plant
with little noise. The aim is to accurately weigh small samples by measuring the current required
to keep the balance arm level, or alternatively just to position the arm at different angles. The
output response to a step in input shown in Fig. 1.5(b) indicates how long it would take for the
oscillations to die away.
8 CHAPTER 1. INTRODUCTION

Step response of the balance arm

1

0.8

0.6

arm position
0.4

0.2

−0.2

−0.4
20 30 40 50 60 70 80 90 100

0.3

0.2

0.1
input

−0.1

−0.2

20 30 40 50 60 70 80 90 100
time ∆T=0.05 seconds

(a) The electromagnetic balance arm (b) The response of the arm to a step changes in input.

Figure 1.5: The electromagnetic balance arm

F LAPPER

Contrary to the balance arm, the flapper in Fig. 1.6 and Fig. 1.7(b) has few dynamics, but signif-
icant low-pass filtered measurement noise. An interesting exercise is to place two flapper units
in close proximity. The air from one then disturbs the flapper from the other which makes an
interacting multivariable plant.

S TEPPER MOTORS

A stepper motor is an example of a totally discrete system.

1.4.2 M ULTI - INPUT AND MULTI - OUTPUT PLANTS

It is possible to construct multivariable interacting plants by physically locating two plants close
to each other. One possibility is to locate two flappers adjacent to each other, another possibility
is one flapper and one balance arm. The extent of the interaction can be varied by adjusting the
relative position of the two plants.

H ELICOPTER

The model helicopter, Fig. 1.8, is an example of a highly unstable, multivariable (3 inputs, 2
outputs) nonlinear, strongly interacting plant. It is a good example of where we must apply
control (or crash and burn). Fig. 1.9(b) shows the controlled response using 2 PID controllers to
control the direction and altitude. Fig. 1.10(a) shows a 3 dimensional view of the desired and
actual flight path. Section 10.2.4 illustrates the control of a different three degree of freedom
helicopter.
1.4. LABORATORY EQUIPMENT FOR CONTROL TESTS 9

position transducer
✙

❖
✲ angle
✲ Flapper arm
Motor ✲
✲
✻ ✻
❄

✻ ✻
power to the motor

auto/man
✌ ◆

✻ ✻ switch On/Off
manual

❄ ❄
Digital to Analogue
A/D
(Fan speed)

Figure 1.6: The flapper wiring

Step response of the flapper

0.35

0.3

0.25
output

0.2

0.15

0.1
0 5 10 15 20 25 30 35 40 45 50

0.4

0.3
input

0.2

0.1

−0.1
0 5 10 15 20 25 30 35 40 45 50
time (sec) dt=0.05

(a) The fan/flapper equipment (b) Step response of the flapper

Figure 1.7: The Fan and flapper.

 10 CHAPTER 1. INTRODUCTION

control inputs
z }| {
top rotor

side rotor

✲ ✶ moveable counter weight

 ✙
 azimuth
plant outputs
 ✲

 ✮
elevation angle
support stand

Figure 1.8: Helicopter plant with 2 degrees of freedom. See also Fig. 1.9(a).

1.5 S LOWING DOWN S IMULINK

For some applications like the development of PID controllers you want to be able to slow down
the S IMULINK simulation to have time to manually introduce step changes, add disturbances,
switch from automatic to manual etc. If left alone in simulation mode, S IMULINK will run as fast
as possible, but it will slow down when it needs to sub-sample the integrator around discontinu-
ities or periods when the system is very stiff.

The S IMULINK Execution Control block allows you to specify that the simulation runs at a multi-
ple of real-time. This is most useful when you want to slow down a simulation, or ensure that it
runs at a constant rate.

The block is available from: http://www.mathworks.com/matlabcentral/fileexchange

The implementation is shown in Fig. 1.11 where the Simulation Execution Control block looks
like an A/D card but in fact is not connected to anything, although it does export some diagnostic
timing information.

An alternative method to slow down a Simulink simulation and force it to run at some set rate is
to use the commercial Humusoft real-time toolbox, but not with the A/D card actually interfaced
to anything.
1.5. SLOWING DOWN SIMULINK 11

(a) The ‘flying’ helicopter balanced using 2 PID controllers since it is openloop unstable and would otherwise crash.

Helicopter
1.5

0.5
output

−0.5

−1

−1.5
0 10 20 30 40 50 60 70 80 90 100

0.5
input

−0.5

−1
0 10 20 30 40 50 60 70 80 90 100
Time (sec)
(b) Multivariable PID control of the helicopter exhibiting mediocre controlled response and severe
derivative kick.

Figure 1.9: Multivariable PID control of an unstable helicopter

12 CHAPTER 1. INTRODUCTION

Helicopter flight path

0.5

0
Up/Down

−0.5

−1
1.5
1
0.5
0 250
200
−0.5
150
−1 100
50
East/West −1.5 0
time

(a) Multivariable helicopter control (b) Model helicopter in the trees

Figure 1.10: Helicopter flying results

Scope1
Simulink
Execution
Control

1
s+1
Signal Transfer Fcn Scope
Generator

Figure 1.11: Real-time Simulink simulations. A parameter in the S IMULINK E XECUTION C ON -

TROL block sets the speed of the simulation.
C HAPTER 2

F ROM DIFFERENTIAL TO DIFFERENCE

EQUATIONS

Excerpt from What a sorry state of affairs, Martin Walker, Guardian Weekly, June 29, 1997.

. . . But then he said something sensible, as he quite often does. . . .

“We need to treat individuals as individuals and we need to address discrete problems for what they are,
and not presume them to be part of some intractable racial issue.” Gingrich, properly understood, is a
national treasure, and not only because he is one of the few Americans who understand the difference
between “discreet” and “discrete”.

2.1 C OMPUTER IN THE LOOP

The cost and flexibility advantages of implementing a control scheme in software rather than
fabricating it in discrete components today are simply too large to ignore. However by inserting a
computer to run the software necessitates that we work with discrete regularly sampled signals.
This added complexity, which by the way is more than compensated by the above mentioned
advantages, introduces a whole new control discipline, that of discrete control.

Fig. 2.1 shows a common configuration of the computer in the loop. For the computer to respond
to any outside events, the signals must first be converted from an analogue form to a digital
signal, say 1 to 5 Volts, which can, with suitable processing be wired to an input port of the
computer’s processor. This device that accomplishes this conversion is called an Analogue to
Digital converter or A/D. Similarly, any binary output from the pins on a processor must first be
converted to an analogue signal using a Digital to Analogue converter, or D/A. In some micro-
controllers, (rather than micro-processors), such as Intel’s 8048 or some versions of the 8051, these
converters may be implemented on the microcontroller chip.

Digital to analogue conversion is easy and cheap. One simply loads each bit across different re-
sistors, and sums the resultant voltages. The conversion is essentially instantaneous. Analogue
to digital is not nearly as easy nor cheap, and this is the reason that the common data acquisition
cards you can purchase for your PC will often multiplex the analogue input channel. There are
various schemes for the A/D, one using a D/A inside a loop using a binary search algorithm.
Obviously this conversion is not instantaneous, although this is not normally considered a prob-
lem for process control applications. Any introductory electrical engineering text such as [173]

13
14 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

Process control computer

+
✲
r(t) + ✲ A/D ✲ computer ✲ D/A ✲ Plant ✲ y(t)
setpoint ✻
− output

Figure 2.1: The computer in the control loop

will give further details on the implementation details.

2.1.1 S AMPLING AN ANALOGUE SIGNAL

It is the A/D converter that is the most interesting for our analysis of the discrete control loop.
The A/D converter will at periodic intervals defined by the computer’s clock sample the continu-
ous analogue input signal. The value obtained is typically stored in a device called a zeroth-order
hold until it is eventually replaced by the new sample collected one sample time later. Given con-
straints on cost, the A/D converter will only have a limited precision, or limited number of bits
in which to store the incoming discrete value. Common A/D cards such as the PCLabs card,
[3], use 12 bits giving 212 or slightly over 4 thousand discretisation levels. The residual chopped
away is referred to as the quantisation error. For a given number of bits, b, used in the converter,
the amplitude quantisation is
δ = 2−b
Low cost analogue converters may only use 8 bits, while digital audio equipment use between
16 and 18 bits.

Fig. 2.2 shows the steps in sampling an analogue signal with a three bit (8 discrete levels), A/D
sampler. The dashed stair plot gives an accurate representation of the sampled signal, but owing
to the quantisation error, we are left with the solid stair plot. You can reproduce Fig. 2.2 in
M ATLAB using the fix command to do the chopping, and stairs to construct the stair plot.
While other types of hold are possible, anything higher than a first-order hold is rarely used.

Sampling a continuous signal

8
Analogue
7 Sampled
quantisized
6

4
Signal

Figure 2.2: Sampling an analogue signal 2

(heavy solid) with a three bit (8 discrete lev- 1

els) A/D converter and zeroth-order hold. 0

The sampled values, •, are chopped to the −1

next lowest discrete integer level giving the −2

0 2 4 6 8 10 12
sampled and quantisied output Sample time, (k∆T)
2.1. COMPUTER IN THE LOOP 15

2.1.2 S ELECTING A SAMPLE RATE

Once we have decided to implement discrete control, rather than continuous control, we must
decided on a reasonable sampling rate. This is a crucial parameter in discrete control systems.
The sample time, (T or sometimes denoted ∆t ), is measured in time units, say seconds or in
industrial applications, minutes. The reciprocal of the sample time is the sample frequency, f ,
and is usually measured in cycles/second or Hertz. The radial or angular velocity (which some
confusingly also term frequency) is denoted ω and is measured in radians/second. The inter-
relationships between these quantities are
   
1 cycles ω radians/s
f= = (2.1)
T second 2π radians/cycle

The faster the sampling rate, (the smaller the sampling time, T ), the better our discretised signal
approximates the real continuous signal. However, it is uneconomic to sample too fast, as the
computing and memory hardware may become too expensive. When selecting an appropriate
sampling interval, or sample rate, we should consider the following issues:

• The maximum frequency of interest in the signal

• The sampling theorem which specifies a lower limit required on the sampling rate to resolve
any particular frequency unambiguously. (See §2.1.3 following.)

• Any analogue filtering that may be required (to reduce the problem of aliasing)

• The cost of the hardware and the speed of the A/D converters.

Ogata discusses the selection of a sample time qualitatively in [150, p38]. However for most
chemical engineering processes, which are dominated by relatively slow and overdamped pro-
cesses, the sample time should lie somewhere between ten times the computational delay of the
hardware tc and some small fraction of the process dominant time constant τ , say
τ
10tc ≤ T ≤ (2.2)
10
For most chemical engineering applications, the computational delay is negligible compared with
the process time constant, (tc → 0), so we often choose T ≈ τ /10. Thus for a simple first order
rise, we would expect to have about 20–30 data samples from 0 to 99%. Some may argue that
even this sampling rate is too high, and opt for a more conservative (larger) sample time down
to τ /6. Note that commonly used guidelines such as presented in Table 22.1 in [182, p535] span a
wide range of recommended sample times.

O VERLY FAST SAMPLING

Apart from the high cost of fast A/D converters, there is another argument against fast sampling.
When one samples a continuous system, the poles of a stable continuous system map to the poles
of a stable discrete system as T goes to zero. However the zeros in the LHP of a continuous
system may not map to zeros inside the unit disk of a discrete system as T tends to zero. This
nonintuitive result could create problems if one samples a system, and then uses the inverse
of this system within a one step controller, since now the zeros outside the unit circle become
unstable poles inside the controller for small sample times. Note that the inverse continuous
system is stable, and the discrete inverse system will be stable for large sample times, and it will
only be unstable for small sample times.
16 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

2.1.3 T HE SAMPLING THEOREM AND ALIASES

“Dinsdale was a gentleman. And what’s more he knew how to treat a female impersonator”
John Cleese in Monty Python (15-9-1970)

The sampling theorem gives the conditions necessary to ensure that the inverse sampling pro-
cedure is a “one to one” relationship. To demonstrate the potential problems when sampling,
consider the case where we have two sinusoidal signals, but at different frequencies.1
   
7 1
y1 = − sin 2π t and y2 = sin 2π t
8 8

If we sample these two signals relatively rapidly, say at 25Hz (T =40ms) then we can easily see
two distinct sine curves. However if we sample at T = 1, we obtain identical points.

Aliasing
1

1 t = [0:0.04:8]'; 0.5
y1= -sin(2*pi*7/8*t);
signal

y2= sin(2*pi/8*t);
plot(t,[y1, y2]) 0

As shown opposite, we note that y1 −0.5

(solid) and y2 (dashed) just happen to
coincide at t = 0, 1, 2, 3, . . . (•). −1
0 2 4 6 8
time (s)

Consequently, at the slower sampling rate of 1 Hz, we cannot distinguish between the two dif-
ferent signals. This is the phenomenon known as aliasing since one of the frequencies “pretends”
to be another. In conclusion, two specific sinusoids of different frequencies can have identical
sampled signals. Thus in the act of taking samples from a continuous measurement, we have lost
some information. Since we have experienced a problem when we sample too slowly, it is rea-
sonable to ask what the minimum rate is so that no aliasing occurs. This question is answered by
the sampling theorem which states:

To recover a signal from its sample, one must sample at least two times a
period, or alternatively sample at a rate twice the highest frequency of interest
in the signal.

Alternatively, the highest frequency we can unambiguously reconstruct for a given sampling
rate, 1/T , is half of this, or 1/(2T ). This is called the Nyquist frequency, fN .

In the second example above when we sampled at 1 Hz, the sampling radial velocity was ωs =
2π = 6.28 rad/s. This was satisfactory to reconstruct the low frequency signal (f1 = 1/8 Hz) since
2ω1 = 1.58 rad/s. We are sampling faster than this minimum, so we can reconstruct this signal.
However for the faster signal (f2 = 7/8 Hz), we cannot reconstruct this signal since 2ω2 = 11.0
rad/s, which is faster than the sampling radial velocity.

1 This example was adapted from Franklin & Powell p81

2.1. COMPUTER IN THE LOOP 17

2.1.4 D ISCRETE FREQUENCY

If we sample the continuous signal,

x(t) = A cos(ωt)
with a sample time of T ,
x(nT ) = A cos(ωnT ) = A cos(Ωn)
where the digital frequency, Ω, is defined as

def f
Ω = ωT = 2πf T = 2π
fs

then the range of analogue frequencies is 0 < f < ∞ while the range of digital frequencies is
limited by the Nyquist sampling limit, fs /2 giving the allowable range for the digital frequency
as
0≤Ω≤π

M ATLAB unfortunately decided on a slightly different standard in the S IGNAL P ROCESSING tool-
box. Instead of a range from zero to π, M ATLAB uses a range from zero to 2 where 1 corresponds
to half the sampling frequency or the Nyquist frequency. See [42].

In summary:

symbol units
sample time T or ∆t s
sampling frequency fs = T1 Hz
angular velocity ω = 2πf rad/s
digital frequency Ω = ωT = 2πf
fs –

where the allowable ranges are:

0 ≤ω < ∞, continuous
0 ≤Ω ≤ π, sampled

and Nyquist frequency, fN , and the dimensionless Nyquist frequency, ΩN are:

fs 1
fN = = [Hz]
2 2T
2πfN
ΩN = =π
fs

It is practically impossible to avoid aliasing problems when sampling, using only digital filters.
Almost all measured signals are corrupted by noise, and this noise has usually some high or even
infinite frequency components. Thus the noise is not band limited. With this noise, no matter
how fast we sample, we will always have some reflection of a higher frequency component that
appears as an impostor or alias frequency.

If aliasing is still a problem, and you cannot sample at a higher rate, then you can insert a low
pass analogue filter between the measurement and the analogue to digital converter (sampler).
The analogue filter or in this case known as an anti-aliasing filter, will band-limit the signal, but
not corrupt it with any aliasing. Expensive high fidelity audio equipment will still use analogue
filters in this capacity. Analogue and digital filters are discussed in more detail in chapter 5.
18 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

Detecting aliases

Consider the trend y(t) in Fig. 2.3 where we wish to estimate the important frequency compo-
nents of the signal. It is evident that y(t) is comprised of one or two dominating harmonics.

3
continuous
Ts=0.7
2
Ts=1.05

output
0

Figure 2.3: Part of a noisy time series with

−1
unknown frequency components. The con-
tinuous underlying signal (solid) is sam- −2
0 1 2 3 4 5 6 7
pled at T = 0.7, (•), and at T = 1.05 s (△). time (s)

The spectral density when sampling at T = 0.7s (• in Fig.2.3) given in the upper trend of Fig. 2.4
exhibits three distinct peaks. These peaks are the principle frequency components of the signal
and are obtained by plotting the absolute value of the Fourier transform of the time signal2 ,
|DF T {y(t)}|. Reading off the peak positions, and for the moment overlooking any potential

2
10 Ts = 0.7 s
power

0
10
fN=0.71

2
10 Ts = 1.05 s
power

0
10
fN=0.48

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

frequency (Hz)

Figure 2.4: The frequency component of a signal sampled at Ts = 0.7s (upper) and Ts = 1.05s
(lower). The Nyquist frequencies for both cases are shown as vertical dashed lines. See also
Fig. 2.5.

problems with undersampling, we would expect y(t) to be something like

y(t) ≈ sin(2π0.1t) + sin(2π0.5t) + sin(2π0.63t)
However in order to construct Fig. 2.4 we had to sample the original time series y(t) possibly
introducing spurious frequency content. The Nyquist frequency fN = 1/(2Ts ) is 0.7143 and is
shown as a vertical dashed line in Fig. 2.4(top). The power spectrum is reflected in this line, but
is not shown in Fig. 2.4.

If we were to re-sample the process at a different frequency and re-plot the power density plot
then the frequencies that were aliased will move in this second plot. The △ points in Fig. 2.3 are
sampled at ∆t = 1.05s with corresponding spectral power plot in Fig. 2.4. The important data
from Fig. 2.4 is repeated below.
2 More about the spectral analysis or the power spectral density of signals comes in chapter 5.
2.2. FINITE DIFFERENCE MODELS 19

Curve Ts (s) fN (Hz) peak 1 peak 2 peak 3

top 0.7 0.7143 0.1 0.50 0.63
bottom 1.05 0.4762 0.1 0.152 0.452

First we note that the low frequency peak (f1 = 0.1Hz) has not shifted from curve a (top) to
curve b (bottom), so we would be reasonably confident that f1 = 0.1Hz and is not corrupted by
the sampling process.

However, the other two peaks have shifted, and this shift must be due to the sampling process.
Let us hypothesize that f2 = 0.5Hz. If this is the case, then it will appear as an alias in curve b
since the Nyquist frequency for curve b (fN (b) = 0.48) is less than the proposed f2 = 0.5, but it
will appear in the correct position on curve a. The apparent frequency fˆ2 (b) on curve b will be

fˆ2 (b) = 2fN (b) − f2 = 2 × 0.4762 − 0.5 = 0.4524

which corresponds to the third peak on curve b. This would seem to indicate that our hypothesis
is correct for f2 . Fig. 2.5 shows this reflection.

2
10 T = 0.7 s
s
power

0
10
f =0.71
N

2
10 T = 1.05 s
s
Figure 2.5: Reflecting the frequency re-
power

sponse in the Nyquist frequency from

0
10 Fig. 2.4 shows why one of the peaks is ev-
f =0.48
N
ident in the frequency spectrum at Ts =
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
frequency (Hz) 1.05, but what about the other peak?

Now we turn our attention to the final peak. f3 (a) = 0.63 and f2 (b) = 0.152. Let us again try
the fact that f3 (a) is the true frequency = 0.63Hz. If this were the case the apparent frequency
on curve b would be fˆ3 (b) = 2fN (b) − f2 = 0.3224Hz. There is no peak on curve b at this
frequency so our guess is probably wrong. Let us suppose that the peak at f3 (a) is the first
harmonic. In that case the true frequency will be f3 = 2fN (a) − fˆ3 (a) = 0.8Hz. Now we check
using curve b. If the true third frequency is 0.8Hz, then the apparent frequency on curve b will
be fˆ3 (b) = 2fN (b) − f3 = 0.153Hz. We have a peak here which indicates that our guess is a good
one. In summary, a reasonable guess for the unknown underlying function is

y(t) ≈ sin(2π0.1t) + sin(2π0.5t) + sin(2π0.8t)

although we can never be totally sure of the validity of this model. At best we could either re-
sample at a much higher frequency, say fs > 10Hz, or introduce an analogue low pass filter to cut
out, or at least substantially reduce, any high frequencies that may be reflected.

2.2 F INITE DIFFERENCE MODELS

To compute complex mathematical operations such as differentiation or integration, on a digital

computer, we must first discretise the continuous equation into a discrete form. Discretising is
where we convert a continuous time differential equation into a difference equation that is then
possible to solve on a digital computer.
20 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

In control applications, we are often required to solve

ordinary differential equations such as
dy
= f (t, y) (2.3)
dt
for y(t). The derivative can be approximated by the sim-
ple backward difference formula
dy yt − yt−T
≈ (2.4)
dt T
where T is the step size in time or the sampling inter-
val. Provided T is kept sufficiently small, (but not zero) Figure 2.6: In Nathaniel Hawthorne’s The
the discrete approximation to the continuous derivative Scarlet Letter, Hester Prynne is forced to
is reasonable. Inserting this approximation, Eqn. 2.4 wear a large, red “A” on her chest when
into our original differential equation, Eqn. 2.3, and re- she is found guilty of adultery and refuses
arranging for yt results in the famous Euler backward to name the father of her illegitimate child.
difference scheme;
yt = yt−T + T f (t − T, yt−T ) (2.5)
Such a scheme is called a recurrence scheme since a previous solution, yt−T , is used to calculate
yt , and is suitable for computer implementation. The beauty of this crude method of differen-
tiation is the simplicity and versatility. Note that we can discretise almost any function, linear
or nonlinear without needing to solve analytically the original differential equation. Of course
whether our approximation to the original problem is adequate is a different story and this com-
plex, but important issue is addressed in the field of numerical analysis.

Problem 2.1 1. What is the finite difference approximation for a third order derivative?
2. Write down the second order finite difference approximation to
d2 y dy
τ2 2
+ 2ζτ + y = ku
dt dt

2.2.1 D IFFERENCE EQUATIONS

Difference equations are the discrete counterpart to continuous differential equations. Often the
difference equations we are interested in are the result of the discretisation of a continuous dif-
ferential equation, but sometimes, as in the case below, we may be interested in the difference
equation in its own right. Difference equations are much easier to simulate and study in a digital
computer than a differential equation for example, since they are already written in a form where
we can just step along following a relation
xk+1 = f (xk )
rather than resorting to integration. This step by step procedure is also known as a mathematical
mapping since we map the old vector xk to a new vector xk+1 .

H ÉNON ’ S CHAOTIC ATTRACTOR

The system of discrete equations known as Hénon’s attractor,

xk+1 = (yk + 1) − 1.4x2k , x0 = 1 (2.6)
yk+1 = 0.3xk , y0 = 1 (2.7)
2.3. THE Z TRANSFORM 21

is an interesting example of a discrete mapping which exhibits chaotic behaviour.

We will start from point x0 = (x0 , y0 ) = (1, 1) and investigate the subsequent behaviour using
S IMULINK. Since Hénon’s attractor is a discrete system, we will use unit delays to shift back
from xk+1 → xk and to better visualise the results we will plot an (x, y) phase plot. The S IMULINK
simulation is given in Fig. 2.7. Compare this S IMULINK construction (which I have drawn de-
liberately to flow from right to left) so that it matches the way we read the system equations,
Eqn. 2.6-Eqn. 2.7.

1.4
Product
Gain2

x 1/z 1
To Workspace x1 Sum1 Constant

XY Graph
y 1/z 0.3
To Workspace2 y1
Gain

Figure 2.7: Hénon’s attractor, (Eqns 2.6–2.7), modelled in S IMULINK. See also Fig. 2.8.

The time response of x and y, (left figures in Fig. 2.8), of this mapping is not particularly inter-
esting, but the result of the phase plot (without joining the dots), (right figure in Fig. 2.8), is an
interesting mapping.

Actually what the chaotic attractor means is not important at this stage, this is after all only a
demonstration of a difference equation. However the points to note in this exercise are that we
never actually integrated any differential equations, but only stepped along using the unit delay
blocks in S IMULINK. Consequently these types of simulations are trivial to program, and run
very fast in a computer. Another example of a purely discrete system is given in later Fig. 6.10.

2.3 T HE z TRANSFORM

The sampling of a continuous function f (t) at regular sampling intervals equally spaced T time
units apart creates a new function f ∗ (t),

∞
X
f ∗ (t) = f (kT ) δ(t − kT ) (2.8)
n=0

where δ(x) is the Dirac delta function which is defined as being a perfect impulse at x = 0,
and zero elsewhere. Note that the actual value of δ(0) is infinity, but the integral of δ(t) is 1.0.
Expanding the summation term in Eqn 2.8 gives

f ∗ (t) = f0 δ(t) + fT δ(t − T ) + f2T δ(t − 2T ) + f3T δ(t − 3T ) + · · ·

22 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

2
1
1 Start point
xk

−1 0.5

k
−2
1

y
0.5 0
k
y

−0.5 −0.5
0 5 10 15 20 25 30 35 −1.5 −1 −0.5 0 0.5 1 1.5
time counter x
k

(a) The time response (b) The (x, y) or phase response

Figure 2.8: The time response (left) and phase plot (right) of Hénon’s attractor as computed by
the S IMULINK block diagram in Fig. 2.7.

Now suppose we wish to know the value of the third sampled function f ∗ (t = 3T ). For simplicity
we will assume the sample time is T = 1.

f ∗ (3) = f0 δ(3) + f1 δ(2) + f2 δ(1) + f3 δ(0) + f4 δ(−1) + · · ·

| {z } | {z }
all zero all zero

All the terms except for the term containing δ(0) are zero, while δ(0) = ∞. Thus the “value”
of the function f ∗ (3) = ∞. Often you will see a graph of f ∗ (t) depicted where the heights of
the values at the sample times are the same as the height of the continuous distribution such as
sketched in Fig. 2.13. Strictly this is incorrect, as it is the ‘strength’ or integral of f ∗ (t) which is the
same as the value of the continuous distribution f (t). I think of the “function” f ∗ (t) as a series of
impulses whose integral is equal to the value of the continuous function f (t).

If we take the Laplace transform of the sampled function f ∗ (t), we get

(∞ )
X
∗
L {f (t)} = L fkT δ(t − kT ) (2.9)
n=0

Now the function fkT is assumed constant so it can be factored out of the Laplace transform,
and the Laplace transform of the δ(0) is simply 1.0. If the impulse is delayed kT units, then the
Laplace transform is e−kT × 1. Thus Eqn. 2.9 simplifies to
∞
X ∞
X
L {f ∗ (t)} = fkT e−kT s = fk z −k
n=0 k=0

where we have defined

def
z = esT (2.10)
∗
In summary, the z-transform of the sampled function f (t) is defined as
∞
X
def
F (z) = Z {f ∗ (t)} = fk z −k (2.11)
k=0

The function F (z) is an infinite series, but can be written in a closed form if f (t) is a rational
function.
2.3. THE Z TRANSFORM 23

2.3.1 z- TRANSFORMS OF COMMON FUNCTIONS

We can use the definition of the z-transform, Eqn. 2.11, to compute the z-transform of common
functions such as steps, ramps, sinusoids etc. In this way we can build for ourselves a table of
z-transforms such as those found in many mathematical or engineering handbooks.

S AMPLED STEP FUNCTION

The unit sampled step function is simply s(kT ) = 1, k ≥ 0. The z-transform of s(kT ) following
the definition of Eqn. 2.11 is
∞
X
S(z) = kn z −k = 1 + 1 · z −1 + 1 · z −2 + 1 · z −3 + · · · (2.12)
k=0

By using the sum to infinity for a geometric series,3 we obtain the closed form equivalent for
Eqn. 2.12 as
1
S(z) = (2.13)
1 − z −1

R AMP AND EXPONENTIAL FUNCTIONS

For the ramp function, x(k) = n for k ≥ 0, and the exponential function, x(k) = e−an , k ≥ 0, we
could go through the same formal analytical procedure, but in this case we can use the ztrans
command from the symbolic toolbox in M ATLAB.

 z-transform
The of the exponential function,
The z-transform of the ramp function, Z {k}, Z e−ak is
using the Symbolic toolbox for M ATLAB is
1 >> syms k T a
>> ztrans(exp(-a*k*T))
1 >> syms k T
ans =
>> ztrans(k*T)
z/exp(-a*T)/(z/exp(-a*T)-1)
ans =
>> simplify(ans)
T*z/(z-1)ˆ2
6 ans =
z*exp(a*T)/(z*exp(a*T)-1)

In a similar manner we can build our own table of z-transforms of common functions.

F INAL AND INITIAL VALUE THEOREMS

Tables and properties of z-transforms of common functions are given in many handbooks and
texts such as [150, p49,67], but two of the more useful theorems are the initial value and the final
3 The sum of a geometric series of n terms is
a (1 − r n )
a + ar + ar 2 + · · · + ar n−1 = ; r 6= 1
1−r
and if |r| < 1, then the sum for an infinite number of terms converges to
a
S∞ =
1−r
24 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

value theorems given in Table 2.1. As in the continuous case, the final value theorem is only
applicable if the transfer function is stable.

Table 2.1: Comparing final and initial value theorems for continuous and discrete systems

Continuous Discrete
Initial value, limt→0 y(t) lims→∞ sY (s) lim
 z→∞ Y (z)


Final value, limt→∞ y(t) lims→0 sY (s) limz→1 1 − z −1 Y (z)

2.4 I NVERSION OF z- TRANSFORMS

The usefulness of the z-transform is that we can do algebraic manipulations on discrete difference
equations in the same way we can do algebraic manipulations on differential equations using the
Laplace transform. Generally we:

1. Convert to z-transforms, then

2. do some relatively simple algebraic manipulations to make it easier to
3. invert back to the time domain.

The final step, the inversion, is the tricky part. The process of obtaining f ∗ (t) back from F (z) is
called inverting the z-transform, and is written as
f ∗ (t) = Z −1 {F (z)} (2.14)
∗
Note that only the sampled function is returned f (t), not the original f (t). Thus the full inver-
sion process to the continuous domain f (t) is a “one to many” operation.

There are various ways to invert z-transforms:

1. Use tables of z-transforms and their inverses in handbooks, or use a symbolic manipulator.
(See §2.4.1.)
2. Use partial-fraction expansion and tables. In M ATLAB use residue to compute the partial
fractions. (See §2.4.2.)
3. Long division for rational polynomial expressions in the discrete domain. In M ATLAB use
deconv to do the polynomial division. (See §2.4.3.)
4. Computational approach while suitable for computers, has the disadvantage that the an-
swer is not in a closed form. (See §2.4.4.)
5. Analytical formula from complex variable theory. (Not useful for engineering applications,
see §2.4.5.)

2.4.1 I NVERTING z- TRANSFORMS WITH SYMBOLICALLY

The easiest, but perhaps not the most instructive, way to invert z-transforms is to use a computer
algebra package or symbolic manipulator such as the symbolic toolbox or M U PAD. One simply
enters the z-transform, then requests the inverse using the iztrans function in a manner similar
to the forward direction shown on page 23.
2.4. INVERSION OF Z -TRANSFORMS 25

>> syms z % Define a symbolic variable z

>> G = (10*z+5)/(z-1)/(z-1/4) % Construct G(z) = (10z + 5)/(z − 1)/(z − 1/4)
3 G =
(10*z+5)/(z-1)/(z-1/4)
>> pretty(G) % check it
10 z + 5
-----------------
8 (z - 1) (z - 1/4)
>> iztrans(G) % Invert the z-transform, Z −1 {G(z)}
ans =
20*kroneckerDelta(n, 0) - 40*(1/4)ˆn + 20

The Kronecker δ function, or kroneckerDelta(n, 0) is a shorthand way of expressing piece-

wise functions in M ATLAB. The expression kroneckerDelta(n, 0) returns 1 when n = 0,
and 0 for all other values.
1
δ(n)

0.5

0
−5 0 5
n

This is a rather messy and needlessly complicated artifact due to the fact that the symbolic tool-
box does not know that n is defined as positive. Also M ATLAB automatically choose n for the
indexing variable, but we might prefer to use k and explicitly inform M ATLAB that k > 0. This
gives us a cleaner solution.

>> syms z
>> syms k positive
>> G = (10*z+5)/(z-1)/(z-1/4);
4

>> y=iztrans(G,k) % Invert to get y[k] = Z −1 {G(z)}

y =
20-40*(1/4)ˆk

9 >> limit(y,k,inf) % Compute final value: y[∞]

ans =
20

The last command computed the steady-state by taking the limit of y[k] as k → ∞.

2.4.2 T HE PARTIAL FRACTION METHOD

Inverting transforms by partial fractions is applicable for both discrete z-transforms and con-
tinuous Laplace transforms. In both cases, we find an equivalent expression for the transform
in simple terms that are summed together. Hopefully we can then consult a set of tables, such
as [150, Table 2.1 p49], to invert each term individually. The inverse of the sum, is the sum of
the inverses meaning that these expressions when summed together give the full inversion. The
M ATLAB function residue forms partial fractions from a continuous transfer function, although
the help file warns that this operation is numerically poorly conditioned. Likewise the routine
residuez (that is residue with a trailing ‘z’) is designed to extract partial fractions in a form
suitable when using z−transforms.
26 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

As an example, suppose we wish to invert

−4s2 − 58s − 24 −4s2 − 58s − 24

G(s) = =
s3 + 2s2 − 24s s(s + 6)(s − 4)
using partial fractions noting that it contains no multiple or complex roots.

It is easy to spot the pole at s = 0, and the two remaining can be found by synthetic division, or
by factorising using the roots([1 2 -24 0]) command. The partial fraction decomposition
of G(s) can be written as

−4s2 − 58s − 24 A B C
G(s) = = + +
s(s + 6)(s − 4) s s+6 s−4
where the coefficients A, B and C can be found either by equating the coefficients or using the
‘cover-up’ rule shown below.
 
−4s2 − 58s − 24  −4s2 − 58s − 24 
A= = 1, B= =3
(s + 6)(s − 4) s=0 s(s − 4) 
s=−6

−4s2 − 58s − 24 
C= = −8
s(s + 6) 
s=4

In M ATLAB we can use the residue command to extract the partial fractions.

>> B = [-4 -58 -24]; % Numerator of G(s) = (−4s2 − 58s − 24)/(s(s + 6)(s − 4))
>> A = [ 1 2 -24 0]; % Denominator of G(s)
>> [r,p,k] = residue(B,A) % Find partial fractions
4 r = % residues (top line)
3
-8
1
p = % factors or poles
9 -6
4
0
k = % No extra parts
[]

The order of the residue and pole coefficient produced by residue is the same for each, so the
partial fraction decomposition is again

−4s2 − 58s − 24 3 −8 1
G(s) = = + +
s3 − 2s − 24s s+6 s−4 s
and we can invert each term individually perhaps using tables to obtain the time domain solution

g(t) = 3e−6t − 8e4t + 1

If the rational polynomial has repeated roots or complex poles, slight modifications to the proce-
dure are necessary so you may need to consult the help file.

S PECIAL CASES FOR z- TRANSFORMS

We can invert z-transforms in a manner similar to Laplace transforms, but if you consult a table
of z-transforms and their inverses in standard control textbooks such as Table 2-1 in [151, p49],
2.4. INVERSION OF Z -TRANSFORMS 27

you will discover that the table is written in terms of factors of the form z/(z + a) or alternatively
1/(1 − z −1) rather than say 1/(z + a). This means that we should first perform the partial fraction
decomposition on X(z)/z, rather than just X(z) as in the continuous case. An outline to follow
is given in Algorithm 2.1.

Algorithm 2.1 Inverting z-transforms by partial fractions

To invert a z transform, X(z), using partial fractions, do the following:

1. Divide X(z) by z obtaining X(z)/z.

2. Find the partial fractions of X(z)/z perhaps using residue.
3. Multiply each partial fraction by z to obtain a fraction of the form z/(z + a)
4. Find the inverse of each fraction separately using tables, and sum together to find x(t).

S YMBOLIC PARTIAL FRACTIONS IN M ATLAB

Suppose we want to invert

10z
G(z) =
z 2 + 4z + 3
to g[k] using partial fractions. Unfortunately there is no direct symbolic partial fraction command,
but Scott Budge from Utah University found the following sneaky trick which involves first inte-
grating, then differentiating the polynomial. This round-about way works because M ATLAB in-
tegrates the expression by internally forming partial fractions and then integrating term by term.
Taking the derivative brings you back to the original expression, but now in partial fraction form.

>> syms z
2 >> G = 10*z/(zˆ2 + 4*z + 3);
>> G= diff(int(G/z)) % Extract partial fractions of G(z)/z
G =
5/(z+1)-5/(z+3)
>> G=expand(G*z)
7 G =
5z 5z
5*z/(z+1)-5*z/(z+3) % Gives us G(z)/z = z+1 − z+3

>> syms k positive; iztrans(G,z,k)

ans =
12 -5*(-3)ˆk+5*(-1)ˆk

While this strategy currently works, we should take note that in future versions of the symbolic
toolbox this may not continue to work. In the above example we also have the option of using
residuez.

2.4.3 L ONG DIVISION

Consider the special case of the rational polynomial where the denominator is simply 1,
c0 + c1 z −1 + c2 z −2 + c3 z −3 + · · ·
F (z) =
1
28 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

The inverse of this transform is trivial since the time solution at sample number k or t = kT is sim-
ply the coefficient of the term z −k . Consequently it follows that f (0) = c0 , f (T ) = c1 , f (2T ) =
c2 , . . . , f (kT ) = ck etc. Thus if we have a rational transfer function, all we must do is to synthet-
ically divide out the two polynomials to obtain a single, albeit possibly infinite, polynomial in
z.

This approach is known as “direct division” [150, p69] and in general will result in an infinite
power series in z. As such it is not a particularly elegant closed-form solution, but for most stable
systems, the power series can be truncated eventually with little error.

So to invert a z-transform using long-division, we:

1. Convert the z-transform to a (possibly infinite) power series by dividing the denominator
into the numerator using long division or deconv.
2. The coefficients of z k are the solutions at the sample times kT .

Long division by hand is tedious and error-prone, but we can use the polynomial division capa-
bilities of M ATLAB’s deconv to do this long division automatically. Since deconv only performs
“integer” division returning a quotient and remainder, we must fool it into giving us the “infi-
nite” series by padding the numerator with zeros.

To invert
z2 + z z2 + z
Y (z) = 2
= 3
(z − 1)(z − 1.1z + 1) z − 2.1z 2 + 2.1z − 1
using deconv for the long division, we pad the numerator on the right with say 5 zeros, (that is,
we multiply by z 5 ), and then do the division.
z7 + z6
= z 4 + 3.1z 3 + 4.41z 2 + 3.751z + 1.7161 + remainder
z3 − 2.1z 2 + 2.1z − 1 | {z }
Q

We are not particularly interested in the remainder polynomial. The more zeros we add, (5 in the
above example), the more solution points we get.

>> Yn = [1,1,0]; % Numerator of G(z)

>> Yd = conv([1,-1],[1 -1.1 1]); % Denominator of G(z) is (z − 1)(z 2 − 1.1z + 1)
3 >> [Q,R] = deconv([Yn,0,0,0,0,0],Yd) % Multiply by z 5 to zero pad, then do long division
Q =
1.0000 3.1000 4.4100 3.7510 1.7161

>> dimpulse(Yn,Yd,7); % Matlab's check

We can also numerically verify the inversion using dimpulse for six samples as shown in Fig. 2.9.

Problem 2.2 1. The pulse transfer function of a process is given by

Y (z) 4(z + 0.3)
= 2
X(z) z − z + 0.4
Calculate the response of y(nT ) to a unit step change in x using long division.
2. Determine the inverse by long division of G(z)
z(z + 1)
G(z) =
(z − 1)(z 2 − z + 1)
2.4. INVERSION OF Z -TRANSFORMS 29

3
Output

1
Figure 2.9: Inverting z-transforms using dimpulse for
0
the first six samples. It is possible also to construct a
discrete transfer function object and then use the generic
0 1 2 3 4 5 6
sample # impulse routine.

2.4.4 C OMPUTATIONAL APPROACH

The computational method is so-called because it is convenient for a computer solution technique
such as M ATLAB as opposed to an analytical explicit solution. In the computational approach we
convert the z-transform to a difference equation and use a recurrence solution.

Consider the following transform (used in the example from page 24) which we wish to invert
back to the time domain.
X(z) 10z + 5 10z + 5 10z −1 + 5z −2
= = 2 = (2.15)
U (z) (z − 1)(z − 0.25) z − 1.25z + 0.25 1 − 1.25z −1 + 0.25z −2
where we have coloured the input blue and the output red. (We will again use this colour con-
vention in chapter 6).

The inverse is equivalent to solving for x(t) when u(t) is an impulse function. The transform can
be expanded and written as a difference equation

xk = 1.25xk−1 − 0.25xk−2 + 10uk−1 + 5uk−2 (2.16)

Since U (z) is an impulse, U (z) = 1 which means that only u(0) = 1 while all other samples
u(k) = 0, k 6= 0. Now we substitute k = 0 to start, and note that u(k) = x(k) = 0 when k < 0 by
definition. We now have enough information to compute x(0).

x0 = 1.25x−1 − 0.25x−2 + 10u−1 + 5u−2

= 1.25 × 0 − 0.25 × 0 + 10 × 0 + 5 × 0
=0

Continuing we substitute k = 1 into Eqn. 2.16 and using our previously computed x0 and noting
that now u0 is nonzero, we can find the next term.

x1 = 1.25x0 − 0.25x−1 + 10u0 + 5u−1

= 1.25 × 0 − 0.25 × 0 + 10 × 1 + 5 × 0
= 10

Just to clarify this recurrence process further, we can try the next iteration of Eqn. 2.16 with k = 2,

x2 = 1.25x1 − 0.25x0 + 10u1 + 5u0

= 1.25 × 10 − 0.2 × 0 + 10 × 0 + 5× 1
= 17.5
30 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

Now we can use the full recurrence relation, (Eqn 2.16), to obtain the solution in a stepwise
manner to build up the solution as shown in Table 2.2. All the terms on the right hand side of
Eqn 2.16 are either known initially (u(t)), or involve past information, hence it is possible to solve.
Note that in this inversion scheme, u(t) can be any known input series.

Table 2.2: The inversion of a z transform using the computation method. Only the first 5 samples
and the final value are calculated. The final value is as calculated using the symbolic manipulator
earlier on page 25.

time (kT ) u(k) x(k)

−2 0 0
−1 0 0
0 1 0
1 0 10
2 0 17.5
3 0 19.375
4 0 19.8438
.. .. ..
. . .
∞ 0 20

M ATLAB is well suited for this type of computation. To invert Eqn. 2.15 we can use the discrete
filter command. We rewrite Eqn. 2.15 in the form of a rational polynomial in z −1 ,

b0 + b1 z −1 + b2 z −2
G(z) =
1 + a1 z −1 + a2 z −2
where the numerator and denominator are entered as row vectors in decreasing powers of z. We
then create an impulse input vector u with say 6 samples and then compute the output of the
system

>> b = [0.0, 10.0, 5.0] % Numerator (10z −1 + 5z −2 )

>> a = [1.0, -1.25, 0.25] % Denominator (1 − 1.25z −1 + 0.25z −2 )
3

>> u = [1, zeros(1,5)] % impulse input, u(k) = δ(k)

u =
1 0 0 0 0 0
>> x = filter(b,a,u) % compute y(k)
8 x =
10.0000 17.5000 19.3750 19.8438 19.9609 19.9902

which gives the same results as Table 2.2.

The control toolbox has a special object for linear time invariant (LTI) systems. We can create a
discrete system with the transfer function command, tf, and by specifying a sampling time, we
imply to M ATLAB that it is a discrete system.

1 >> sysd = tf([10 5],[1 -1.25 0.25],1) % System of interest (10z + 5)/(z 2 − 1.25z + 0.25)

Transfer function:
10 z + 5
-------------------
6 zˆ2 - 1.25 z + 0.25
2.4. INVERSION OF Z -TRANSFORMS 31

Sampling time: 1

>> y=impulse(sysd) % compute the impulse response

11 y =
0
10.0000
17.5000
19.3750
16 19.8438
19.9609

Finally, the contour integral method can also be used to invert z-transforms, [150, §3–3], but
Seborg et al maintains it is seldom used in engineering practice, [182, p571] because it is fraught
with numerical implementation difficulties. These difficulties are also present in the continuous
equivalent; some of which are illustrated in the next section.

2.4.5 N UMERICALLY INVERTING THE L APLACE TRANSFORM

Surprisingly, the numerical inversion of continuous transfer functions is considered far less impor-
tant than the computation of the inverse of discrete transfer functions. This is fortunate because
the numerical inversion of Laplace transforms is devilishly tricky. Furthermore, it is unlikely that
you would ever want to numerically invert a Laplace transform in control applications, since
most problems involver little more than a rational polynomial with a possible exponential term.
For these problems we can easily factorise the polynomial and use partial fractions or use the
step orimpulse routines for continuous linear responses.

However in the rare cases where we have a particularly unusual F (s) which we wish to convert
back to the time domain, we might be tempted to use the analytical expression for the inverse
directly
Z σ+j∞
1
f (t) = F (s) est ds (2.17)
2πj σ−j∞
where σ is chosen to be larger than the real part of any singularity of F (s). Eqn. 2.17 is sometimes
known as the Bromwich integral or Mellin’s inverse formula.

The following example illustrates a straight forward application of Eqn. 2.17 to invert a Laplace
transform. However be warned that for all but the most well behaved rational polynomial exam-
ples this strategy is not to be recommended as it results in severe numerical roundoff error due
to poor conditioning.

Suppose we try to numerically invert a simple Laplace transform,

3 3

F (s) = ←→ f (t) = 1 − e−5t
s(s + 5) 5
with a corresponding simple time solution. We can start by defining an anonymous function of
the Laplace transform to be inverted.

Fs = @(s) 3.0./(s+5)./s; % Laplace function to invert F (s) = 3/(s(s + 5)) . . .

Fsexp = @(s,t) Fs(s).*exp(s*t) % and associated integrand F (s)est

Now since the largest singularity of F (s) is 0, we can choose the contour path safely to be say
σ = 0.1. We can also approximate the infinite integration interval by reasonably large numbers.
It is convenient that the M ATLAB routine integral works directly in the complex domain.
32 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

0.7 2

0.6 1.5 n o
0.5 L−1 √ 1
1
s 2 +1

0.4
f(t)

f(t)
0.5
0.3
0
0.2 n o
3
L−1 s(s+5) −0.5
0.1

0 −1
0.02 2

1
error

error
0
0

−0.02 −1
0 2 4 6 8 0 5 10 15 20
t t

(a) Inverting
n o a benign transfer function, (b) Inverting
 a challenging transfer function,
3
L−1 s(s+5) L −1 √ 1
s2 +1

Figure 2.10: Numerically inverting Laplace transforms by direct evaluation of the Bromwich
integral. Since this strategy is suspectable to considerable numerical errors, it is not to be recom-
mended.

c = 0.1; % Value σ > than any singularities of F (s)

a = c-100j; b = c+200j; % Contour path approx: σ − j∞ to σ + j∞
3

t=linspace(0,8); % time range of interest

% Numerically approximate Eqn. 2.17.
ft = 1./(2*pi*1j)*arrayfun(@(t) integral(@(x) Fsexp(x,t),a,b), ti);
plot(t,ft) % Compare results in Fig. 2.10(a).

Due to small numerical roundoff however, the returned solution has a small imaginary compo-
nent, which we can in this instance safely ignore. In this rather benign example, where we know
the analytical solution, we can validate the accuracy as shown in Fig. 2.10(a), which it has to be
admitted, is not that wonderful. We should also note in the simplistic implementation above we
have not adequately validated that our algorithmic choices such as the finite integration limits
are appropriate, so we should treat any subsequent result with caution. In fact if you repeat this
calculation with a larger value of σ, or a smaller integration range such as say 1 − 10j to 1 + 10j
then the quality of the solution drops alarmingly. One way to improve the integration accuracy
is to use the waypoints option in the integral routine.

Over the years there have been many attempts to improve the inversion strategy as reviewed
in [1], but nearly all of the proposals run into the problem√of precision. For example, Fig. 2.11
shows the problem when attempting to invert F (s) = 1/ s2 + 1 using the Gaver-Stehfest ap-
proximation algorithm with a varying number of terms. The Gaver-Stehfest family of schemes
are well known to produce inaccurate results for underdamped systems, but when we increase
the number of terms in an attempt to improve the accuracy, the results deteriorate due to numer-
ical round off. Directly applying the Bromwich integral to this transform as shown in Fig. 2.10(b)
is no better. While the strategy suggested in [1] seemingly circumvents the problem of precision,
it does so by brute-force using multi-precision arithmetic which is hardly elegant.
2.4. INVERSION OF Z -TRANSFORMS 33

No. of terms = 8
1

0.5

−0.5
0 2 4 6 8 10
No. of terms = 16
1

0.5

−0.5
0 2 4 6 8 10
No. of terms = 20
1

0.5

−0.5
0 2 4 6 8 10
No. of terms = 26
1

0
Figure 2.11: Demonstrating the precision problems
when numerically inverting the Laplace transform using
the Gaver-Stehfest algorithm with a varying number of
−1
0 2 4 6 8 10 terms. The exact inversion is the red solid line while the
t approximate inversion is given by –◦–.

An alternative numerical algorithm is the inverse Laplace transform function invlap function
contributed by Karl Hollenbeck from the Technical University of Denmark. You can test this
routine using one of the following Laplace transform pairs in Table 2.3. The first two are standard,
those following are more unusual. A more challenging collection of test transforms is given in
[1].

Table 2.3: A collection of Laplace transform pairs suitable for testing numerical inversion strate-
gies.

Description F (s) f (t) Comment

Easy test 1/s2 t
1 1


Control TF 2 e−t − e−3t
(s + 2)(s + 3)
1
Oscillatory √ J0 (t) Bessel function, see Fig. 2.10(b) & Fig. 2.11.
2
s +1
1 −1/s √
Nasty √ e √1 cos( 4t), See Fig. 2.12.
s πt

n o
The numerical solution to L−1 √1 e−1/s is compared with the analytical solution from Table 2.3
s
in Fig. 2.12 using a modest tolerance of 10−2 . The routine splits the time axis up into decades,
34 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

and inverts each separately which is why we can see some numerical noise starting at t = 10.

1.5

n o
exp(−1/s)

f(t)
0.5 L−1 √
s

0
−0.5
10

−2
10
error

−4
10

Figure 2.12: Testing the numerical inver-

√ −6
sion of the Laplace transform e−1/s / s us- 10
0 5 10 15 20 25 30 35 40
ing the invlap routine. time

All of these numerical inversion strategies, starting from the direct application of the Bromwich
integral in Eqn. 2.17, the Gaver-Stehfest algorithm and variants and the implementation invlap
require some care to get reasonable results. Unfortunately in practical cases it can be difficult to
choose suitable algorithmic tuning constants, which mean that it is difficult to generate any sort
of error bound on the computed curve.

2.5 D ISCRETISING WITH A SAMPLE AND HOLD

In §2.1.1 we saw how the continuous analogue signal must be sampled to be digestible by the
process control computer. To analyse this sampling operation, we can approximate the real sam-
pler with an ideal sampler in cascade with a hold element such as given in Fig. 2.13.

✒ f ∗ (t)
f (t) ✲ ✲ Hold ✲ fh (t)
ideal sampler

✻ ✻ ✻
✻✻✻
✻
✻
✻

✲ ✲ ✲
analogue signal sampled signal sample & held signal

Figure 2.13: Ideal sampler and zeroth-order hold

2.5. DISCRETISING WITH A SAMPLE AND HOLD 35

When we sample the continuous function f (t), we get a sampled function f ∗ (t), that is a series
of spikes that exist only at the sampling instant as shown in the middle plot of Fig. 2.13. The
sampled function, f ∗ (t), is undefined for the time between the sampling instants which is incon-
venient given that this is the bulk of the time. Obviously one sensible solution is to retain the
last sampled value until a new one is collected. During the sample period, the sampled value
f ∗ (t) is the same as the last sampled value of the continuous function; f (kT + t) = f (kT ) where
0 ≤ t < T shown diagrammatically in the third plot in Fig. 2.13. Since the last value collected is
stored or held, this sampling scheme is referred to as a zeroth-order hold. The zeroth-order refers
to the fact that the interpolating function used between adjacent values is just a horizontal line.
Higher order holds are possible, but the added expense is not typically justified even given the
improved accuracy.

We can find the Laplace transform of the zeroth-order hold element by noting that the input is an
impulse function, and the output is a rectangular pulse of T duration, and y(t) high,

output ✻
Rectangular pulse
y(t) ✙

✛ ✲ time, t
0 sample time, T

Directly applying the definition of the Laplace transform,

Z ∞ Z T
−st
L{zoh} = y(t) e dt = yt 1 × e−st dt
0 0
1 − e−sT
= yt ×
| {zs }
zoh

we obtain the transform for the zeroth-order hold. In summary, the transformation of a continu-
ous plant G(s) with zero-order hold is
 
G(s)
G(z) = (1 − z −1 ) Z (2.18)
s

Note that this is not the same as simply the z-transform of G(s).

Suppose we wish to discretise the continuous plant

8
G(s) = (2.19)
s+3
with a zero-order hold at a sample time of T = 0.1 seconds. We can do this analytically using ta-
bles, or we can trust it to M ATLAB. First we try analytically and apply Eqn 2.18 to the continuous
process. The difficult bit is the z-transformation. We note that
   
G(s) 8 3
Z = Z (2.20)
s 3 s(s + 3)

and that this is in the form given in tables of z-transforms such as [150, p50], part of which is
repeated here
36 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

X(s) X(z)
a (1−e−aT )z −1
s(s+a) (1−z −1 )(1−e−aT z −1 )

Inverting using the tables, and substituting in Eqn 2.18 gives

   
G(s) 8 3
G(z) = (1 − z −1 ) · = (1 − z −1 ) · Z
s 3 s(s + 3)
8 (1 − e−3T )z −1
= (1 − z −1 ) ·
3 (1 − z −1 )(1 − e−3T z −1 )
−0.3 −1
8 (1 − e )z
= ·
3 (1 − e−0.3 z −1 )
0.6912z −1
= , for T = 0.1 (2.21)
1 − 0.7408z −1
We can repeat this conversion numerically using the continuous-to-discrete function, c2d. We
first define the continuous system as an LTI object, and then convert to the discrete domain using
a zeroth order hold.

>>sysc = tf(8,[1 3]) % Continuous system Gc (s) = 8/(s + 3)

3 Transfer function:
8
-----
s + 3

8 >> sysd = c2d(sysc,0.1,'zoh') % convert to discrete using a zoh

Transfer function:
0.6912
----------
13 z - 0.7408

Sampling time: 0.1

which once again is Eqn. 2.21. In fact we need not specify the ’zoh’ option when constructing a
discrete model using c2d or even in S IMULINK, since a zeroth-order hold is employed by default.

Methods for doing the conversion symbolically are given next in section 2.5.1.

2.5.1 C ONVERTING L APLACE TRANSFORMS TO z- TRANSFORMS

In many practical cases we may already have a continuous transfer function of a plant or con-
troller which we wish to discretise. In these circumstances we would like to convert from a con-
tinuous transfer function in s, (such as say a Butterworth filter) to an equivalent discrete transfer
function in z. The two common ways to do this are:

1. Analytically by first inverting the Laplace transform back to the time domain, and then take
the (forward) z-transform.
def 
G(z) = Z L−1 {G(s)} (2.22)

2. or by using an approximate method known as the bilinear transform.

2.5. DISCRETISING WITH A SAMPLE AND HOLD 37

The bilinear method whilst approximate has the advantage that it just involves a simple algebraic
substitution for s in terms of z. The bilinear transform is further covered in §2.5.2.

The formal method we have already seen previously, but since it is such a common operation, we
can write a symbolic M ATLAB script to do it automatically for us for any given transfer function.
However we should be aware whether the transformation includes the sample and hold, or if
that inclusion is left up to the user.

Listing 2.1 converts a Laplace expression to a z-transform expression without a zeroth-order hold.
Including the zeroth-order hold option is given in Listing 2.2.

Listing 2.1: Symbolic Laplace to z-transform conversion

function Gz = lap2ztran(G)
% Convert symbolically G(s) to G(z)
syms t k T z
Gz = simplify(ztrans(subs(ilaplace(G),t,k*T),k,z));
5 return

Listing 2.2: Symbolic Laplace to z-transform conversion with ZOH

function Gz = lap2ztranzoh(G)
% Convert symbolically G(s) to G(z) with a ZOH.
syms s t k T z
Gz = simplify((1-1/z)*ztrans(subs(ilaplace(G/s),t,k*T),k,z));
5 return

We can test the conversion routines in Listings 2.1 and 2.2 on the trial transfer function G(s) =
1/s2 .

>> Gz =lap2ztran(1/sˆ2); % Convert G(s) = 1/s2 to the discrete domain G(z).

>> pretty(Gz)
T z
--------
5 2
(z - 1)

>> Gz = lap2ztranzoh(1/sˆ2);% Do the conversion again, but this time include a ZOH.
>> pretty(Gz)
10 2
(z + 1) T
1/2 ----------
2
(z - 1)

2.5.2 T HE BILINEAR TRANSFORM

The analytical ‘one step back–one step forward’ procedure while strictly correct, is a little tedious,
so a simpler, albeit approximate, way to transform between the Laplace domain and the z-domain
is to use the bilinear transform method, sometimes known as Tustin’s method. By definition
def
z = esT or z −1 = e−sT , (from Eqn. 2.10). If we substituted directly natural logs would appear in
38 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

the rational polynomial in z,

ln(z)
s= (2.23)
T
making the subsequent analysis difficult. For example the resulting expression could not be
transformed into a difference equation which is what we desire.

We can avoid the troublesome logarithmic terms by using a Padé approximation for the expo-
nential term as
2 − sT
e−sT = z −1 ≈ (2.24)
2 + sT
or alternatively
2 (1 − z −1 )
s≈ (2.25)
T (1 + z −1 )
This allows us to transform a continuous time transfer function G(s) to a discrete time transfer
function G(z),
G(z) = G(s)|s= 2 1−z−1 (2.26)
T 1+z −1

Eqn. 2.26 is called the bilinear transform owing to the linear terms in both the numerator and
denominator and it has the advantage that a stable continuous time filter will be stable in the dis-
crete domain. The disadvantage is that the algebra required soon becomes unwieldy if attempted
manually. Other transforms are possible and are discussed in §4-2 p308 of Ogata [150]. Always
remember that this transformation is approximate, being equivalent to a trapezoidal integration.

Here we wish to approximately discretise the continuous plant

1
G(s) =
(s + 1)(s + 2)

at a sample time of T = 0.1 using the bilinear transform. The discrete approximate transfer
function is obtained by substituting Eqn. 2.25 for s and simplifying.

1 
G(z) = 
(s + 1)(s + 2) s= 2 1−z−1
T 1+z −1

1 1 T 2 (z + 1)2
=   =
2 1−z −1
+1 2 1−z −1
+2 2 (2T 2 + 6T + 4)z 2 + (4T 2 − 8)z + 2T 2 − 6T + 4
T 1+z −1 T 1+z −1

0.0022z 2 + 0.0043z + 0.0022

= for T = 0.1
z 2 − 1.723z + 0.74

Since this transformation is just an algebraic substitution, it is easy to execute it symbolically in

M ATLAB.

1 >> syms s T z
>> G = 1/(s+1)/(s+2);
>> Gd = simplify(subs(G,s,2/T*(1-1/z)/(1+1/z)))
Gd =
1/2*Tˆ2*(z+1)ˆ2/(2*z-2+T*z+T)/(z-1+T*z+T)
6 >> pretty(Gd)
2 2
T (z + 1)
1/2 -------------------------------------
(2 z - 2 + T z + T) (z - 1 + T z + T)
2.5. DISCRETISING WITH A SAMPLE AND HOLD 39

Alternatively we could use M ATLAB to numerically verify our symbolic solution.

>> G = zpk([],[-1 -2],1)

Zero/pole/gain:
1
-----------
5 (s+1) (s+2)

>> T=0.1;
>> Gd = c2d(G,T,'tustin')

10 Zero/pole/gain:
0.0021645 (z+1)ˆ2
---------------------
(z-0.9048) (z-0.8182)

15 Sampling time: 0.1

>> tf(Gd)

Transfer function:
0.002165 zˆ2 + 0.004329 z + 0.002165
20 ------------------------------------
zˆ2 - 1.723 z + 0.7403

Sampling time: 0.1

The bilinear command in the S IGNAL P ROCESSING toolbox automatically performs this map-
ping from s to z. This is used for example in the design of discrete versions of common analogue
filters such as the Butterworth or Chebyshev filters. These are further described in §5.2.3.

Problem 2.3 1. Use Tustin’s method (approximate z-transform) to determine the discrete time
response of
4
G(s) =
(s + 4)(s + 2)

to a unit step change in input by long division. The sample time T = 1 and solve for 7 time
steps. Compare with the exact solution.

T HE FREQUENCY RESPONSE CHARACTERISTICS OF A HOLD ELEMENT

The zeroth order hold element modifies the transfer function, so consequently influences the
closed loop stability, and frequency characteristics of the discrete process. We can investigate this
influence by plotting the discrete Bode and Nyquist plots for a sampled process with and without
a zeroth order hold.

Suppose we have the continuous process

1.57
Gp (s) = (2.27)
s(s + 1)

which is sampled at a frequency of ω = 4 rad/s. (This corresponds to a sample time of T = π/2

seconds.) First we will ignore the zeroth order hold and transform Eqn 2.27 to the z domain by
40 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

using tables such as say Table 2–1, p49 #8 in DCS.

(1 − e−aT )z −1
Gp (z) = 1.57 (2.28)
(1 − z −1 )(1 − e−aT z −1 )
1.243z
= (2.29)
(z − 1)(z − 0.208)
To get the discrete model with the zeroth order hold, we again use Eqn 2.18, and the tables. Doing
this, in a similar manner to what was done for Eqn 2.21, we get
1.2215z + 0.7306
Gh0 Gp (z) = (2.30)
(z − 1)(z − 0.208)
Now we can plot the discrete frequency responses using M ATLAB to duplicate the figure given
in [111, p649] as shown in Fig. 2.14.

In the listing below, I first construct the symbolic discrete transfer function, substitute the sam-
pling time, then convert to a numeric rational polynomial. This rational polynomial is now in a
suitable form to be fed to the control toolbox routines such as bode.

syms s T z
2 G = 1.57/(s*(s+1))
Gd = lap2ztran(G) % convert to a z−transform without ZOH
Gd2 = subs(Gd,T,pi/2)
[num,den] = numden(Gd2) % extract top & bottom lines

7 B= sym2poly(num); A = sym2poly(den); % convert to polynomials

Gd = tf(B,A,pi/2) % construct a transfer function

To generate the discretisation including the zeroth-order hold, we could use the symbolic lap2ztranzoh
routine given in Listing 2.2, or we could use the built-in c2d function.

1 [num,den] = numden(G) % extract numerator & denominator

Bc= sym2poly(num); Ac = sym2poly(den);
Gc = tf(Bc,Ac)

Gdzoh = c2d(Gc,pi/2,'zoh')
6

bode(Gc,Gd,Gdzoh)
legend('Continuous','No ZOH','with ZOH')

In this case I have used vanilla Bode function which will automatically recognise the difference
between discrete and continuous transfer functions. Note that it also automatically selects both
a reasonable frequency spacing, and inserts the Nyquist frequency at fN = 1/2T = 1/π Hz or
ωN = 2 rad/s.

These Bode plots shown in Fig. 2.14, or the equivalent Nyquist plots, show that the zeroth order
hold destabilises the system. The process with the hold has a larger peak resonance, and smaller
gain and phase margins.

Of course we should compare the discrete Bode diagrams with that for the original continuous
process. The M ATLAB bode function is in this instance the right tool for the job, but I will con-
struct the plot manually just to demonstrate how trivial it is to substitute s = iω, and compute
the magnitude and phase of G(iω) for all frequencies of interest.
2.5. DISCRETISING WITH A SAMPLE AND HOLD 41

Bode Diagram

100

50
Magnitude (dB)

−50

−100
−90
Continuous
No ZOH
with ZOH
Phase (deg)

−135

−180
Figure 2.14: The Bode dia-
gram showing the difference
−225 between the continuous plant,
−2 −1 0 1 2
10 10 10 10 10 a discretisation with and with-
Frequency (rad/sec) out the zeroth-order hold.

num = 1.57; % Plant of interest G(s) = 1.57/(s2 + s + 0)

2 den = [1 1 0];
w = logspace(-2,1)'; % Frequency range 10−2 < ω < 101 rad/s.
iw = 1i*w; % s = jω
G = polyval(num,iw)./polyval(den,iw); % G(s = jω)
loglog(w,abs(G)) % Plot magnitude |G(iω)|, see Fig. 2.15.
7 semilogx(w,angle(G)*180/pi) % and phase φ(iω)

We can compare this frequency response of the continuous system with the discretised version
in Fig. 2.14.

5
10

0
AR

10

−5
10
−2 −1 0 1
10 10 10 10

−50
Phase lag φ [deg]

−100

−150

−200
Figure 2.15: A frequency response plot of
10
−2
10
−1
10
0
10
1 G(s) where s = jω constructed manually
frequency ω [rad/s] without using the bode command.
42 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

2.6 D ISCRETE ROOT LOCUS DIAGRAMS

The root locus diagram is a classical technique used to study the characteristic response of trans-
fer functions to a single varying parameter such as different controller gains. Traditionally the
construction of root locus diagrams required tedious manual algebra. However now modern
programs such as M ATLAB can easily construct reliable root locus diagrams which makes the
design procedure far more attractable.

For many controlled systems, the loop is stable for small values of controller gain Kc , but unstable
for a gain above some critical value Ku . It is therefore natural to ask how does the stability (and
response) vary with controller gain? The root locus gives this information in a plot form which
shows how the poles of a closed loop system change as a function of the controller gain Kc .
Given a process transfer function G and a controller Kc Q, the closed loop transfer function is the
familiar
Kc Q(z)G(z)
Gcl = (2.31)
1 + Kc Q(z)G(z)

where the stability is dependent on the denominator (characteristic equation) 1 + Kc Q(z)G(z),

which is of course itself dependent on the controller gain Kc . Plotting the roots of

1 + Kc Q(z)G(z) = 0 (2.32)

as a function of Kc creates the root locus diagram.

In this section4 we will plot the discrete root locus diagram for the process

0.5(z + 0.6)
Gp (z) = (2.33)
z 2 (z − 0.4)

which is a discrete description of a first order process with dead time. We also add a discrete
integral controller of the form
z
Gc (z) = Kc
z−1
We wish to investigate the stability of the characteristic equation as a function of controller gain.

With no additional arguments apart from the transfer function, rlocus will draw the root locus
selecting sensible values for Kc automatically. For this example, I will constrain the plots to have
a square aspect ratio, and I will overlay a grid of constant shape factor ζ and natural frequency
ωn using the zgrid function.

Q = tf([1 0],[1 -1],-1); % controller without gain, Q(z)

G = tf(0.5*[1 0.6],[1 -0.4 0 0],-1); % plant, Gp (z)
3

Gol = Q*G; % open loop Q(z)Gp (z)

zgrid('new'); xlim([-1.5 1.5]); axis equal
rlocus(Gol) % plot root locus, see Fig. 2.16.
rlocfind(Gol)

Once the plot such as Fig. 2.16 is on the screen, I can use rlocfind interactively to establish the
values of Kc at different critical points in Fig. 2.16. Note that the values obtained are approximate.

4 Adapted from [71, p124]

2.6. DISCRETE ROOT LOCUS DIAGRAMS 43

Pole description Name ≈ Kc

border-line stability ultimate gain 0.6
cross the ζ = 0.5 line design gain 0.2
critically damped breakaway gain 0.097
overdamped, ωn = 36◦ 0.09

Root Locus

1.5

1
0.5π/T Ultimate
0.6π/T 0.4π/T gain
0.7π/T 0.1
0.3π/T
0.2
0.3
0.8π/T 0.4 0.2π/T Design gain
0.5 0.5
0.6
0.7
0.9π/T 0.8 0.1π/T
Imaginary Axis

0.9
π/T
0
π/T

0.9π/T 0.1π/T
breakaway
−0.5 0.8π/T 0.2π/T gain

0.7π/T 0.3π/T
0.6π/T 0.4π/T
0.5π/T
−1

−1.5

−1.5 −1 −0.5 0 0.5 1 1.5

Real Axis

Figure 2.16: The discrete root locus of Eqn. 2.33 plotted using rlocus. See also Fig. 2.17 for the
resultant step responses for various controller gains.

Once I select the gain that corresponds to the ζ = 0.5 crossing, I can simulate the closed loop step
response.

Gol.Ts = 1; % Set sampling time T = 1

Kc = 0.2 % A gain where we expect ζ = 0.5
3 step(Kc*Gol/(1 + Kc*Gol),50) % Simulate closed loop, see Fig. 2.17.

A comparison of step responses for various controller gains is given in Fig. 2.17. The curve in
Fig. 2.17 with Kc = 0.2 where ζ ≈ 0.5 exhibits an overshoot of about 17% which agrees well with
the expected value for a second order response with a shape factor of ζ = 0.5. The other curves
do not behave exactly as one would expect since our process is not exactly a second order transfer
function.
44 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

1
Kc=0.1
0
2

1
Kc=0.2
0
2

Figure 2.17: Various discrete closed loop re-

sponses for gains Kc = 0.6, 0.2 and 0.1. Note that 1
with Kc = 0.6 we observe a response on the sta- Kc=0.6
bility boundary, while for Kc = 0.1 we observe
0
a critically damped response as anticipated from 0 10 20 30 40 50
the root locus diagram given in Fig. 2.16. time [s]

2.7 M ULTIVARIABLE CONTROL AND STATE SPACE ANALYSIS

In the mid 1970s, the western world suffered an ‘oil shock’ when the petroleum producers and
consumers alike realised that oil was a scarce, finite and therefore valuable commodity. This had
a number of important implications, and one of these in the chemical processing industries was
the increased use of integrated heat plants. This integration physically ties separate parts of the
plant together and demands a corresponding integrated or ‘plant-wide’ control system. Clas-
sical single-input/single-output (SISO) controller design techniques such as transfer functions,
frequency analysis or root locus were found to be deficient, and with the accompanying devel-
opment of affordable computer control systems that could administer hundreds of inputs and
outputs, many of which were interacting, nonlinear and time varying, new tools were required
and the systematic approach offered by state-space analysis became popular.

Many physical processes are multivariable. A binary distillation column such as given in Fig. 2.18
typically has four inputs; the feed flow and composition and reflux and boil-up rates, and at least
four outputs, the product flow and compositions. To control such an interacting system, mul-
tivariable control such as state space analysis is necessary. The Wood-Berry distillation column
model (Eqn 3.24) and the Newell & Lee evaporator model are other examples of industrial pro-
cess orientated multivariable models.

State space analysis only considers first order differential equations. To model higher order sys-
tems, one needs only to build systems of first order equations. These equations are conveniently
collected, if linear, in one large matrix. The advantage of this approach is a compact represen-
tation, and a wide variety of good mathematical and robust numerical tools to analyse such a
system.

A FEW WORDS OF CAUTION

The state-space analysis is sometimes referred to as the “modern control theory”, despite the fact
that is has been around since the early 1960s. However by the end of the 1970s, the promise of
2.7. MULTIVARIABLE CONTROL AND STATE SPACE ANALYSIS 45

Feed composition
✲ ✲ distillation composition, xD
(
Disturbances
✲ ✲ distillate rate, D
Feed rate
Distillation ✲ tray 5 temperature, T5
column ✲ tray 15 temperature, T15
reflux(rate
Manpipulated ✲ ✲ Bottoms rate, B
variables ✲ ✲ Bottoms comp., xB
reboiler heat

Inputs Outputs

Figure 2.18: A binary distillation column with multiple inputs and multiple outputs

the academic advances made in the previous decade were turning out to be ill-founded, and it
was felt that this new theory was ill-equipped to cope with the practical problems of industrial
control. In many industries therefore, ‘Process Control’ attracted a bad smell. Writing a decade
still later in 1987, Morari in [138] attempts to rationalise why this disillusionment was around at
the time, and whether the subsequent decade of activity alleviated any of the concerns. Morari
summarises that commentators such as [70] considered theory such as linear multivariable con-
trol theory, (i.e. this chapter) which seemed to promise so much, actually delivered very little
‘and had virtually no impact on industrial practice’. There were other major concerns such as
the scarcity of good process models, the increasing importance of operating constraints, operator
acceptance etc, but the poor track record of linear multivariable control theory landed top billing.
Incidentally, Morari writing almost a decade later still in 1994, [141], revisits the same topic giving
the reader a nice linear perspective.

2.7.1 S TATES AND STATE SPACE

State space equations are just a mathematical equivalent way of writing the original common dif-
ferential equation. While the vector/matrix construction of the state-space approach may initially
appear intimidating compared with high order differential equations, it turns out to be more con-
venient and numerically more robust to manipulate these equations when using programs such
as M ATLAB. In addition the state-space form can be readily expanded into multivariable systems
and even nonlinear systems.

The state of a system is the smallest set of values such that knowledge of these values, knowledge
of any future inputs and the governing dynamic equations, it is possible to predict everything
about the future (and past) output of the system. The state variables are often written as a column
vector of length n and denoted x. The state space is the n dimensional coordinate space that all
possible state vectors must lie in.

The input to a dynamic system is the vector of m input variables, u, that affect the state vari-
ables. Historically control engineers further subdivided the input variables into those they could
easily and consciously adjust (known as control or manipulated variable inputs), and those vari-
ables that may change, but are outside the engineer’s immediate sphere of influence (like the
weather) known as disturbance variables. The number of input variables need not be the same as
the number of state variables, and indeed m is typically less than n.
46 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

Many mathematical control text books follow a standard nomenclature. Vectors are written using
a lower case bold font such as x, and matrices are written in upper case bold such as A.

The state space equations are the n ordinary differential equations that relate the state derivatives
to the state themselves, the inputs if any, and time. We can write these equations as
ẋ = f (x, u, t) (2.34)

where f (·) is a vector function meaning that both the arguments and the result are vectors. In
control applications, the input is given by a control law which is often itself a function of the
states,
u = h(x) (2.35)
which if we substitute into Eqn. 2.34, we get the closed loop response
ẋ = f (x, h(x), t)
For autonomous closed loop systems there is no explicit dependence on time, so the differential
equation is simply
ẋ = f (x) (2.36)
For the continuous linear time invariant case, Eqn. 2.34 simplifies to
ẋ = Ax + Bu (2.37)
or in discrete form at time t = kT
xk+1 = Φxk + ∆uk (2.38)
where the state vector x is an (n × 1) vector, or alternatively written as x ∈ ℜn , the control input
vector has m elements, or u ∈ ℜm . The model (transition) matrix is (n × n) or A, Φ ∈ ℜn×n , and
the control or input matrix is (n × m) or B, ∆ ∈ ℜn×m .

Block diagrams of both the continuous and discrete formulations of the state-space model are
shown in Fig. 2.19 and both are conceptually very similar. Such a form is suitable to implement
state-space systems in a simulator such as S IMULINK for example. Fig. 2.19 also introduces a
convention we will use later in system identification, where I colour the input signals blue, and
the states red.

ẋ R x(t) xk xk+1
u(t) ✲ B ✲+ ✲ ✲ uk ✲ ∆ ✲+ ✲ z −1 ✲
✻ ✻

✛ A ✛ ✛ Φ ✛

Continuous plant Discrete plant

Figure 2.19: A block diagram of a state-space dynamic system, (a) continuous system: ẋ = Ax +
Bu, and (b) discrete system: xk+1 = Φxk + ∆uk . (See also Fig. 2.20.)

The output or measurement vector, y, are the variables that are directly measured from the operat-
ing plant, since often the true states cannot themselves be directly measured. These outputs are
related to the states by
y = g(x) (2.39)
2.7. MULTIVARIABLE CONTROL AND STATE SPACE ANALYSIS 47

If the measurement relation is linear, then

y = Cx (2.40)

where the r element measurement vector is ℜr , and the measurement matrix is sized C ∈ ℜr×n .

Sometimes the input may directly affect the output bypassing the process, in which case the full
linear system in state space is described by,

ẋ = Ax + Bu
(2.41)
y = Cx + Du

Eqn. 2.41 is the standard starting point for the analysis of continuous linear dynamic systems, and
are shown in block diagram form in Fig. 2.20. Note that in the diagram, the internal state vector
(bold lines), has typically more elements than either the input or output vectors (thin lines). The
diagram also highlights the fact that as an observer to the plant, we can relatively easily measure
our outputs (or appropriately called measurements), we presumably know our inputs to the
plant, but we do not necessarily have access to the internal state variables since they are always
contained inside the dashed box in Fig. 2.20. Strategies to estimate these hidden internal states is
known as state estimation and are described in section 9.5.

✲ D

ẋ R x(t) ❄
u(t) ✲✲ B ✲+ ✲ ✲ C ✲+ ✲ y(t)
✻

✛ A ✛

Figure 2.20: A complete block diagram of a continuous state-space dynamic system with outputs
and direct measurement feed-through, Eqn. 2.41.

2.7.2 C ONVERTING DIFFERENTIAL EQUATIONS TO STATE - SPACE FORM

If we have a collection or system of interlinked differential equations, it is often convenient and

succinct to group them together in a vector/matrix collection. Alternatively if we start with an
nth order differential equation, for the same reasons it is advisable to rewrite this as a collection
of n first order differential equations known as Cauchy form.

Given a general linear nth order differential equation

y (n) + a1 y (n−1) + a2 y (n−2) + · · · + an−1 ẏ + an y = b0 u(n) + b1 u(n−1) + · · · + bn−1 u̇ + bn u (2.42)

by inspection, we can create the equivalent transfer function as

Y (s) b0 sn + b1 sn−1 + · · · + bn−1 s + bn

= n (2.43)
U (s) s + a1 sn−1 + · · · + an−1 s + an
48 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

We can cast this transfer function (or alternatively the original differential equation) into state
space in a number of equivalent forms. They are equivalent in the sense that the input/output
behaviour is identical, but not necessarily the internal state behaviour.

The controllable canonical form is

      
ẋ1 0 1 0 ··· 0 x1 0
 ẋ2   0
   0 1 ··· 0 


 x2  
  0 

 ..   .. .. .. .. ..   ..   .. 
 . = . . . . .   . + . u (2.44)
       
 ẋn−1   0 0 0 ··· 1   xn−1   0 
ẋn −a0 −a1 −a2 · · · −an−1 xn 1
 
x1
h

i x2 

.. .. .  .. 
y= b n − an b 0 . bn−1 − an−1 b0 . · · · .. b1 − a1 b0 
 .  + b0 u

(2.45)
 xn−1 
xn
which is useful when designing pole-placement controllers. The observable canonical form is
      
ẋ1 0 0 · · · 0 −an x1 b n − an b 0
 ẋ2   1 0 · · · 0 −an−1   x2   bn−1 − an−1 b0 
      
 ..  =  .. . . .. .. ..   ..  +  .. u (2.46)
 .   . . . . .   .   . 
ẋn 0 0 · · · 1 −a1 xn b 1 − a1 b 0
 
x1
 x2 
 
y = 0 0 · · · 0 1  ...  + b0 u
  
(2.47)
 
 xn−1 
xn
and is useful when designing observers.

If the transfer function defined by Eqn. 2.43 has real and distinct factors,
Y (s) b0 sn + b1 sn−1 + · · · + bn−1 s + bn
=
U (s) (s + p1 )(s + p2 ) · · · (s + pn )
c1 c2 cn
= b0 + + + ···+
s + p1 s + p2 s + pn
we can derive an especially elegant state-space form
      
ẋ1 −p1 0 x1 1
 ẋ2  
   −p 2 


 x2  
  1 

 ..  =  ..   .. + .. u (2.48)
 .   .   .   . 
ẋn 0 −pn xn 1
 
x1
 x2

 

y= c1 c2 · · · cn  ..  + b0 u (2.49)
 . 
xn
that is purely diagonal. As evident from the diagonal system matrix, this system is totally decou-
pled and possess the best numerical properties for simulation. For systems with repeated roots,
or complex roots, or both, the closest we can come to a diagonal form is the block Jordan form.

We can interconvert between all these above forms using the M ATLAB canon command.
2.7. MULTIVARIABLE CONTROL AND STATE SPACE ANALYSIS 49

>> G = tf([5 6],[1 2 3 4]); % Define transfer function G = (5s2 + 6s)/(s3 + 2s2 + 3s + 4)
2

>> canon(ss(G),'compan') % Convert TF to the companion or observable canonical form

a =
x1 x2 x3
7 x1 0 0 -4
x2 1 0 -3
x3 0 1 -2

b =
12 u1
x1 1
x2 0
x3 0

17 c =
x1 x2 x3
y1 0 5 -4

d =
22 u1
y1 0

These transformations are discussed further in [150, p515]. M ATLAB uses a variation of this form
when converting from a transfer function to state-space in the tf2ss function.

2.7.3 I NTERCONVERTING BETWEEN STATE SPACE AND TRANSFER FUNCTIONS

The transfer function form and the state-space form are, by in large, equivalent, and we can con-
vert from one representation to the other. Starting with our generic state-space model, Eqn. 2.37,

ẋ = Ax + Bu (2.50)

with initial condition x(0) = 0 and taking Laplace transforms gives

sx(s) = Ax(s) + Bu(s)

x(s) (sI − A) = Bu(s)
x(s) = (sI − A)−1 B u(s) (2.51)
| {z }
Gp (s)

where Gp (s) is a matrix of expressions in s and is referred to as the transfer function matrix. The
M ATLAB command ss2tf (state-space to transfer function) performs this conversion.

State-space to transfer function conversion

Suppose we wish to convert the following state-space description to transfer function form.
   
−7 10 4 −1
ẋ = x+ u (2.52)
−3 4 2 −3
50 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

From Eqn. 2.51, the transfer function matrix is defined as

Gp (s) = (sI − A)−1 B

    −1  
1 0 −7 10 4 −1
= s −
0 1 −3 4 2 −3
 −1  
(s + 7) −10 4 −1
= (2.53)
3 (s − 4) 2 −3

Now at this point, the matrix inversion in Eqn. 2.53 becomes involved because of the presence of
the symbolic s in the matrix. Inverting the symbolic matrix expression gives the transfer function
matrix.
  
1 s − 4 −10 4 −1
Gp (s) =
(s + 2)(s + 1) −3 s + 7 2 −3
 
4 s + 26
− 2
 s+2 s + 3s + 2 
= 2

3(6 + s) 
− 2
s+2 s + 3s + 2
We can directly apply Eqn. 2.51 using the symbolic toolbox in M ATLAB.

>> syms s
2 >> A = [-7 10; -3 4]; B = [4 -1; 2 -3];
>> G = (s*eye(2)-A)\B % Pulse transfer function Gp (s) = (sI − A)−1 B.
G =
[ 4/(s+2), -(s+26)/(2+sˆ2+3*s)]
[ 2/(s+2), -3*(s+6)/(2+sˆ2+3*s)]

The method for converting from a state space description to a transfer function matrix described
by Eqn. 2.51 is not very suitable for numerical computation owing to the symbolic matrix in-
version required. However [176, p35] describes a method due to Faddeeva that is suitable for
numerical computation in M ATLAB.

Faddeeva’s algorithm to convert a state-space description, ẋ = Ax + Bu into transfer function

form, Gp (s).

1. Compute a, the characteristic equation of the n × n matrix A. (Use poly in M ATLAB.)

2. Set En−1 = I.
3. Compute recursively the following n − 1 matrices

En−1−k = AEn−k + an−k I, k = 1, 2, · · · , n − 1

4. The transfer function matrix is then given by

sn−1 En−1 + sn−2 En−2 + · · · + E0

Gp (s) = B (2.54)
an sn + an−1 sn−1 + · · · + a1 s + a0

We can repeat the state-space to transfer function example, Eqn. 2.52,

     
ẋ1 −7 10 4 −1
= x+ u
ẋ2 −3 4 2 −3
2.7. MULTIVARIABLE CONTROL AND STATE SPACE ANALYSIS 51

given on page 49 using Faddeeva’s algorithm. We start by computing the characteristic polyno-
mial of A
a(s) = s2 + 3s + 2 = (s + 1)(s + 2)
and with n = 2 we can compute
   
1 0 5 −1
E1 = , E0 =
0 1 2 2

Following Eqn. 2.54 gives a matrix of polynomials in s,

    
−1 1 1 0 −4 10
(sI − A) = 2 s +
s + 3s + 2 0 1 −3 7
 
1 s−4 10
=
(s + 1)(s + 2) −3 s+7

which we can recognise as the same expression as we found using the symbolic matrix inverse
on page 50. All that remains is to post multiply by B to obtain the pulse transfer function matrix.

We can implement Faddeeva’s algorithm using the symbolic toolbox as

syms s
a = poly(A); n = size(A,1);
as = poly2sym(a,'s'); % Construct characteristic polynomial, a(s).
4 E = eye(n,n); I = eye(n,n);
Gp = sˆ(n-1)*E; % En−1 = I
for i=n-2:-1:0
E = A*E + a(i+2)*I;
Gp = Gp + sˆi*E;
9 end
Gp = simplify(Gp*B/as); % Eqn. 2.54.

Finally, we can also use the control toolbox for the conversion from state-space to transfer func-
tion form. First we construct the state-space form,

A = [-7 10; -3,4]; B = [4 -1; 2 -3]; % Continuous state-space example from Eqn. 2.52.
C = eye(size(A));
sys = ss(A,B,C,[]) % construct system object

where in this case if leave out the D matrix, M ATLAB assumes no direct feed though path. We
now can convert to the transfer function form:

>> sys_tf = minreal(tf(sys)) % Convert to transfer function & remember to cancel common
factors
2

Transfer function from input 1 to output...

4
#1: -----
s + 2
7

2
#2: -----
s + 2

12 Transfer function from input 2 to output...

52 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

-s - 26
#1: -------------
sˆ2 + 3 s + 2

17 -3 s - 18
#2: -------------
sˆ2 + 3 s + 2

Once again we get the same pulse transfer function matrix. Note however that it is also possible
to use ss2tf.

T RANSFER FUNCTION TO STATE SPACE

To convert in the reverse direction, we can use the M ATLAB function tf2ss (transfer function to
state space) thus converting an arbitrary transfer function description to the state space format of
Eqn 2.37. For example starting with the transfer function

(s + 3)(s + 4) s2 + 7s + 12
G(s) = = 3
(s + 1)(s + 2)(s + 5) s + 8s + 17s + 10
we can convert to state-space form using tf2ss.

1 >>num = [1 7 12]; % Numerator: s2 + 7s + 12

>>den = [1 8 17 10]; % Denominator: s3 + 8s + 17s + 10
>>[A,B,C,D] = tf2ss(num,den) % convert to state-space

which returns the following four matrices

  
−8 −17 −10 1
A= 1 0 0 , B= 0 
0 1 0 0
 
C = 1 7 12 , D=0

This form is a variation of the controllable canonical form described previously in §2.7.2. This
is not the only state-space realisation possible however as the C ONTROL T OOLBOX will return a
slightly different ABCD package

>> Gc = tf([1 7 12],[1 8 17 10]) % Transfer function in polynomial form

2 Transfer function:
sˆ2 + 7 s + 12
-----------------------
sˆ3 + 8 sˆ2 + 17 s + 10

7 >>[A,B,C,D] = ssdata(Gc) % extract state-space matrices

A =
-8.0000 -4.2500 -1.2500
4.0000 0 0
0 2.0000 0
12 B =
2
0
0
C =
2.7. MULTIVARIABLE CONTROL AND STATE SPACE ANALYSIS 53

17 0.5000 0.8750 0.7500

D =
0

Alternatively we could start with zero-pole-gain form and obtain yet another equivalent state-
space form.

1 G = zpk([-3 -4],[-1 -2 -5],1) % Transfer function model in factored format

Gss = ss(G)

A later section, (§2.7.4) shows how we can convert between equivalent dynamic forms. The four
matrices in the above description form the linear dynamic system as given in Eqn. 2.41. We can
concatenate these four matrices into one large block thus obtaining a shorthand way of storing
these equations

# of inputs
n m
✛ ✲✛ ✲
✻
A B n # of states
G= (n × n) (n × m)
❄
C D ✻p # of measurements
(p × n) (p × m) ❄

which is often called a packed matrix notation of the quadruplet (A,B,C,D). If you get confused
about the appropriate dimensions of A, B etc, then this is a handy check, or alternatively you
can use the diagnostic routine abcdchk. Note that in practice, the D matrix is normally the
zero matrix since the input does not generally effect the output immediately without first passing
through the system dynamics.

2.7.4 S IMILARITY TRANSFORMATIONS

The state space description of a differential equation is not unique. We can transform from one
description to another by using a linear invertible transformation such as x = Tz. Geometrically
in 2 dimensions, this is equivalent to rotating the axes or plane. When one rotates the axes, the
inter-relationship between the states do not change, so the transformation preserves the dynamic
model.

Suppose we have the dynamic system

ẋ = Ax + Bu (2.55)

which we wish to transform in some manner using the non-singular transformation matrix, T,
where x = Tz. Naturally the reverse transformation z = T−1 x exists because we have restricted
ourselves to consider only the cases when the inverse of T exists. Writing Eqn. 2.55 in terms of
54 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

our new variable z we get

d(Tz)
= ATz + Bu
dt
ż = T−1 ATz + T−1 Bu (2.56)
Eqn 2.56 and Eqn 2.55 represent the same dynamic system. They have the same eigenvalues
hence the similarity transform, but just a different viewpoint. The mapping from A to T−1 AT
is called a similarity transform and preserves the eigenvalues. These two matrices are said to be
similar. The proof of this is detailed in [87, p300] and [150, p513–514].

The usefulness of these types of transformations is that the dynamics of the states are preserved
(since the eigenvalues are the same), but the shape and structure of the system has changed.
The motivation is that for certain operations (control, estimation, modelling), different shapes
are more convenient. A pure, (or nearly), diagonal shape of the A matrix for example has much
better numerical properties than a full matrix. This also has the advantage that less computer
storage and fewer operations are needed.

To convert a system to a diagonal form, we use the transformation matrix T, where the columns
of T are the eigenvectors of A. Systems where the model (A) matrix is diagonal are especially
easy to manipulate. We can find the eigenvectors T and the eigenvalues e of a matrix A using
the eig command, and construct the new transformed, and hopefully diagonal system matrix,

[T,e] = eig(A) % find eigenvectors & eigenvalues of A

V = T\A*T % New system matrix, T−1 AT, or use canon

Another use is when testing new control or estimation algorithms, it is sometimes instructive to
devise non-trivial systems with specified properties. For example you may wish to use as an
example a 3 × 3 system that is stable and interacting, and has one over-damped mode and two
oscillatory modes. That is we wish to construct a full A matrix with specified eigenvalues. We
can use the similarity transformations to obtain these models.

Other useful transformations such as the controllable and observable canonical forms are covered
in [150, §6–4 p646 ]. The M ATLAB function canon can convert state-space models to diagonal
or observable canonical form (sometimes known as the companion form). Note however the
help file for this routine discourages the use of the companion form due to its numerical ill-
conditioning.

M ATLAB comes with some utility functions to generate test models such as ord2 which gener-
ates stable second order models of a specified damping and natural frequency, and rmodel is a
flexible general stable random model builder or arbitrary order.

2.7.5 I NTERCONVERTING BETWEEN TRANSFER FUNCTIONS FORMS

The previous section described how to convert between different representations of linear dy-
namic systems such as differential equations, transfer functions and state-space descriptions.
This section describes the much simpler task of converting between the different ways we can
write transfer functions.

Modellers tend to think in continuous time systems, G(s), and in terms of process gain and time
constants, so will naturally construct transfer functions in factored form as
K Πi (αi s + 1)
(2.57)
Πj (τj s + 1)
2.7. MULTIVARIABLE CONTROL AND STATE SPACE ANALYSIS 55

where the all the variables of interest such as time constants τj are immediately apparent. On
the other hand, system engineers tend to think in terms of poles and zeros, so naturally construct
transfer functions also in factored form, but with a different scaling

K ′ Πi (s + zi )
(2.58)
Πj (s + pi )

so where now the poles, pj , and zeros zi are immediately apparent. This is the form that M ATLAB
uses in the zeros-pole-gain format, zpk.

Finally the hardware engineer would prefer to operate in the expanded polynomial form, (par-
ticularly in discrete cases), where the transfer function is of the form

bm sm + bm−1 sm−1 + · · · + b0
(2.59)
sn + an−1 sn−1 + · · · + a0

This is the form that M ATLAB uses in the transfer-function format, tf. Note that the leading
coefficient in the denominator is set to 1. As a M ATLAB user, you can define a transfer function
that does not follow this convention, but M ATLAB will quietly convert to this normalised form if
you type something like G = tf(zpk(G)).

The inter-conversions between the forms is not difficult; between expressions Eqn. 2.57 and
Eqn. 2.58 simply require some adjusting of the gains and factors, while to convert from Eqn. 2.59
requires one to factorise polynomials.

For example, the following three transfer function descriptions are all equivalent

2 (10s + 1)(−3s + 1) −1.5 (s − 0.3333)(s + 0.1) −60s2 + 14s + 2

= =
(20s + 1)(2s + 1)(s + 1) (s + 1)(s + 0.5)(s + 0.05) 40s3 + 62s2 + 23s + 1
| {z } | {z } | {z }
time constant zero-pole-gain expanded polynomial

It is trivial to interconvert between the zero-pole-gain, Eqn. 2.58, and the transfer function for-
mats, Eqn. 2.59, in M ATLAB, but it is less easy to convert to the time constant description. List-
ing 2.3 extracts from an arbitrary transfer function form the time constants, τ , the numerator time
constants, α, and the plant gain K.

Listing 2.3: Extracting the gain, time constants and numerator time constants from an arbitrary
transfer function format
Gplant = tf(2*mconv([10 1],[-3 1]),mconv([20 1],[2 1],[1 1])) % TF of interest
G = zpk(Gplant); % Convert TF description to zero-pole-gain
3

Kp = G.k;
p = cell2mat(G.p); z = cell2mat(G.z); % Extract poles & zeros
delay0 = G.iodelay ; % Extract deadtime (if any)

8 tau = sort(-1./p,'descend'); % Convert (& sort) to time constants, 1/(τj s + 1)

ntau = sort(-1./z,'descend'); % Convert to numerator time constants, (αi s + 1)
K = Kp*prod(tau)/prod(ntau); % Adjust plant gain

We could of course use the control toolbox functions pole and zero to extract the poles and
zeros from an arbitrary LTI model.
56 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

2.7.6 T HE STEADY STATE

The steady-state, xss , of a general nonlinear system ẋ = f (x, u) is the point in state space such
that all the derivatives are zero, or the solution of

0 = f (xss , u) (2.60)

If the system is linear, ẋ = Ax+ Bu, then the steady state can be evaluated algebraically in closed
form
xss = −A−1 Bu (2.61)

Consequently to solve for the steady-state one could invert the model matrix A but which may be
illconditioned or computationally time consuming. Using M ATLAB one should use xss = -A\B*u.

If A has no inverse, then no (or alternatively infinite) steady states exist. An example of a process
that has no steady state could be a tank-flow system that has a pump withdrawing fluid from
the outlet at a constant rate independent of liquid height say just exactly balancing an input flow
shown in the left hand schematic in Fig. 2.21. If the input flow suddenly increased, then the level
will rise until the tank eventually overflows. If instead the tank was drained by a valve partially
open at a constant value, then as the level rises, the increased pressure (head) will force more
material out through the valve, (right-hand side of Fig. 2.21.) Eventually the system will rise to a
new steady state. It may however overflow before the new steady state is reached, but that is a
constraint on the physical system that is outside the scope of the simple mathematical description
used at this time.

flow in

❄
❄
✻
✻
h h
❄ ❄

✲ flow = constant ✲ flow = f(height)

constant flow pump restriction

Figure 2.21: Unsteady and steady states for level systems

If the system is nonlinear, there is the possibility that multiple steady states may exist. To solve
for the steady state of a nonlinear system, one must use a nonlinear algebraic solver such as
described in chapter 3.

Example The steady state of the differential equation

d2 y dy
2
+7 + 12y = 3u
dt dt
2.8. SOLVING THE VECTOR DIFFERENTIAL EQUATION 57

where dydt = y = 0 at t = 0 and u = 1 , t ≥ 0 can be evaluated using Laplace transforms and the
final value theorem. Transforming to Laplace transforms we get

s2 Y (s) + 7sY (s) + 12Y (s) = 3U (s)

Y (s) 3 3
= 2 =
U (s) s + 7s + 12 (s + 3)(s + 4)

while for a step input

1 3
Y (s) = ·
s s2 + 7s + 12
The final value theorem is only applicable if the system is stable. To check, we require that the
roots of the denominator, s2 + 7s + 12, lie in the left hand plane,
√
−6 ± 72 − 4 × 12
s= = −4 and − 3
2
Given that both roots have negative real parts, we have verified that our system is stable and we
are allowed to apply the final value theorem to solve for the steady-state, y(∞),
3 3
lim y(t) = lim sY (s) = lim = = 0.25
t→∞ s→0 s→0 s2 + 7s + 12 12
Using the state space approach to replicate the above, we first cast the second order differential
equation into two first order differential equations using the controllable canonical form given in
Eqn. 2.45. Let z1 = y and z2 = dydt , then
   
0 1 0
ż = z+ u
−12 −7 3

and the steady state is now

    
−1 −7/12 −1/12 0 0.25
zss = −A Bu = − 1=
1 0 3 0

Noting that z1 = y, we see that the steady state is also at 0.25. Furthermore, the derivative term
(z2 ) is zero, which is as expected at steady state.

2.8 S OLVING THE VECTOR DIFFERENTIAL EQUATION

Since we can solve differential equations by inverting Laplace transforms, we would expect to be
able to solve state-space differential equations such as Eqn. 2.37 in a similar manner. If we look
at the Laplace transform of a simple linear scalar differential equation, ẋ = ax + bu we find two
terms,

sX(s) − x0 = aX(s) + bU (s)

x0 b
X(s) = + U (s) (2.62)
s−a s−a
One of these terms is the response of the system owing to the initial condition, x0 eat , and is
called the homogeneous solution, and the other term is due to the particular input we happen
to be using. This is called the particular integral, and we must know the form of the input, u(t),
before we solve this part of the problem. The total solution is the sum of the homogeneous and
particular components.
58 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

T HE HOMOGENEOUS SOLUTION

For the moment we will consider just the homogeneous solution to our vector differential equa-
tion. That is, we will assume no driving input, or u(t) = 0. (In the following section we will add
in the particular integral due to a non-zero input.)

Our vector differential equation, ignoring any input, is simply

ẋ = Ax, x(t = 0) = x0 (2.63)

Taking Laplace transforms, and not forgetting the initial conditions, we have

sx(s) − x0 = A x(s)
x(s) (sI − A) = x0 (s)
−1
x(s) = (sI − A) x0

Finally inverting the Laplace transform back to the time domain gives
n o
−1
x(t) = L−1 (sI − A) x0 (2.64)

Alternatively we can solve Eqn. 2.63 by separating the variables and integrating obtaining

x(t) = eAt x0 (2.65)

where the exponent of a matrix, eAt , is itself a matrix of the same size as A. We call this matrix
exponential the transition matrix because it transforms the state vector at some initial time x0 to
some point in the future, xt . We will give it the symbol Φ(t).

The matrix exponential is defined just as in the scalar case as a Taylor series expansion,

A2 t2 A3 t3
eAt or Φ(t) = I + At +
def
+ + ···
2! 3!
∞
X Ak tk
= (2.66)
k!
k=0

although this series expansion method is not recommended as a reliable computational strategy.
Better strategies are outlined on page 62.
n o
Comparing Eqn. 2.64 with Eqn. 2.65 we can see that the matrix eAt = L−1 (sI − A)
−1
.

So to compute the solution, x(t), we need to know the initial condition and a strategy to numeri-
cally compute a matrix exponential.

T HE PARTICULAR SOLUTION

Now we consider the full differential equation with nonzero input ẋ = Ax + Bu. Building on
the solution to the homogeneous part in Eqn. 2.65, we get
Z t
x(t) = Φt x0 + Φt−τ Buτ dτ (2.67)
0

where now the second term accounts for the particular input vector u.
2.8. SOLVING THE VECTOR DIFFERENTIAL EQUATION 59

Eqn. 2.67 is not particularly useful as written as both terms are time varying. However the con-
tinuous time differential equation can be converted to a discrete time differential equation that
is suitable for computer control implementation provided the sampling rate is fixed, and the in-
put is held constant between the sampling intervals. We would like to convert Eqn. 2.37 to the
discrete time equivalent, Eqn. 2.38, repeated here

xk+1 = Φxk + ∆uk (2.68)

where xk is the state vector x at time t = kT where T is the sample period. Once the sample pe-
riod is fixed, then Φ and ∆ are also constant matrices. We have also assumed here that the input
vector u is constant, or held, over the sample interval, which is the norm for control applications.

So starting with a known xk at time t = kT , we desire the state vector at the next sample time,
xk+1 , or
Z (k+1)T
A T
xk+1 = e xk + e A (k+1)T
e−Aτ Buτ dt (2.69)
kT

But as we have assumed that the input u is constant using a zeroth-order hold between the
sampling intervals kT and (k + 1)T , Eqn. 2.69 simplifies to
Z T
xk+1 = eAT xk + eAλ Buk dλ (2.70)
0

where λ = T − t. For convenience, we can define two new matrices as

Φ = eAT
def
(2.71)
Z T !
∆=
def
e Aλ
dλ B (2.72)
0

which gives us our desired transformation in the form of Eqn. 2.68. In summary, to discretise
ẋ = Ax + Bu at sample interval T , we must compute a matrix exponential, Eqn. 2.71, and
integrate a matrix exponential, Eqn. 2.72.

Note that Eqn. 2.70 involves no approximation to the continuous differential equation provided
the input is constant over the sampling interval. Also note that as the sample time tends to zero,
the state transition matrix Φ tends to the identity, I.

If the matrix A is non-singular, then

 
∆ = eAT − I A−1 B (2.73)

is a simpler expression than the more general Eqn. 2.72.

Double integrator example. The single-input/single output double integrator system, G(s) =
1/s2 , can be represented in continuous state space using two states as
   
0 0 1
ẋ = x+ u (2.74)
1 0 0
 
y= 0 1 x (2.75)

At sample time T , the state transition matrix is from Eqn. 2.71,

   
0 0 1 0
Φ = eAT = exp =
T 0 T 1
60 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

and the control matrix is given by Eqn. 2.72

Z !
T
∆= e Aλ
dλ B
0
Z T   ! 
1 0 1
= dλ
0 λ 1 0
   
T 0  1  T
=  T2  =  T2 
T 0
2 2
For small problems such as this, the symbolic toolbox helps do the computation.

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

>> syms T lambda % Symbolic sample time T and λ
>> Phi = expm(A*T) % Φ(T )
Phi =
5 [ 1, 0]
[ T, 1]
>> Delta = int(expm(A*lambda),lambda,0,T)*B % ∆(T )
Delta =
[ T]
10 [ 1/2*Tˆ2]

Symbolic example with invertible A matrix. We can discretise the continuous state-space sys-
tem    
−1.5 −0.5 1
ẋ = x+ u
1 0 0
analytically at a sample time T by computing matrix exponentials symbolically.

>> A = [-3/2, -1/2; 1, 0]; B = [1;0];

>> syms T
>> Phi = expm(A*T)
5 Phi =
[ -exp(-1/2*T)+2*exp(-T), exp(-T)-exp(-1/2*T)]
[ -2*exp(-T)+2*exp(-1/2*T), 2*exp(-1/2*T)-exp(-T)]

Since A is invertible in this example, we can use the simpler Eqn. 2.73

 
>> Delta = (Phi - eye(2))*A\B % ∆ = eAT − I A−1 B
Delta =
3 [ 2/(exp(-1/2*T)*exp(-T)-exp(-T)-exp(-1/2*T)+1)*(exp(-T)-exp(-1/2*T))]
[ -2*(-exp(-1/2*T)+2*exp(-T)-1)/(exp(-1/2*T)*exp(-T)-exp(-T)-exp(-1/2*T)+1)]
>> simplify(Delta)
ans =
[ 2*exp(-1/2*T)/(exp(-T)-1)]
8 [ -2*(2*exp(-1/2*T)+1)/(exp(-T)-1)]

Of course it is evident from the above example that the symbolic expressions for Φ(T ) and ∆(T )
rapidly become unwieldy for dimensions much larger than about 2. For this reason, analytical
2.8. SOLVING THE VECTOR DIFFERENTIAL EQUATION 61

expressions are of limited practical worth. The alternative numerical schemes are discussed in
the following section.

2.8.1 N UMERICALLY COMPUTING THE DISCRETE TRANSFORMATION

Calculating numerical values for the matrices Φ and ∆ can be done by hand for small dimensions
by converting to a diagonal or Jordan form, or numerically using the exponential of a matrix.
Manual calculations are neither advisable nor enjoyable, but [19, p35] mention that if you first
compute
Z T
AT 2 A2 T 3 Ak−1 T k
Ψ= eAτ dτ = IT + + + ···+ + ···
0 2! 3! k!
then
Φ = I + AΨ and ∆ = ΨB (2.76)
A better approach, at least when using M ATLAB, follows from Eqn. 2.71 and Eqn. 2.72 where

dΦ
= ΦA = AΦ (2.77)
dt
d∆
= ΦB (2.78)
dt
These two equations can be concatenated to give
    
d Φ ∆ Φ ∆ A B
= (2.79)
dt 0 I 0 I 0 0
R R
which is in the same form as da/dt = ab. Rearranging this to da/a = b dt leads to the analytical
solution     
Φ ∆ A B
= exp T (2.80)
0 I 0 0
enabling us to extract the required Φ and ∆ matrices provided we can reliably compute the
exponential of a matrix. The M ATLAB C ONTROL T OOLBOX routine to convert from continuous
to discrete systems, c2d, essentially follows this algorithm.

We could try the augmented version of Eqn. 2.80 to compute both Φ and ∆ with one call to the
matrix exponential function for the example started on page 60.

>> A = [-3/2, -1/2; 1, 0]; B = [1;0];

2 >> syms T

>> [n,m] = size(B); % extract dimensions  

A B
>> Aa = [A,B;zeros(m,n+m)] % Augmented A matrix, Ã =
0 0
Aa =
7 -1.5 -0.5 1.0
1.0 0.0 0.0
0.0 0.0 0.0
>> Phi_a = expm(Aa*T) compute exponential, Φ̃ = exp(ÃT )
Phi_a =
12 [ -exp(-1/2*T)+2*exp(-T), exp(-T)-exp(-1/2*T), -2*exp(-T)+2*exp(-1/2*T)]
[-2*exp(-T)+2*exp(-1/2*T), 2*exp(-1/2*T)-exp(-T),2*exp(-T)-4*exp(-1/2*T)+2]
[ 0, 0, 1]
>> Phi = Phi_a(1:n,1:n)
Phi =
62 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

17 [ -exp(-1/2*T)+2*exp(-T), exp(-T)-exp(-1/2*T)]
[ -2*exp(-T)+2*exp(-1/2*T), 2*exp(-1/2*T)-exp(-T)]
>> Delta = Phi_a(1:n,n+1:end)
Delta =
[ -2*exp(-T)+2*exp(-1/2*T)]
22 [ 2*exp(-T)-4*exp(-1/2*T)+2]

R ELIABLY COMPUTING MATRIX EXPONENTIALS

There are in fact several ways to numerically compute a matrix exponential. Ogata gives three
computational techniques, [150, pp526–533], and a paper by Cleve Moler (one of the original
M ATLAB authors) and Charles Van Loan is titled Nineteen dubious ways to compute the exponential
of a matrix5 , and contrasts methods involving approximation theory, differential equations, matrix
eigenvalues and others. Of these 19 methods, two found their way into M ATLAB, namely expm,
which is the recommended default strategy, and expm1 intended to compute ex − 1 for small x .

The one time when matrix exponentials are trivial to compute is when the matrix is diagonal.
Physically this implies that the system is totally decoupled since the matrix A is comprised of
only diagonal elements and the corresponding exponential matrix is simply the exponent of the
individual elements. So given the diagonal matrix
 
λ1 0 0
D =  0 λ2 0 
0 0 λ3

then the matrix exponential is simply

 
eλ1 0 0
D
e =  0 eλ2 0 
0 0 eλ3

which is trivial and reliable to compute. So one strategy then is to transform our system to a
diagonal form (if possible), and then simply find the standard scalar exponential of the individual
elements. However some matrices, such as those with multiple eigenvalues, is is impossible to
convert to diagonal form, so in those cases the best we can do is convert the matrix to a Jordan
block form as described in [150, p527], perhaps using the jordan command from the Symbolic
toolbox.

However this transformation is very sensitive to numerical roundoff and for that reason is not
used for serious computation. For example the matrix
 
1 1
A=
ǫ 1

with ǫ = 0 has a Jordan form of  

1 1
0 1
but for ǫ 6= 0, then the Jordan form drastically changes to the diagonal matrix
 √ 
1+ ǫ 0√
0 1− ǫ
5 SIAM Review, vol. 20, A 1978, pp802–836 and reprinted in [159, pp649–680]. In fact 25 years later, an update was

published with some recent developments.

2.8. SOLVING THE VECTOR DIFFERENTIAL EQUATION 63

In summary, for serious numerical calculation we should use the matrix exponential function
expm. Remember not to confuse finding the exponent of a matrix, expm, with the M ATLAB func-
tion exp which simply finds the exponent of all individual elements in the matrix.

2.8.2 U SING M ATLAB TO DISCRETISE SYSTEMS

All of the complications of discretising continuous systems to their discrete equivalent can be
circumvented by using the M ATLAB command c2d which is short for continuous to discrete.
Here we need only to pass the continuous system of interest, and the sampling time. As an
example we can verify the conversion of the double integrator system shown on page 2.8.

G = tf(1,[1 0 0]) % Continuous system G(s) = 1/s2 in transfer function form

Gc = ss(G) % Convert to continuous state-space
3 Gd = c2d(Gc,2) % Convert to discrete state-space with a sample time of T = 2

a =
x1 x2
x1 1 0
8 x2 2 1

b =
u1
x1 2
13 x2 2

c =
x1 x2
y1 0 1
18

d =
u1
y1 0

23 Sampling time (seconds): 2

Discrete-time state-space model.

Unlike in the analytical case presented previously, here we must specify a numerical value for
the sample time, T .

Example: Discretising an underdamped second order system with a sample time of T = 3

following Eqn. 2.80 using M ATLAB to compute the matrix exponential of the augmented system.

1 [A,B,C,D] = ord2(3, 0.5); % generate a second order system

Ts = 3.0; % sample time

[na,nb] = size(B);
X = expm([A,B; zeros(nb,na+nb)]*Ts); % matrix exponential
6 Phi = X(1:na,1:na); % Pull off blocks to obtain Φ & ∆.
Del = X(1:na,na+1:end);

We can use M ATLAB to numerically verify the expression found for the discretised double inte-
grator on page 59 at a specific sample time, say T = 4.
64 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

The double integrator in input/output

form is
1 >>Gc = tf(1,[1,0,0])% G(s) = 1/s2
G(s) = 2 >>Gc_ss = ss(Gc);% state space
s
3 >>dt = 4; % sample time
which in packed state space notation is
  >>Gd_ss = c2d(Gc_ss,dt)
  0 0 1 a =
A B
= 1 0 0  x1 x2
C D
0 1 0 8 x1 1.000 0
x2 4.000 1.000
Using a sample time of T = 4, we can com-
pute b =
    u1
1 0 1 0 13 x1 4.000
Φ= =
T 1 4 1 x2 8.000
   
T 4
∆ = T 2 = 16 c =
2 2 x1 x2
18 y1 0 1.000
We can verify this using c2d in M ATLAB as
demonstrated opposite. d =
u1
y1 0

Once the two system matrices Φ and ∆ are known, solving the difference equation for further
values of k just requires the repeated application of Eqn. 2.38. Starting from known initial condi-
tions x0 , the state vector x at each sample instant is calculated thus

x1 = Φx0 + ∆u0
x2 = Φx1 + ∆u1 = Φ2 x0 + Φ∆u0 + ∆u1
x3 = Φx2 + ∆u2 = Φ3 x0 + Φ2 ∆u0 + Φ∆u1 + ∆u2
.. ..
. .
n−1
X
xn = Φn x0 + Φn−1−k ∆uk (2.81)
k=0

The M ATLAB function ltitr which stands for linear time invariant time response or the more
general dlsim will solve the general linear discrete time model as in 2.81. In special cases such
as step or impulse tests, you can use dstep, or dimpulse.

Problem 2.4 1. Evaluate An where  

1 1
A=
1 0
for different values of n. What is so special about the elements of An ? (Hint: find the
eigenvalues of A.)
2. What is the determinant of An ?
3. Write a couple of M ATLAB m-functions to convert between the A, B, C and D form and the
packed matrix form.
4. Complete the state-space to transfer function conversion analytically started in §2.7.3,
Eqn. 2.53. Compare your answer with using ss2tf.
2.8. SOLVING THE VECTOR DIFFERENTIAL EQUATION 65

A submarine example A third-order model of a submarine from [63, p416 Ex9.17] is

   
0 1 0 0
ẋ =  −0.0071 −0.111 0.12  x +  −0.095  u (2.82)
0 0.07 −0.3 0.072
where the state vector x is defined as
 
x⊤ = θ dθ
dt α
The state θ is the inclination of the submarine and α is the angle of attack above the horizontal.
The scalar manipulated variable u is the deflection of the stern plane. We will assume that of the
three states, we can only actually measure two of them, θ and α, thus
 
1 0 0
C=
0 0 1
In this example, we consider just three states x, one manipulated variable u, and two outputs, y.

The stability of the open-loop system is given by the eigenvalues of A. In M ATLAB, the command
eig(A) returns  
−0.0383 + 0.07i
eig(A) =  −0.0383 − 0.07i 
−0.3343
showing that all eigenvalues have negative real parts, indicating that the submarine is oscillatory,
though stable. In addition, the complex conjugate pair indicates that there will be some oscillation
to the response of a step change in stern plane.

We can simulate the response of our submarine system (Eqn 2.82) to a step change in stern plane
movement, u, using the step command. The smooth plot in Fig. 2.22 shows the result of this
continuous simulation, while the ‘staired’ plot shows the result of the discrete simulation using
a sample time of T = 5. We see two curves corresponding to the two outputs.

Listing 2.4: Submarine simulation

A = [0,1,0;-0.0071,-0.111,0.12;0,0.07,-0.3]; B = [0,-0.095,0.072]';
C = [1 0 0;0 0 1];
3

Gc = ss(A,B,C,0) % Continuous plant, Gc (s)

dt = 5; % sample time T = 5 (rather coarse)

Gd = c2d(Gc,dt) % Create discrete model
8

step(Gc,Gd); % Do step response and see Fig. 2.22.

Fig. 2.22 affirms that the openloop process response is stable, supporting the eigenvalue analysis.
Notice how the step command automatically selected appropriate time scales for the simulation.
How did it do this?

The steady state of the system given that u = 1 is, using Eqn 2.61,

1 >> u = 1;
>>xss = -A\B*u; % Steady-states given by xss = −ABu.
>>yss = C*xss
yss = % See steady-state outputs in Fig. 2.22.
-9.3239
6 0.2400
66 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

0
G
c
G
−5 d

1
x
−10

−15

0.4

2 0.3

0.2
x

0.1
Figure 2.22: Comparing the continuous
0
and discrete step responses of the subma- 0 20 40 60 80
rine model. time

which corresponds to the final values given in Fig. 2.22.

We can verify the results from the c2d routine using Eqns 2.71 and 2.72, or in this case since A is
invertible, Eqn. 2.73.

>> Phi = expm(A*5.0); % Don’t forget the ‘m’ in expmatrix()

 
>> Delta = (Phi-eye(size(A)))*(A\B) % ∆ = eAT − I A−1 B

This script should give the same Φ and ∆ matrices as the c2d function since A is non-singular.

Given a discrete or continuous model, we can compute the response to an arbitrary input se-
quence using lsim:

t = [0:dt:100]'; % time vector

U = sin(t/20); % Arbitrary input U (t)
3

lsim(Gc,U,t) % continuous system simulation

lsim(Gd,U) % discrete system simulation

We could explicitly demand a first-order-hold as opposed to the default zeroth order by setting
the options for c2d.

optn = c2dOptions('Method','foh') % Ask for a first-order hold

Gd2 = c2d(Gc,5,optn) % Compare the step response of this with a zeroth-order hold

Problem 2.5 By using a suitable state transformation, convert the following openloop system

ẋ1 = x2
ẋ2 = 10 − 2x2 + u
2.8. SOLVING THE VECTOR DIFFERENTIAL EQUATION 67

where the input u is the following function of the reference r,

 
u = 10 + 9r − 9 4 x

into a closed loop form suitable for simulation in M ATLAB using the lsim function. (ie: ẋ =
Ax + Br) Write down the relevant M ATLAB code segment. (Note: In practice, it is advisable to
use lsim whenever possible for speed reasons.)

Problem 2.6 1. Write an m-file that returns the state space system in packed matrix form that
connects the two dynamic systems G1 and G2 together,

Gall = G1 · G2

Assume that G1 and G2 are already given in packed matrix form.

2. (Problem from [185, p126]) Show that the packed matrix form of the inverse of Eqn. 2.41,
u = G−1 y, is  
A − BD−1 C −BD−1
G−1 = (2.83)
D−1 C D−1
assuming D is invertible, and we have the same number of inputs as outputs.

2.8.3 T HE DISCRETE STATE SPACE EQUATION WITH TIME DELAY

The discrete state space equation, Eqn. 2.37, does not account for any time delay or dead time
between the input u and the output x. This makes it difficult to model systems with delay.
However, by introducing new variables we can accommodate dead time or time delays. Consider
the continuous differential equation,

ẋt = Axt + But−θ (2.84)

where θ is the dead time and in this particular case, is exactly equal to 2 sample times (θ = 2T ).
The discrete time equivalent to Eqn 2.84 for a two sample time delay is

xk+3 = Φxk+2 + ∆uk (2.85)

which is not quite in our standard state-space form of Eqn. 2.38 owing to the difference in time
subscripts between u and x.

Now let us introduce a new vector of state variables, z where

 
xk
def
zk =  xk+1  (2.86)
xk+2

is sometimes known as a shift register or tapped delay line of states. Using this new state vector,
we can write the dynamic system, Eqn. 2.85 compactly as
   
xk+1 Ixk+1
zk+1 =  xk+2  =  Ixk+2  (2.87)
xk+3 Φxk+2 + ∆uk

or in standard state-space form as

   
0 I 0 0
zk+1 = 0 0 I  z k +  0  uk (2.88)
0 0 Φ ∆
68 CHAPTER 2. FROM DIFFERENTIAL TO DIFFERENCE EQUATIONS

This augmented system, Eqn. 2.88 (with 2 units of dead time) is now larger than the original
Eqn. 2.85 given that we have 2n extra states. Furthermore, note that the new augmented transi-
tion matrix in Eqn. 2.88 is now no longer of full rank (provided it was in the first place). If we
wish to incorporate a delay θ, of β sample times (θ = βT ), then we must introduce β dummy
state vectors into the system creating a new state vector z with dimensions (nβ × 1);
 ⊤
zk = xk xk+1 xk+2 · · · xk+β

The augmented state transition matrix Ξ and augmented control matrix Ψ are of the form
   
0 I 0 ··· 0 0
 0 0 I ··· 0   0 
   
 0 0 0 ··· 0   0 
   
Ξ= . . . . , and Ψ =  0  (2.89)
 .. .. .. . . ... 
  
   .. 
 0 0 0 ··· I   . 
0 0 0 ··· Φ ∆

Now the discrete time equation is

zk+1 = Ξzk + Ψuk (2.90)
Since all the introduced states are unobservable, we have
 
xk = I 0 0 · · · 0 zk (2.91)

and that the output equation is amended to

 
yk = C I 0 0 ··· 0 zk (2.92)

If the dead time is not an exact multiple of the sample time, then more sophisticated analysis is
required. See [71, p174]. Systems with a large value of dead time become very unwieldy as a
dummy vector is required for each multiple of the sample time. This creates large systems with
many states that possibly become numerically ill-conditioned and difficult to manipulate.

Problem 2.7 1. Simulate the submarine example (Eqn 2.82) but now introducing a dead time
of 2 sample time units.
2. Explain how the functions step and dstep manage to guess appropriate sampling rates,
and simulation horizons. Under what circumstances will the heuristics fail?

2.9 S TABILITY

Stability is a most desirable characteristic of a control system. As a control designer, you should
only design control systems that are stable. For this reason, one must at least be able to analyse
a potential control system in order to see if it is indeed theoretically stable at least. Once the
controller has been implemented, it may be too late and costly to discover that the system was
actually unstable.

One definition of stability is as follows:

A system of differential equations is stable if the output is bounded for all

time for any given bounded input, otherwise it is unstable. This is referred to
as bounded input bounded output (BIBO) stability.
2.9. STABILITY 69

System Stability

✙ s

Transfer function Nonlinear

■ q
✙ Lyapunov
non-poly ✌