The Use of Mathematica in Control Engineering

1
The Use of Mathematica inControl Engineering
Neil Munro Control Systems Centre UMIST Manchester
, England.

• Linear Model Descriptions
• Linear Model Transformations
• Linear System Analysis Tools
• Design/Synthesis Techniques
• Pole Assignment
• Model-Reference Optimal Control
• PID Controller
• Concluding Remarks

2
Linear Model Descriptions

• The Control System Professional currently
  provides
• several ways of describing linear system models
  e.g.
• For systems described by the state-space equations

where y is a vector of the system outputs
u is a vector of the system inputs and x
is a vector of the system state-variables

1. For systems described by transfer-function
   relationships

where s is the complex variable
3
Examples -
ss StateSpace0, 0, 1, 0,0, 0, 0, 1, -(a
b), 0, ab, 0,0, -(a b), 0, ab, 0, 0,0,
0,1, 0,0, 1,-b, -a, 1, 1
4
tf TransferFunctions,1/(s-a),1/(s-b)
5
New Data Formats have been implemented, for these
objects, which are fully editable, as follows -
ss StateSpace0, 0, 1, 0,0, 0, 0, 1, -(a
b), 0, ab, 0,0, -(a b), 0, ab, 0, 0,0,
0,1, 0,0, 1,-b, -a, 1, 1 now results
in the composite data matrix
6
tfsys TransferFunctions, ((s2)(s3))/(s1)2
,1/(s1)2, (s2)/(s1)2,(s1)/((s1)2(s3)),
1/(s2),1/(s1) which now yields
7

• New Data Objects
• Four new data objects have been introduced
  namely,
• Rosenbrocks system matrix in polynomial form
• Rosenbrocks system matrix in state-space form
• The right matrix-fraction description of a system
• The left matrix-fraction description of a system.

8
The system matrix in polynomial form provides a
compact description of a linear dynamical system
described by arbitrary ordered differential
equations and algebraic relationships, after the
application of the Laplace transform with zero
initial conditions namely
where ?, u and y are vectors of the Laplace
transformed system variables. This set of
equations can equally be written as
9
The system matrix in polynomial form is then
defined as
When a system matrix in polynomial form is being
created, it is important to note that the
dimension of the square matrix T(s) must be
adjusted to be r, where r is ? the degree of
DetT(s).
If the system description is known in state-space
form then a special form of the system matrix can
be constructed, known as the system matrix in
state-space form, as shown below
10

• Matrix Fraction Forms
• Given a system description in transfer
  function matrix form G(s), for certain analysis
  and design purposes e.g. the H-? approach to
  robust control system design it is often
  convenient to express this model in the form of a
  left or right matrix-fraction description e.g.
• A left matrix-fraction form of a given transfer
  function matrix G(s) might be
• 2 A right matrix fraction form of a given
  transfer function matrix G(s) might be

11
Linear Model Transformations
G(s)
A, B, C, D
G(s)
System Matrix P(s) in polynomial form
G(s)
System Matrix P(s) in state-space form
A, B, C, D
System Matrix P(s) in polynomial form
T(s),U(s),V(s),W(s)
Least Order
P1(s) ? P2(s)
12
New Data Transformations
All data formats are fully editable
tfsys TransferFunctions, ((s2)(s3))/(s1)2
,1/(s1)2, (s2)/(s1)2,(s1)/((s1)2(s3)),
1/(s2),1/(s1)
ss StateSpacetfsys
rff RightMatrixFractiontfsys
rff RightMatrixFractionss,s
13
New Data Transformations
tfsys TransferFunctions,
(s1)/(s22s1),(s2)/(s1)
ps SystemMatrixtfsys,TargetForm?RightFraction
rff RightMatrixFractionps
TransferFunction
ps SystemMatrixtfsys,TargetForm ?LeftFraction
lf LeftFractionForm

14
New Data Transformations
tfsys TransferFunctions, (s1)/(s22s1),
(s2)/(s1)
rff RightMatrixFractiontfsys
dt ToDiscreteTimetfsys,Sampled-gt20//Simplify
15
New Data Transformations
RightMatrixFraction
SystemMatrixdt,TargetForm-gtRightFraction
SystemMatrixdt
16
A Least-Order form of a System-Matrix in
polynomial form is one in which there are no
input-decoupling zeros and no output-decoupling
zeros, and would yield a minimal state-space
realization, when directly converted to
state-space form. For example, the polynomial
system matrix
is not least order, since T(s) and U(s) have
3 input-decoupling zeros at s 0, 0, -1
i.e. T(s) U(s) has rank ? 4 at these values of
s.
17
Hence
18
System Analysis
Controllabless Observabless
ss A, B, C, D
Controllableps Observableps
Controllableps Observableps
MatrixLeftGCDT(s) U(s) MatrixRightGCDT(s)
V(s)
SmithFormT(s) U(s) McMillanFormG(s)
Decoupling Zeros
Decoupling Zeros
19
Controllability and Observability
In the same way that the controllability and
observability of a system described by a set of
state space equations can be determined in the
Control System Professional by entering the
commands Controllabless and Observabless where
ss is a StateSpace object. These tests can now
also be directly applied to a system matrix
object by entering the commands Controllablesm
and Observablesm where sm is a SystemMatrix
object in either polynomial form or state space
form.
20
Preliminary Analysis

• Reduction of State-Space Equations
• Given a system matrix in state-space form
• then an input-decoupling zeros algorithm,
• implemented in Mathematica, reduces P(s) to
• The completely controllable part is then given by

21
Gasifier
Model Format
A is 25 x 25 B is 25 x 6 C is 4 x 25 D is 4 x 6
Inputs- 1 char 2 air
3 coal 4 steam
5 limestone 6 upstream disturbance
Outputs- 1 gas cv 2 bed mass
3 gas pressure 4
gas temperature
22
Preliminary Analysis
The original 25th order system is numerically
very ill conditioned. The eigenvalues cover
a significant range in the complex plane, ranging
from -0.00033 to -33.1252. The condition number
is 5.24 x 1019. At w 0 the maximum and minimum
singular values are 147500 and 50,
respectively. The Kalman controllability and
observability tests yield a rank of 1, and the
controllability and observability gramians are -
23
Preliminary Analysis
Application of the decoupling zeros algorithm to
sI-A, B yielded
Dimensions of
Dimensions of
Dimensions of
indicating that the system had 7 input-decoupling
zeros, which was confirmed by transforming A and
B to spectral form.
24
Coprime Factorizations
25
Smith and McMillan Forms
The Smith form of a polynomial matrix and the
McMillan form of a rational polynomial matrix are
both important in control systems
analysis. Consider an ? x m polynomial matrix
N(s), then the Smith form of N(s) is defined
as S(s) L(s)N(s)R(s)
and L(s) and R(s) are unimodular matrices.
26
Smith and McMillan Forms
Consider now an ? x m rational polynomial matrix
G(s), and let G(s) N(s)/d(s) where d(s) is
the monic least common denominator of G(s), then
the McMillan form of G(s) is defined as
where M(s) is the result of dividing the Smith
form of N(s) by d(s), and cancelling out all
common factors
27
Design Methods
Synthesis Methods
Pole Assignment
PID Controller
Optimal Control
Nyquist Array
Model Ref. Opt. Control
Robust NA
Nonlinear Systems
Model-Order Reduction
28
Design/Synthesis Methods
Methods implemented are- 1 Pole Assignment -
Some Observations 2 Model-Reference Optimal
Control 3 The Systematic Design of PID
Controllers 4 Uncertain Nonlinear Systems 5
Robust Direct Nyquist Array Design Method 6
Model-Order Reduction
29
Pole Assignment

• We consider four main types of approaches
• ACKERMANNS FORMULA
• SPECTRAL APPROACH
• MAPPING APPROACH
• EIGENVECTOR METHODS

Control Systems Centre - UMIST
30
Ackermanns Formula
Here ? is the controllability matrix of A,b,
and pc(s) is the desired closed-loop system
characteristic polynomial.
Spectral Approach
where
Here, ?i and ?i are the open-loop system and
desired closed-loop system poles, respectively,
and the vi are the associated reciprocal
eigenvectors.
Control Systems Centre - UMIST
31
Mapping Approach
The state-feedback matrix is given as
where ? is the controllability matrix of A, b
Here, the ai and ?i are the coefficients of the
open-loop and closed-loop system characteristic
polynomials, respectively.
Control Systems Centre - UMIST
32
Eigenvector Methods
It is also possible to determine the
state-feedback pole assignment compensator as
where the mi are randomly chosen scalars and the
uci are the closed-loop system eigenvectors
calculated from
Selecting mi 1, for example, the
state-feedback compensator can be found as
Control Systems Centre - UMIST
33
Comparison of Dyadic Methods under Numeric
Considerations
Control Systems Centre - UMIST
34
Comparison of Dyadic Methods under Symbolic
Considerations
Control Systems Centre - UMIST
35
Model-Reference Optimal Control
Model Reference LQR System
The resulting optimal feedback controller Ko can
be partitioned as
where K11 and K21 operate on the reference-model
state vector xM and K12 and K22 operate on the
system state vector x.
36
Model Reference LQR feedback paths.
Model Reference LQR System Closed-Loop.
37
(No Transcript)
38
 Unit-step on reference input 1
Unit-step on reference input 2
 Unit-step on reference input 3
Unit-step on reference input 4
39
The PID Controller

• In recent years, several researchers have been
  re-examining the PID controller to determine the
  limiting Kp, Ki, and Kd parameter values to
  guarantee a stable closed-loop system namely,
• Keel and Bhattacharyya
• Ho, Datta, and Bhattacharyya
• Shafei and Shenton
• Astrom and Hagglund
• Munro and Soylemez

40
The PID Controller
41
 
Test compensator arrangement
 
 
 
Test compensator space 
42
 
 
The Nyquist plot for Kp 0.5 and Ki 0.5
The admissible PI compensator space
43
(No Transcript)
44
Design Requirements

• Stability
• Performance
• Robustness
• Simplicity
• Transparency

• increasing
• difficulty

45
Acknowledgements My thanks to Dr Igor Bakshee of
Wolfram Research for his interest and help in
carrying out this work.
46
D-Stability
Control Systems Centre UMIST
? Sin(?)
Im
?
D
Re
d
47
The Nyquist Plot Approach
Control Systems Centre UMIST

• Here, we detect 5 axis crossings,
  (-2,2,2,-1,-1), where the last is due to the
  infinite arc, on the right, due to the pole at
  the origin.

48
The Nyquist Plot Approach
The resulting stability boundary is
Note that the origin is not included in the
region because the basic system is unstable.
49
The Nyquist Plot Approach
Control Systems Centre UMIST

• Here, with Kp 5 and Ki 18 the system is
  stable, even with an additional gain of k
  1.3134
• yielding closed-loop poles
• -0.2519 5.4879i
• -1.2320 1.5258i
• -0.0161 0.4510i

50
Diagonal Dominance Concepts
Various definitions of Diagonal Dominance exist,
namely - Rosenbrocks row/column form
R Limebeers Generalised
Diagonal Dominance L Bryant Yeungs
Fundamental Dominance Y where the
conservatism of the resulting dominance criterion
reduces as Y lt L lt R
Mathematica Code
51
Nyquist Array Example
52
Gasifier
Model Format
A is 25 x 25 B is 25 x 6 C is 4 x 25 D is 4 x 6
Inputs- 1 char 2 air
3 coal 4 steam
5 limestone 6 upstream disturbance
Outputs- 1 gas cv 2 bed mass
3 gas pressure 4
gas temperature
53
Combined Sequential Loop Closing and Diagonal
Dominance Method

• This approach is a new combination of Bryants
  Sequential Loop Closing Approach with MacFarlane
  and Kouvaritakis ALIGN Algorithm, Edmunds
  Scaling and Normalization Technique, and
  Rosenbrocks Diagonal Dominance.
• It is particularly appropriate in cases where a
  simple controller structure is desired.
• Advantages
• It can be implemented by closing one loop at a
  time.
• Usually, the resulting control scheme is quite
  simple and can be easily realized in practice.

54
Achieving Diagonal Dominance

• Normalization
• Generates the input-output scaling to be applied
  to the system in order to minimize interaction.
• Determines the best input-output pairing for
  control purposes.
• Produces good diagonal dominance properties at
  low and intermediate frequencies.
• Results are obtained by using simple, wholly real
  permutation matrices.
• High frequency decoupling
• Aims at improving the transient response of the
  system.
• Emphasis is on frequencies close to the
  bandwidth, around which interaction is most
  severe.
• Results are obtained by making use of wholly real
  matrices.

55
Preliminary Analysis
The original 25th order system is numerically
very ill conditioned. The eigenvalues cover
a significant range in the complex plane, ranging
from -0.00033 to -33.1252. The condition number
is 5.24 x 1019. At w 0 the maximum and minimum
singular values are 147500 and 50,
respectively. The Kalman controllability and
observability tests yield a rank of 1, and the
controllability and observability gramians are -
56
Preliminary Analysis
Application of the decoupling zeros algorithm to
sI-A, B yielded
Dimensions of
Dimensions of
Dimensions of
indicating that the system had 7 input-decoupling
zeros, which was confirmed by transforming A and
B to spectral form.
57
(No Transcript)
58
(No Transcript)
59
(No Transcript)
60
(No Transcript)
61
Design Procedure - 1

• The Nyquist Array after an initial output scaling
  of diag0.00001 , 0.001 , 0.001 , 0.1 looks like
  

62
Design Procedure - 2

• The Nyquist Array after swapping the first two
  outputs (calorific value of fuel gas and bedmass)
  and closing the bedmass/char off-take loop is

63
Design Procedure - 3

• The Nyquist Array of the 3 x 3 subsystem after
  normalisation and high frequency decoupling at w
  0.001 rad/sec is (where the outputs are
  pressure, temperature and calorific value of fuel
  gas)

64
Design Summary

• Implement PI controller on bedmass/char-extraction
  loop.
• Scale inputs and outputs, to normalize them.
• Use ALIGN Algorithm for the remaining 3-input
  3-output subsystem.
• Design a PI controller for the fast Calorific
  Value Loop.
• Design a PI controller for the fast Pressure
  Loop.
• Design a Lag-Lead controller for the remaining
  slow Temperature Loop.

The control scheme resulting from this approach
is as follows
65
Controller
Constant Pre-compensator Constant
Post-compensator
Dynamic Controller

66
Model Simplification
67
Model Simplification
Root-locus diagram of g1,1(s)
68
Theorem By using just the first input of a
given MIMO system, it is almost always possible
to arbitrarily assign ?1 self conjugate poles of
the system, and make these poles uncontrollable
from the other inputs, provided that the
system A, b1,C has ?1 controllable and
observable poles, where b1 is the first column of
the input matrix (B), where
This result can be compared with a previous
result developed by Munro and Novin-Hirbod (1979)
for the case of dynamic output feedback, where
the degree r of the necessary compensator is
given by
Control Systems Centre - UMIST