DESIGN AND ANALYSIS OF MULTIVARIABLE PREDICTIVE CONTROL
APPLIED TO AN OILWATERGAS SEPARATOR: A POLYNOMIAL APPROACH
By
GIOVANI CAVALCANTI NUNES
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
2001
Copyright 2001
by
Giovani Cavalcanti Nunes
To my beloved daughters Julia and Giovanna and to my mother Nisete.
ACKNOWLEDGMENTS
I would like to thank Professor Oscar Crisalle for his help. His insight and
knowledge of Process Control were fundamental in the development of this dissertation.
I consider myself a lucky person for finding such a wonderful advisor.
I would also like to thank professor Santiago for his immense help and support in
solving all kinds of problems that naturally arise to a foreign student in graduate school.
I gratefully acknowledge the support received from my employee Petrobras for
this work. The conditions were excellent. I especially thank the engineer and leader Jose
Antonio de Figueiredo, head of the Division of Offshore Basic Engineering, at the time,
to whom this thesis owes its existence. His energy and ability to motivate people are
outstanding.
Professor Antonio Coelho, a colleague and a friend, is another person I must
mention for his help and solidarity during the long hours of work in the control
laboratory. Excitement and enthusiasm were always present when Antonio "Banderas"
was around.
I would finally like to thank the big group of wonderful friends I had the chance
to meet while in Gainesville. I would like to list all their names here but this is
unfeasible. Kerry Alahar, however, must be mentioned for our friendship proved to be a
special one. Our common interests, experience and burdens during graduate school
created bonds that will last forever.
TABLE OF CONTENTS
page
A C K N O W L E D G M E N T S ........................................................................ .................... iv
A B ST R A C T ................................v ii............................
CHAPTERS
1 IN TR O D U C T IO N .................................. ........................................ ................. ...
2 L ITE R A TU R E R E V IE W .................................................... .................................. 7
2.1 H history of M odel Predictive Control .................................................. .... .......... 7
2.2 Stability ....................... ........... ....... ............. ..... .. ...................... 9
2.3 Model Predictive Control: An Optimal Control Problem ................................ .... 10
3 FORMULATION AND SOLUTION OF THE CONTROL PROBLEM ................ 12
3.1 Objective Function of Multivariable Predictive Control.................. ................ .... 12
3.2 Derivation of the Predictive Control Law .................................... ..............16
3.3 D erivation of the Constant Forcing Function................................... .................... 19
3.4 F inal C control L aw E quation ............................................................................ .... 2 1
3 .5 C onclu sion.................................................. ....... ......... ...... 22
4 ANALYSIS OF CLOSEDLOOP BEHAVIOR ................................................... 23
4.1 Stability Analysis of Closedloop Multivariable Systems ........................................23
4.2 Closedloop Transfer Function M atrices.................. ............................................. 24
4.3 MFD of Transfer Function Matrices ................................ .............. 26
4.4 Poles and Zeros of a Transfer Matrix................................ .................. 28
4.5 Perform ance A analysis ................................................... ................................. 38
4 .6 Im plem entation ............................................... .............. ............ .. ...... .... 42
4.7 C conclusion ............................. .............. ...... 44
5 THREEPHASE SEPARATOR ........................................ .......................... 45
5.1 D description of the Separator................................................. .......................... 46
5.2 Rigorous M odel ....... ...... ................. .............................. .............47
5.3 Sim plified M odel .......... .... ................................... ............ .. 48
5.4 F low E quations.............................................. ....... ......... ...... 5 1
5.5 D differential Equations ............................................ .... ...... ........ ........ .... 51
5.6 C onclu sion.................................................. ....... ......... ...... 52
5.7 Nomenclature ............................................ 53
6 DESIGN OF A PREDICTIVE CONTROLLER FOR THE SEPARATOR ............. 55
6.1 Discrete Linear M odel of the Separator ................ ....... ........ .................... 55
6.2 D esign of Predictive C ontroller........................................................ .... ................ 57
6.3 Case 1: First Tuning .. ................. ........................ .... .......... .. 58
6.4 Case 1: Second Tuning.................... ......................................... .......................... 67
6.5 C ase 2 .................................................................................. ......... ...............................75
6.6 C onclu sion ............................................ 8 1
7 SIM ULATION AND RESULTS ....................................................... .................. 82
7.1 Results of Case 1 ............................................ .......... 82
7.2 R results of C ase 2 ................ ............................................ ......... ... 86
7 .3 C onclu sion.................................................. ....... ......... ...... 89
8 C O N C LU SIO N S .............................. ..................................................... 90
APPENDICES
A. DETAILED DERIVATION OF CONTROL LAW EQUATION........................... 92
A. 1 Derivation of Matrix Diophantine Equation ..................................................92
A.2 Derivation of Predictor Equation .......................... ......... .. ...... ......... 95
A.3 Simplification of Predictor Equation .................................. .............. 97
A.4 Derivation of Constant Forcing Response Equation............................................. 100
A. 5 Substitution of Constant Forcing Response in the Control Law Equation............. 100
B. DERIVATION OF PRODUCT KP AND KF....................................... ............... 104
B P rodu ct of K P ............................................ ......... ................ ............................... 104
B .2 Product of K F .................. .......... .. ................. .. ......... ........ ..... 106
C. OBTAINING A LEFT MATRIX DESCRIPTION .............................................. 107
D. CLOSEDLOOP TRANSFER MATRICES.................... .......... ............. 110
D 1 Disturbances in the Closedloop........................................................... .............. 110
D.2 Closedloop Transfer M atrices of Predictive Control............................................. 112
E. OBTAINING M ATRIX K ............................ ................................. .............. 113
L IST O F R E F E R E N C E S ................................................................................................ 115
BIO GRAPH ICAL SKETCH .................................................. ............. 118
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
DESIGN AND ANALYSIS OF MULTIVARIABLE PREDICTIVE CONTROL
APPLIED TO AN OILWATERGAS SEPARATOR: A POLYNOMIAL APPROACH
By
Giovani Cavalcanti Nunes
August 2001
Chairman: Tim Anderson
Major Department: Chemical Engineering
This dissertation uses a polynomialoperator technique to study stability and
performance of unconstrained multivariable predictive control. A simple and direct way
to determine stability of the closedloop system is developed.
It is shown that to guarantee stability of the closedloop two transfer matrices
must be stable. These are represented as fraction descriptions, a ratio of "numerator"
and "denominator" polynomial matrices from which the poles can be determined.
Because both transfer matrices possess the same denominator matrix the position of the
roots of its determinant give sufficient conditions for stability of the closedloop.
Furthermore it is shown that if a coprime fraction description is done for the process
transfer matrix then it is necessary and sufficient to check if the roots of the determinant
of the denominator matrix lie inside the unit circle. This technique avoids the inversion
of transfer matrices which is a numerically difficult task allowing the tuning of
multivariable systems with many inputs and outputs.
Performance is also studied and it is proven that the system has zero offset
response to step changes in the reference, a property known to be valid for the single
input singleoutput case. For systems with an equal number of inputs and outputs the
"inversion of the plant" is also proven. In this case the weights on the input are zero and
the solution to the optimization problem results in a controller that inverts the plant and
the output matches the reference.
The use of a multivariable predictive control for an oilwatergas separator is
studied. The nonlinear model of the plant is linearized around a steady state and different
predictive controllers are designed for it. The controller responds positively to changes in
the parameters and performance objectives are pursued. Results show agreement with
the simulations done for the linear model, it is concluded that predictive control is a
successful control strategy for oilwatergas separators.
This method becomes an important tool for the analysis of predictive controllers,
allowing the study of the effect of tuning parameters on the behavior of the system when
the constraints are removed.
CHAPTER 1
INTRODUCTION
The process industry is characterized by ever tighter product quality
specifications, increasing productivity demands, new environmental regulations and fast
changes in the economical market. Over the last two decades predictive control has
proven to be a very successful controller design strategy, both in theory and practice.
The main reason for this acceptance is that it provides high performance controllers that
can easily be applied to difficult highorder and multivariable processes. Process
constraints are handled in a simple and systematic way. However a general stability and
robustness theory is lacking.
Predictive control is a class of control strategies based on the explicit use of a
process model to generate the predicted values of the output at future time instants, which
are then used to compute a sequence of control moves that optimize the future behavior
of a plant. Predictive control is rather a methodology than a single technique. The
difference in the various methods is mainly the way the problem is translated into a
mathematical formulation.
The explicit use of a model is the main difference between predictive control and
the classical PID controller. Its advantage is that the behavior of the controller can be
studied in detail, simulations can be done and performance can be evaluated. One of the
drawbacks is the need of an appropriate model of the process. The benefits obtained are
affected by the discrepancies existing between the real process and the model. According
to some researchers 80% of the work done is in modeling and identification of the plant.
The results almost always show that the effort is paid back in short time. Another
drawback is that although the resulting control law is easy to implement and requires
little computation, its derivation is more complex than that of the PID. If the process
dynamics do not change, the derivation can be done beforehand, but in the adaptive
control case all the computation has to be carried out at each sampling time. When
constraints are considered, the amount of computation is even higher.
Predictive control was pioneered simultaneously by Richalet et al. (1978) and
Cutler and Ramaker (1980). The use of finiteimpulse response models and finitestep
response models, which are easily obtained for open loop stable processes, partly
explains the acceptance in the process industry.
'"'" (t+i t "'  .. .
y(t)
y(t) j (t +i t) "
u(t) u(t + i)
t3 t2 t1 t t+ 1 t+ 2 ... t+Ny
Figure 1: Predictive control strategy
The methodology of all predictive controllers consists in predicting, at each time
t, the future outputs for a determined horizon Ny, called the prediction horizon. This
prediction of the outputs, t(t +i I t) for i = 1... Ny, is based on the model of the process
and depends on the known values of past inputs and outputs up to instant t.
The set of future control signals is calculated by optimizing a given criterion
(called objective function or performance index) in order to keep the process as close as
possible to the reference trajectory r(t+i), which can be the set point or an
approximation of it. Different algorithms present different forms of objective functions,
which usually take the form of a quadratic function of the errors between the predicted
output signal and the predicted reference trajectory. Some algorithms use the states of
the process as opposed to the outputs. The control effort is included in the objective
function in most cases. Weights are used to adjust the influence of each term in the
equation. The solution to the problem is the future control sequence that minimizes the
objective function equation. For that, a model is used to predict future outputs or states.
A typical objective function equation of a singleoutput singleinput process is
Ny2 Nu
J= (i)(r(t+i) (t+i It))2 + i)Au(t+i1)2
i=Nyl I=1
It is a quadratic function of future inputs, u(t + i), and the error between future values of
reference r(t+i) and predicted outputs, t(t+i t). Weights 0 and k are used to adjust
the influence of the error and inputs respectively.
The sequences of predicted outputs and future inputs are limited to horizons Ny2
and Nu respectively. The limitation on the sequence of inputs, u(t+i) from i=1 to Nu,
comes from the assumption that the control action is constant after Nu steps ahead. The
prediction horizon, on the other hand, limits the sequence of predicted output considered
in the objective function equation. The control horizon has to be smaller than the
prediction horizon. Weights and horizons are tuning parameters of the controller.
The optimization of the objective function requires a prediction of the future
outputs. The predicted output is the addition of two signals:
(t +i I t) = y0(t + i I t)+ GAu(t +i)
the constant forcing response and forced response. The constant forcing response,
y(t+ilt), corresponds to the prediction of the evolution of the process under the
consideration that future input values will be constant. The forced response, GAu(t+i)
where G is the dynamic matrix of the process, corresponds to the prediction of the output
when the control sequence is made equal to the solution of the minimization of the
objective function.
The expression for the predicted outputs can be substituted in the objective
function and the solution to the minimization problem leads to the desired control
sequence. Once the control sequence has been obtained only the first control move is
implemented. Subsequently the horizon is shifted and the values of all sequences are
updated and the optimization problem is solved once again. This is known as the
Receding Horizon Principle and is adopted by all predictive control strategies. It is not
advisable to implement the entire sequence over the next Nu intervals because it is
impossible to perfectly estimate the unavoidable disturbances that cause the actual output
to differ from the predictions made. Furthermore the operation might decide to change
the set point over the next Nu intervals.
The various predictive control algorithms only differ themselves in the model
used to represent the process, the model for the noise and the objective function to be
minimized. The models used can be Impulse/Step response models, transfer function
models or state space models.
* Impulse/Step response models: Predictive control has part of its roots in the process
industry where the use of detailed dynamical models is uncommon. Getting reliable
dynamical models of these processes on the basis of physical laws is difficult and the
first models used were Impulse/Step response models. The models are easily
obtained by rather simple experiments and give good results. The disadvantage is
that it requires a large amount of parameters.
* Transfer function models: For some processes good models can be obtained on the
basis of physical laws or by parametric system identification. In this case a transfer
function model is preferred. Less parameters are used than in the Impulse/Step
response models.
* State Space models: These models are the most general description for a linear time
invariant system and being so are the most elaborate and difficult to generate.
Generalized Predictive Control (GPC) is a class of predictive control, proposed by
Clarke et al. (1987). It provides an analytical solution (in the absence of constraints) and
can deal with unstable and nonminimum phase plants. It uses transfer function model
for the plant leading to lowerorder representations.
Different transfer function models are used by GPC such as Controller Auto
Regressive Integrated Moving Average (CARIMA), Controller AutoRegressive Moving
Average (CARMA), etc. In this dissertation a Deterministic AutoRegressive Moving
Average (DARMA) model is considered
A(q')y(t) = B(q')u(t 1) (1.1)
which can be obtained by a left matrix fraction description (see Appendix C) of the
transfer function model.
Most industrial plants have many outputs and manipulated variables. In certain
cases a manipulated variable mainly affects the corresponding controlled variable and
each of the inputoutput pairs can be considered as a singleinput singleoutput (SISO)
plant and controlled by independent loops. In cases where the interaction between the
different variables is not negligible the plant must be considered as a multipleinputs and
multipleoutputs (MIMO) process. These interactions may result in poor performance
and even instability if the process is treated as multiple SISO loops.
In practice all processes are subject to restrictions and these can be considered in
the objective function equation as constraints on the inputs and outputs. In many
industrial plants the control system will operate close to the boundaries due to operational
conditions. This has lead to the widespread belief that the solution to the optimization
problem lies at the boundaries. On the other hand the design of multivariable predictive
controllers requires the specification of horizons and weights for all inputs and outputs.
Therefore a big number of parameters must be selected, a task that may prove
challenging even for systems of modest size. Badly tuned predictive controllers tend to
take the plant to extreme and sometimes unstable conditions, frequently unnoticed
because of physical constraints such as valve opening, maximum flow rate, etc., masking
the instability problem and fooling operational groups into concluding that the plant has
been optimized.
Many different predictive control algorithms have treated multivariable systems
in some way or another but none has established the mathematical basis for a
comprehensive analysis of the system, such as performance or stability.
In this dissertation a practical method for analyzing the stability of classical
multivariable predictive control systems is pursued. The focus is on unconstrained
systems in order to allow the identification of the destabilizing choices of tuning
parameters, avoiding the masking effect introduced by the presence of hard constraints.
CHAPTER 2
LITERATURE REVIEW
Model predictive control (MPC) has had an enormous impact in the process
industry over the last 15 years. It is an effective means of dealing with multivariable
constrained control problems and the many reports of industrial applications confirm its
popularity.
2.1 History of Model Predictive Control
The development of modern control concepts can be traced back to the work of
Kalman in the early 1960s with the linear quadratic regulator (LQR) designed to
minimize an unconstrained quadratic objective function of states and inputs. The infinite
horizon endowed the LQR algorithm with powerful stabilizing properties. However it
had little impact on the control technology development in the process industries. The
reason for this lies in the absence of constraints in its formulation, the nonlinearities of
the real systems, and above all the culture of the industrial process control community at
the time, in which instrument technicians and control engineers either had no exposure to
optimal control concepts or regarded them as impractical. Thus the early proponents of
MPC for process control proceeded independently, addressing the needs of the industry.
In the late 1970s various articles reported successful applications of model
predictive control in the industry, principally the ones by Richalet et al. (1978) presenting
Model Predictive Heuristic Control (later known as Model Algorithmic Control (MAC))
and those of Cutler and Ramaker (1980) with Dynamic Matrix Control (DMC). The
common theme of these strategies was the idea of using a dynamic model of the process
(impulse response in the former and step response in the later) to predict the effect of the
future control actions, which were determined by minimizing the predicted error subject
to operational restrictions. The optimization is repeated at each sampling period with
updated information from the process. These formulations were algorithmic and also
heuristic and took advantage of the increasing potential of digital computers at the time.
Stability was not addressed theoretically and the initial versions of MPC were not
automatically stabilizing. However, by focusing on stable plants and choosing a horizon
large enough compared to the settling time of the plant, stability is achieved after playing
with the weights of the cost function.
Later on a second generation of MPC such as quadratic dynamic matrix control
(QDMC; Garcia, Morshedi, 1986) used quadratic programming to solve the constrained
openloop optimal control problem where the system is linear, the cost quadratic, the
control and state constraints are defined by linear inequalities.
Another line of work arose independently around adaptive control ideas
developing strategies essentially for monovariable processes formulated with transfer
function models (for which less parameters are required in the identification of the
model) and Diophantine equation was used to calculate future input. The first initiative
came from Astron et al. (1970) with the Minimum Variance Control where the
performance index to be minimized is a quadratic function of the error between the most
recent output and the reference (i.e. the prediction horizon Ny=l). In order to deal with
nonminimum phase plants a penalized input was placed in the objective function and
this became the Generalized Minimum Variance (GVM) control. To overcome the
limitation on the horizon, Peterka (1984) developed the PredictorBased SelfTuning
Control. Extended Prediction SelfAdaptive Control (EPSAC) by De Keyser et al.
(1985) proposes a constant control signal starting from the present moment while using a
suboptimal predictor instead of solving a Diophantine equation. Later on the input was
replaced by the increment in the control signal to guarantee a zero steadystate error.
Based on the ideas of GVM Clarke et al. (1987) developed the Generalized Predictive
Control (GPC) and is today one of the most popular methods. A closed form for the GPC
is given by Soeterboek. Statespace versions of unconstrained GPC were also developed.
2.2 Stability
Stability has always been an important issue for those working with predictive
control. Due to the finite horizon stability is not guaranteed and is achieved by tuning the
weights and horizons. Mohtadi proved specific stability theorems of GPC using state
space relationships and studied the influence of filter polynomials on robustness
improvement. However a general stability property for predictive controllers, in general,
with finite horizons was still lacking. This lead researchers to pursue new predictive
control methods with guaranteed stability in the 1990s. With that purpose a number of
design modifications have been proposed since then including the use of terminal
constraints (Kwon et al., 1983; Meadows et al. 1995), the introduction of dualmode
designs (Mayne and Michalska, 1993) and the use of infinite prediction horizons
(Rawlings and Muske, 1993), among others. Clarke and Scattolini (1991) and Mosca et
al. (1990) independently developed stable predictive controllers by imposing endpoint
equality constraints on the output after a finite horizon. Kouvaritakis et al. (1992)
presented a stable formulation for GPC by stabilizing the process prior to the
minimization of the objective function. Many of these techniques are specialized for
statespace representations of the controlled plant, and achieve stability at the expense of
introducing additional constraints and modifying the structure of the design.
Practitioners, however, avoid changing the structure of the problem and prefer to achieve
stability by tuning the controller. For that a good doses of heuristics is used.
2.3 Model Predictive Control: An Optimal Control Problem
Recently a theoretical basis for MPC has started to emerge. Researchers have
revisited the LQR problem arguing that model predictive control essentially solves
standard optimal control problems with receding horizon, ideas that can be traced back to
the 1960s (Garcia et al., 1989). MPC is characterized as a mathematical programming
problem since the slow dynamics of process industry plants allow online solution to
openloop problems in which case the initial state is the current state of the system being
controlled. Determining the feedback solution, on the other hand, requires the solution of
the HamiltonJacobiBellman (a dynamic programming problem) differential or
difference equation which turns out to be more difficult. Riccati equation appears as a
particular case to some optimal control problems. It is shown (Mayne et al., 2000) that
the difference between MPCs approach and the use of dynamic programming is purely
one of implementation. This line of research has an early representative in the work of
Kwon and Pearson (1983) and Keerthi and Gilbert (1988) and has recently gained
popularity through multiple advocates, such as the work of Muske and Rawlings (1993).
More recently two major approaches to the control problem are very popular. The
first one employs the optimal cost function, for a fixed horizon, as a Lyapunov function
and the second approach takes advantage of the monotonicity property of a sequence of
the optimal cost function for various horizons. Note that for linear systems the presence
of hard constraints makes the design of the controller a nonlinear problem so that the
natural tool for establishing stability is Lyapunov theory.
Obtaining a formulation which can relate the tuning parameters of MPC to
properties such as stability and performance is the major goal of this line of work and
recent advances show promising results. Solving the dynamic programming problem,
however, is not practical and this can be accounted for the resistance of the industry to
adopt these new methods.
CHAPTER 3
FORMULATION AND SOLUTION OF THE CONTROL PROBLEM
In this chapter the predictive control law is derived for an unconstrained MIMO
system with m controlled outputs and p inputs. The problem of minimizing the objective
function is solved analytically using least squares. The formulation is done in vector
matrix format. A final expression, function of future values of reference and output, is
reached. Diophantine equation is used to generate the predicted values of the output.
3.1 Objective Function of Multivariable Predictive Control
Unconstrained multivariable predictive controllers are designed to minimize the
objective function:
Ny1 NY2
J (i)(r (t+i)(t+i I t))2 + 2 (i)(r( +i) (t + I t))2 +...
=1 =1
Nym Nu1
+) (i)(r,(t+i)_ (t+i t)) + (i) (t+i1) (3.1)
=1 =1
Nu2 Nup
+1 2 (i)AU2 (t +i1)2 +... + P (i)Ap (t + i 1)2
I=1 I=1
where:
y, is thejth output
r is the desired trajectory for thejth output
u, is the th input
J are the weights on the error between the desired trajectory and the output
k, are the weights on the control signal
Ny, are the prediction horizons
Nu, are the control horizons
The general aim of the objective function is to make the future output, on the
considered horizon, follow a desired reference signal and, at the same time, penalize the
control effort (Au) in order to avoid unfeasible inputs. Many formulations of predictive
control consider an initial and a final prediction horizon to account for the dead times of
the process. In this dissertation the initial prediction horizon considered is always 1 and
if dead time is to be considered then the weights can be appropriately set to zero.
Similarly to the SISO case, where the prediction horizon has to be greater than the
control horizon, the sum of all prediction horizons has to be greater than the sum of all
control horizons; Ny 2 Nu where
m p
Ny A NyJ and Nu AL Nu,
J=1 f=1
Note that no control horizon should be greater than the max prediction horizon:
Nu, < max Ny, for allj
The many different indices required in the derivations that follow makes notation
an important issue. The following rules are applied hereon:
polynomials are written in upper case letters
vectors in bold lower case letters
matrices in bold upper case letters and in some exceptional cases bold lower case
letters.
Vectors and matrices are created to represent the time expansion of the equations.
Whenever the sequence of values from 1 to i has to be represented in a vector then t + i
will be written in between parenthesis. Its absence is indicative of a sequence of values
from 1 to the given horizon. As an example, let us represent the values of .5 from step 1
to i and also from step 1 to Nyi:
.A (t +1 t) 7 (t +1 t)
1 ,(t+2 t) ,(t + 2 1t)
y (t+i t)= y = y (t+Nyz\t)=
(t +i t) J (t+Ny,1t)
Notice the difference between the scalar j, (t + i I t) and the vector j, (t + i I t) The former
indicates the value of A, at time step t + i while the latter denotes the vector of values of
.5, from time step 1 to i. Subscripts are used in matrices for the same purpose.
The objective function can be written in the more elegant vectormatrix format.
This will simplify the notation and help visualization of the solutions.
Claim 3.1. Equation (3.1) can be written in compact vectormatrix notation as
J= (r y)TQ(r ) + Au AAu (3.2)
where the augmented vectors y and r e 9Ny and u e 9" are given by
y2 r, M2
y= r= u= (3.3)
and where matrices Q~e 9Ny Ny and A E 9^Nu"N are given by
O ... iA, O ... O
S 0 0. 0 ... 0
D = A= (3.4)
0 0 ... ... AP
where in turn
,j (t +1 1 t) r (t + 1)
,Q(t + 2 1 t) r(t + 2)
, (t + Ny, t) r (t + Nyj)
0 ... 0
01(2) ... 0
0 .. O (Ny, )
u,(t)
u, (t +1)
U j,= +
u (t + Nut 1)
0 ... 0
kj(2) ... 0
: : (3.6)
0 ... (Nu, )
and y, and r e 9NyJ, u, e ~i"', E e h' yj' and A, e ~~i ,XNu,
Proof. From definitions (3.5) to (3.6), the summations in (3.1) can be written as
Ny1
(i)((t+i)_(t +i I t))2 =(r j),), (r, _ )
1=1
Ny2
2 (i)(r, (t + i) ,(t + i t))2
1=1
= (r2 j,)T Q (r2 y2)
Ny,
(i)(r[(t + i)y (t + i t))2 = ( j)T (r )
1=1
and
Nu1
k,) (i)Au1 (t + i)2
= Au1 A Au1
Nu2
2(i)Au, (t + i)2 = Au2A Au2
1=0
Nup
k p(i)Aup (t +i)2 = AuT A Au
I=0
j (1)
0
0
(3.5)
,(1)
0
A,=
0
substituting the above equations in (3.1) and using the definitions of vectors and
matrices in (3.3) to (3.4) it follows that the objective function equation can be written
as J = (r )')/(r .) + Au'AAu. Q.E.D.
3.2 Derivation of the Predictive Control Law
The problem of minimizing J with relationship to Au requires an expression for
the predicted outputs as a function of future inputs.
Claim 3.2. Future values of the output, y, are the sum of the constant forcing response,
y0, and the response to future inputs
GAu:
y = yO +GAu
where
0 G
y11
y G
y, G,,
where in turn
y (t+1 t) g0 0
Y J G *
y (t+N2yt) G Ny= N 2
and Ge 9`Nu", y~O ^9V, YJ E 9^Y. Each dynamic matrix
corresponding control and prediction horizons, i.e. GJ, e 9N "N
Proof Consider the future values of the first output of the system
)1 = yo +G1Au1 +G12Au2 +...+G ,Aup
0
0
gNy,Nu d to
is truncated to the
(3.8)
(3.7)
 GmP
*  Gp
where y, is the constant forcing response and G,, are the dynamic matrices relating
input to output 1. Similar equations can be written for the remaining outputs of the
MIMO system as follows
2 = y2 + 21Au, + G2 AU2 +... + G2p Au
my = y + GmAu, + G2Au, +... + GmpAup
The above equations can be put together as
Y y1 GI G,, ... GIp Au,
y y G G2, *... G21 Au,
2 + (3.9)
m yO Gm Gm2 . Gp Aup
from where we conclude that
y = yO +GAu Q.E.D.
This is the equation needed for an analytical solution of the objective function
equation since it shows, explicitly, the influence of future inputs in the future outputs.
The solution is a function of the future reference trajectory and the constant forcing
response as seen in the next claim.
Claim 3.3. Given the objective function equation J=(r y)'(r )+ Au'AAu and
the prediction equation for future outputs, y = yO + GAu, then the optimal control
sequence is given by
Au = (GTOG + A)'G TO(r yO) (3.10)
Proof Substitute equation (3.7) in the objective function equation to get
J = (r yO GAu) ((r y GAu) + Au'AAu
Let D = RTR
J = (r yO GAu)TRR(r y GAu) + Au AAu
then
J = (Rr Ry' RGAu)T (Rr Ry' RGAu) + Au AAu
differentiating with respect to Au
jJ
= 2(Rr Ry RGAu)'RG + 2Au'A
3Au
equating the derivative to zero in order to get the minimum of the function then
2(Rr Ry RGAu) RG + 2AuA = 0
and
Au'A = (r y GAu)TR RG
AuA = (r y')T GAuT G'TOG
Au' (A + G'TG) = (r yO)TOG
transposition yields
(A+GTOG)T Au= GOT (r yO)
Note that (A+ G'TOG) and D are symmetric matrices (making the transposition
irrelevant). Thus the optimal control sequence is given by
Au = (GTOG + A)'G TO(r yO) Q.E.D.
From the receding horizon principle only the part corresponding to the first
control move of matrix (G'rOG+A)'G'r is used. That way matrix
(G'rOG + A)'GTO can be replaced by a smaller dimension matrix.
Claim 3.4. The control law for the multivariable predictive control is given by
Au(t) = K(r y) (3.11)
where matrix K e 9~p is defined as
k[1 k ... kTm
k k T ...
K 21 22 2m (3.12)
k k ...k
p1 p2 pm
and its elements
k =[k, k2 ... kN,, (3.13)
are the rows of (G'ROG + A)'GT corresponding to the first control move.
See Appendix E for the proof.
3.3 Derivation of the Constant Forcing Function
In order to solve the control law equation, u = K(r yO), an expression for future
values of the constant forcing response, y0, is needed. This can be done via the solution
of a series of Diophantine equations. For SISO cases the developments are widely
documented in the literature (Clarke et al., 1987; Clarke and Mohtade, 1989; Crisalle et
al., 1989). The derivation for the MIMO systems will be done here. The approach
consists of first extending the SISO results to the MISO plants, and then applying the
MISO approach to each of the outputs of the MIMO plant. The individual MISO
predictors are then stacked into augmented vectors and matrices as needed to build the
final MIMO operators sought, as was done in section 3.2 in the derivation of the
predictive control law. The many different equations for the MIMO case make the
derivation quite cumbersome and only the vectorial equations will be shown here while a
detailed derivation is found in Appendix A. The definitions for the terms of the equations
that follow will not be repeated here and the reader should refer to the appendix
whenever necessary.
Consider the multivariable system defined by the DARMA model of equation
(1.1), A(q')y(t)= B(q')u(t1). The predicted values of the output can be generated
using Matrix Diophantine equation,
E (q')AA(q')+ q' (q) = I, (A.9)
where E (q1) and 1F(q1) are unique polynomial matrices defined in Appendix A. For
that, multiply equation (1.1) by E (q')A to get
E, (q')AA(q')y(t) = E, (q')B(q')Au(t 1)
From Matrix Diophantine equation we see that E,(q')AA(q') =I q 'F(q') and
consequently
[I, q'F (q') y(t) = E,(q')B(q')Au(t1) (3.14)
Simple manipulation leads to y(t) = E (q)B(q')Au(t 1)+ q' (q')y(t). Multiply by
q' to get an expression for predicted values of the output
y(t + i) = E (q')B(q')Au(t + i1)+ F (q')y(t) (3.15)
Recognizing that
E (q')B(q') = G, (q') + q'"(q') (3.16)
then
(3.17)
j(t+ i) = F, (q')y(t) + G, (q')Au(t + i 1)+ q'P (q')Au(t)
Finally, verifying that the predicted values of the output, )(t+i), are the unforced
response of the system, yo(t+i), when future values of the input are kept constant, i.e.
Au(t+i1)=O for i>0, then
yO (t + i) = F (q')y(t)+ q', (q')Au(t) (3.18)
or if we consider the whole prediction and control horizons
yO = F(q')y(t)+q'P(q')Au(t)
3.4 Final Control Law Equation
The control law equation (3.11) can be written, according to the notation adopted,
in terms of a hypothetical sample time t + i as
Au(t)= K, [r(t +i) yo(t + i)]
and substituting the constant forcing function of equation (3.18) then
Au(t)= K, [r(t + i) F, (q')y(t)q' 1, (q')Au(t)]
which can be rewritten as
[I+ q'K,P,(q')] Au(t)= Kq'r(t) KF, (q')y(t)
Define
R (q) = [I+ qyK P(q)] A (3.19)
S,(q') = K,F,(q') (3.20)
T(q)= K,(q)q1 (3.21)
to get
R, (q')u(t) = T (q)r(t) S, (q )y(t) (3.22)
Considering the prediction horizons it yields
R(q')u(t) = T(q)r(t) S(q')y(t) (3.23)
Note that R(q') and S(q) are polynomial matrices in backward shiftoperator while
matrix T(q) is in forward shiftoperator.
A Matlab code was created to generate matrices R(q'), S(q') and T(q). Once
the set of parameters is specified, namely prediction horizons, control horizons, weights X
and 4, then matrices R(q'), S(q') and T(q) are generated for a given linear model of
a process plant.
3.5 Conclusion
In this chapter a concise mathematical development of the control law for the
multivariable predictive control is shown. The appropriate definition of vectors and
matrices allowed a mathematical formalism similar to the SISO case and a design
procedure is derived from it. A Matlab code has been developed for the design of
multivariable unconstrained predictive controllers.
CHAPTER 4
ANALYSIS OF CLOSEDLOOP BEHAVIOR
In this chapter the closedloop system is analyzed for stability and performance.
For that a practical method to determine the closedloop poles of multivariable predictive
control is searched. Proof of properties known to be valid for the SISO case are pursued
for the multivariable case.
4.1 Stability Analysis of Closedloop Multivariable Systems
The issue of ensuring stability of a closedloop feedback controller is of utmost
importance to control system design. Albeit the increasing popularity of predictive
control over the last decades a general method for the analysis of stability of the resulting
closedloop is still not available for multivariable systems. Although unconstrained
systems present analytical solutions to the problem of minimizing the objective function,
obtaining the closedloop transfer matrices is a major challenge. Numerical difficulties in
inverting transfer matrices have limited the size of systems for which the poles can be
calculated.
A stable system is a dynamic system with a bounded response to a bounded input.
Many descriptions of stability exist, namely external stability, internal stability, marginal
stability and exponential stability. In this dissertation external (also known as BIBO for
bounded input bounded output) stability is analyzed.
During the design phase of predictive control no a priori guarantee of stability
exists for the resulting closedloop system. Stability of a system is determined by the
position of its poles in the complex plane, which can be determined, among other ways,
by analyzing the individual transfer functions of the closedloop matrix. Matrix fraction
descriptions can also be used to decompose rational matrices into 'numerator' and
'denominator' matrices allowing for the determination of poles and zeros of multivariable
systems. Under certain circumstances the latter method proves to be better.
Since the aim of this chapter is to study stability of the closedloop system, in the
sequel, the complex variable z is substituted for the forward shift operator q. This
substitution is allowed by the analogy between shiftoperator theory and transfer function
theory for linear systems with zero initial conditions. Hence matrix R(z') is substituted
for matrix R(q') and so on. Furthermore, for simplicity of notation, the arguments z'
and z are omitted whenever the interpretation of the transfer function in question is
deemed unambiguous.
4.2 Closedloop Transfer Function Matrices
Figure 4.1 shows the block diagram of the multivariable predictive control applied
to a plant H(z).
Figure 4.1 Closedloop System
Let H(z) be the process plant. Using blockdiagram algebra operations, the
dynamics of the system shown in Figure 4.1 can be captured by the closedloop transfer
matrices
G, (z) = [R(z)+ S(z)H(z)]1 T(z) (4.1)
and
G,, (z) = H(z)[R(z) + S(z)H(z)]1 T(z) (4.2)
Transfer function matrix T(z) is a polynomial matrix applied to the future reference
trajectory. Hence w(z) is a bounded signal for a bounded reference signal r(z) and if
all possible disturbances and noises are considered then it is seen (refer to Appendix D)
that for the analysis of stability it suffices to analyze the transfer functions relating the
control signal u(z) to w(z) and the output y(z) to w(z); G,,(z) and GY (z)
respectively.
G, (z) = [R(z)+ S(z)H(z)]1 (4.3)
and
G, (z)= H(z)[R(z)+S(z)H(z)]1 (4.4)
THEOREM 4.1 Transfer function matrix R(z) is biproper.
Proof Recall that a proper transfer matrix whose inverse is also proper is called biproper.
From the definitions of K and P(z) its immediate to verify that the product z'KP(z) is
a proper transfer matrix. The addition of the identity matrix to it, I + z'KP(z), causes
the determinant to be biproper and as a consequence R(z) = I+z'KP(z)]A is a
biproper transfer matrix. Q.E.D.
This property of R(z) makes R(z)' biproper and allows the implementation of
the controller.
THEOREM 4.2 Let H(z) be a proper transfer matrix, then G (z) is biproper and
G,(z) is proper.
Proof Matrix S(z) is the product of K and F(z). From the definitions of K and F(z) it
is concluded that S(z) is a proper transfer matrix. Once H(z) is proper then the product
S(z)H(z) is also proper. Consequently the sum R(z)+S(z)H(z) is biproper since
R(z) is, according to Theorem 4.1, biproper. Thus G,(z) is biproper. Bearing in mind
that Gy,(z) = H(z)G,(z) then it is concluded that G,,(z) is proper. Q.E.D.
4.3 MFD of Transfer Function Matrices
A transfer matrix may be represented as the "ratio" of two polynomial matrices as
G(z) = NR(z)D,(z) = D (z)N,(z)
where both matrices have the same number of columns. For an mxp matrix G(z) then
D(z) is m x m and N(z) is p x m. The first representation is a right matrix fraction
description, MFD, while the second is called a left MFD of G(z). Since duality exists
between these descriptions, for the most part, only right MFDs will be mentioned hereon.
The subscript R shall be omitted and restored only when necessary.
Example 4.1 Consider
1 0.5
G(z) = z z
0.04 0.05
which can be decomposed in a right MFD as:
which can be decomposed in a right MFD as:
N(z +15 0.51 zz+2.5 0
0.04z 0.05 2z z
An MFD is not unique and an infinity of others can be obtained by choosing any non
singular polynomial W(z) such that
N(z)= N(z)W (z) D(z)= D(z)W (z)
are polynomial matrices. Then
G(z)= N(z)D'(z)= N(z)D'(z)
since
N(z)= N(z)W(z) D(z)= D(z)W(z)
Matrix W(z) is right divisor of N(z) and D(z).
Some polynomial matrices have inverses that are also polynomial. These are
called unimodular and possess a determinant that is independent of z (a nonzero
constant). When all right divisors of two matrices are unimodular these matrices are
called right coprime and the MFD is called irreducible. Analogous to the scalar case,
coprimeness of two polynomial matrices means that no cancellation of common terms is
possible.
DEFINITION 4.1 A right MFD G(z)= N(z)D (z) is called irreducible if N(z) and
D(z) are right coprime. 0
No direct method for the generation of irreducible MFDs is known. The MFD
should be tested for coprimeness and, if necessary, both matrices should be "divided" by
their greatest common divisor. Bezout Identity can be used to determine if two matrices
are coprime as N(z) and D(z) will be right coprime if and only if there exist polynomial
matrices X(z) and Y(z) such that
X(z)N(z)+ Y(z)D(z) = I
The construction of a greatest common divisor is done by finding a unimodular
matrix U(z) such that at least p of the bottom rows on the right hand side are identically
zero:
U11(z) U,12z) D(z) C (z)
U21(Z) Uzz(Z) N(z)]
LU,(z) U22(z)]_LN(z)] L 0
Proof of this is done by Kaylath.
Example 4.2 Consider the MFD of Example 4.1. Matrices D(z) and N(z) are not
coprime. Their greatest common divisor is
1 0.5
C(z)=
0C=1 z+2.5
Make N = NC1 and D = DC1 then
z+1.5 0.57
N(z) = N(z)C (z)= 0.04z 0.021
0.04z 0.02_
z(z+2.5) 0.5zO
D(z) = D(z)C'(z)= z(z 2.5) 05z
which are now coprime. Note that this new MFD still represents G(z) correctly.
1 0.5
G(z z+1.5 0.5][z(z+2.5) 0.5z z z
G(z)= N(z)D'(z)= '= m
L 0.04z 0.021L 2z 0 0.04 0.05
4.4 Poles and Zeros of a Transfer Matrix
Poles and zeros of transfer matrices can be defined and determined in several
ways one of which consists in decomposing the transfer matrix into a canonical form
known as the SmithMcMillan form. Every rational transfer matrix G(z) can be reduced
to its SmithMcMillan form M(z) as
M(z) = diag ,..., 0, ..,
V (Z) 2 (Z) r (Z)
where e,(z) are called the invariant factors of M(z). Each pair {e,(z), W (z)} is coprime.
Note that multivariable systems can have poles and zeros at the same location once e,
and W,, for i # j, do not have to be coprime. Multiplication by unimodular matrices
does not affect the basic properties of transfer matrices. Therefore matrices that differ
only by unimodular factor, called equivalent matrices, have the same SmithMcMillan
form.
Example 4.3 Consider a process H(z)
z+1.5 0.5
z(z 1) ( 1)
H(z) =
0.04 0.05
(z 1) ( 1)
The SmithMcMillan form of H(z) is given by
0
M(z) z(z 1) 2.
z+2.5
0
z1
from where we conclude that the poles and zeros are p ={0,1,1} and z {2.5}. The
relative position of poles and zeros is seen along the diagonal of M(z). Thus it is best to
rewrite the SmithMcMillan form as
M(z) =diag (z) z 2.5
z(z1)' z1
The SmithMcMillan form tells us more than we need to know. For the analysis
of stability the relative position of poles and zeros, unless cancellations occur, is
irrelevant. Ultimately it is the presence of unstable poles in any of the closedloop
transfer matrices that will determine instability. The SmithMcMillan form is also very
elaborate computationally and not practical for real systems.
Another way of determining the poles of a system is by direct verification of the
poles of the individual transfer functions in the transfer matrix. This, however, requires
the calculation of the inverse of R+ SH which is, in general, a computationally difficult
task. The round off errors associated to the many different convolutions, required during
the process of inversion of polynomial matrices, makes the cancellation of common terms
unfeasible. Even for systems of modest sizes, e.g. systems with two inputs and two
outputs, the necessary cancellations, in general, do not occur. The use of Matlab
commands, such as minreal, which reduces a realization to a minimal realization, solved
the problem partially but cannot be considered a reliable procedure. Thus, from a
practical point of view, a solution to the design problem that would avoid or reduce the
number of inversion of matrices is preferred.
Matrix fraction descriptions may be used, instead, to determine poles of
multivariable systems and are more suitable for computational purposes.
PROPERTY 4.1 Let C denote the complex plane. Let G(z) be a proper transfer
function matrix and N(z), D(z) be an irreducible right MFD of G(z) such that
G(z)= N(z)D'(z)
where D(z) is nonsingular. Then p e C is a pole of G(z) if and only ifp is a root of the
determinant of D(z). 0
For square transfer matrices the zeros can also be determined from the MFD of
G(z) as they are the roots of the determinant of N(z). For nonsquare systems, however,
the concept of zeros is unclear.
The closedloop transfer matrices can be represented as MFDs. Let H = BA' be
a right MFD of H. Then G = (R + SBA')' and G, = A(RA + SB)'. Let
Q = RA + SB (4.5)
then
G,, = AK' (4.6)
the same substitutions in G, give Gy, = BA(R + SBA')' and Gy = B(RA + SB).
Finally
Gy, = B 1' (4.7)
Thus in order to be able to determine the poles of the closedloop system the right
MFDs of G, and G, need to be reduced to lowest terms. For that, let C,, be the
greatest right common divisor of A and Q,, i.e. C,, = gcrd(A, 0), then A = AC, and
K = QC,,. Substituting in equation (4.6) we get
G,, = A, (4.8)
Similarly let C, = gcrd(B, Q), then B = B C,, Q = .Cyw and
GyM = B D (4.9)
We now have the necessary conditions to determine the poles of the closedloop
transfer matrices.
THEOREM 4.3 Consider the closedloop system of Figure 4.1 with H(z) proper.
Thenpy is a pole of G,(z) if and only ifpy is a root of the determinant of 0.(z) and p,
is a pole of G ,(z) if and only if pu is a root of the determinant of .,(z). If H(z) is
square and nonsingular then Zy is a zero of G,,(z) if and only if Zy is a root of the
determinant of B,,(z) and z, is a zero of G (z) if and only if z, is a root of the
determinant of A, (z).
Proof Proof follows from Property 4.1. Q.E.D.
Note that the determination of the poles G,(z) or G ,(z) does not require the inversion
of 4 (z) or Q(z). Instead we only need to calculate their determinants.
Example 4.4 Consider a process H(z)
z+1.5 0.5
z(z 1) ( 1)
H(z)= =
0.04 0.05
(z 1) ( 1)
with irreducible right MFD H(z) = B(z)A1 (z)
Fz+1.5 0.5 z(z 1) 0
0.04z 0.05 0 (z 1)
let the predictive controller designed for this plant be such that Q(z) is
(z)=I z(z + 2.5) 07
2z zJ
The roots of the determinants of B(z), A(z) and Q(z) are {2.5}, {0, 1, 1} and {2.5, 0,
0,} respectively. Thus G,(z) is unstable once its zeros, {0, 1, 1}, cannot cancel the
unstable pole of Q(z) located at 2.5. Refer to examples 1 and 2 to see that matrices
B(z) and Q(z) are not coprime.
Matrix 2(z1) is always polynomial when expressed in backward shiftoperators.
In forward shiftoperators, however, it will normally be a rational matrix where the
denominators are powers of z. Because the poles of a system are defined in forward
shiftoperators, this conversion must be executed before the poles are calculated.
Multiplication by a diagonal matrix with appropriate powers of z will do this
transformation but the numerator matrix must be considered also so that no poles at the
origin are lost. Although the backward and forward shiftoperator versions of a transfer
function are mathematically equivalent changing the arguments may be misleading.
Therefore the symbol G(z) will be used to refer to the forward shiftoperator version of
G(z1).
DEFINITION 4.2 The conjugate of a polynomial M(z1) = mo + mz +... + mz'l in
backward shiftoperator is defined as MA(z) =M(z')zl =mz nAz" + mzz' +...+mn .
The degree of the polynomial in backward shiftoperator is the degree of its conjugate
polynomial. 0
DEFINITION 4.3 The conjugate of a polynomial matrix M(z1) in backward shift
operator is defined as M(z) = M(z')zM where znM =diag{znM ,znM',...,znM} is a
diagonal matrix and z"nM is defined as the maximum degree of all polynomials along the
jth column of M(z1) )
Note that for matrices with polynomials of different degree along its columns this
procedure places excess zs in the lower degree polynomials. This is not the case of
Q2(z') where degree of all polynomials along its columns are the same. This can be seen
from the definitions of R, A, S and B. Thus K2(z1) possesses a conjugate matrix
Q(z) = Q(z')z"n where no excess zs are present.
When transforming backward shiftoperators MFDs into forward shiftoperators
MFDs the respective z"M of numerator and denominator matrices must be reduced to
lowest terms to guarantee coprimeness.
DEFINITION 4.4 The degree of a transfer function
B(z')_ b,+b,z +...+bznB
A(z') a +az1 +...+aznz
is the max{nB,nA}. 0
DEFINITION 4.5 The conjugate of a right MFD H(z) =B(z')A'(z') is
H(z)=B(z)A'(z) where A(z)=A(z')z'H and B(z)=B(z')znH. Matrix z"H is a
diagonal matrix defined as z'" =diag{zn"H "H ...,znH"' where nH, is the maximum
degree of all transfer functions along thejth column of H(z) )
From Definition 4.5 it is trivial to see that B(z')A (z) = B(z)A(z)
In the case of the closedloop transfer matrices, where numerator and denominator
matrices do not have the same degree, a coprime factorization of the individual z"
matrices is necessary. Thus for G,,(z') we have that
G, (z') = A(z')K(z) = A(z)znH (z(z)znK)1 = A(z)z nHIzn1(z)
Multiply matrices z"""z to get Z mzn where the common terms have been
cancelled. Then
G(z') = A(z')i'(z')= A(z)z' (L(z)H)1
G,(z')= A(z' )z"nHz ((2z'1) nH )1
define A(z)= A(z )z"zmn and Q(z)= z(z')znz i
Gu(z)= A(z)Q(z)1
and similarly
G,(z1) = B(z1)1(z) = (z)z)nna ((z)zf)1
G, (z)= B(z')z"zn n(a(z )nznH)1
define B(z)= B(z )z"zn
Gy. (z)= h(z)21(z)
and the determinants are given by:
det(2(z)) = det(2(z')) det(znl) det(zl )
det(B(z)) = det(B(z')) det(zn) det(,ln )
det(A(z)) = det(A(z1)) det(zn") det(,ln )
Note that for nonsquare systems the determinant of det(Q(z)) does not exist. For the
analysis of stability and performance the above development shows that the backward
shiftoperator polynomial matrices contain all the relevant information.
DEFINITION 4.6 Let N(z')D (z1) be a coprime MFD of G(z1). The finite roots
of det(D1 (z')) are the relevant poles of G(z1). 0
Example 4.5 Consider the process of Example 4.4
1 '+ z21.5 0.5z1
H(z) = z1 lz
0.04z1 0.05z1
1 zz 1 z1
36
A^1) z1 0 z [ 0]
0 1z1 0 :
z1 + Z21.5 0.5z1
S 0.04z1 0.05z1
1+)2.5z1 0 2
2z' 1 _0 1 _
then
A(z) =A(z 1)zH Iz1 0 z2 ] = [Z(z) o)]
0 z' 0 z 0 (z1)
(z=B(z) '+z21.5 0.5z1z 2 0Fz+1.5 0.5
(z 0.04z1 0.05z 0 z 0.04z 0.05
Q(1) (Z1)Z 1+2.5z1 0 0 z+2.5 0
2z' 1 0 1 2 1
1 0
z_
mi 1 0 1 I 0
Z I o =[ ( 11 0 0 Z(Z 1) 2
Z) =(Z)Zmln 10.4z 1 0 .04z 0.05
det(a(z'))det(z"a) det(z) =(1+ 2.5z)(z)(z2) = (z + 2.5)z2
det(B(z)) det(z") det(z) = (0.03z2 + 0.075z3)(z3)(1)= 0.03z + 0.075
det(A(z')) det(z ) det(z ) = (1 z)2 (z3)(1) = (z 1)2 z
if we use the backward shiftoperator to obtain the determinants
det(Q(z1)) =1+ 2.5z 2.5
0.03z + 0.075
det(B(z')) = 0.03z +0.075z 3
z
det(A(z')) = (1 ')2 (= z
The particular format of the closedloop transfer matrices and the possibility of
determination of the poles of the system lead to the following theorems:
THEOREM 4.4 Let H be a proper transfer function. The closedloop system will be
stable if the roots of det( Q,~) and det( Q, ) are within the unit circle.
Proof From Theorem 4.3 the poles of transfer matrices G, and G, are given by the
roots of det(,) and det(9,). Thus if these roots are within the unit circle then the
system is stable. Q.E.D.
However if A and B are coprime then these matrices cannot cancel the same
unstable poles of G, and G,.
THEOREM 4.5 Let H be a proper transfer function. Let H = BA1 be a coprime
matrix fraction description of H. The closedloop system will be stable if and only if the
roots of det( Q) are within the unit circle.
Proof The closedloop transfer matrices G,, and G, have, before any cancellation of
common terms, the same denominator matrix 2. Therefore, if the roots of the det(Q)
are within the unit circle, then ,W and , have their roots within the unit circle and
according to Theorem 4.3 the system is stable. Conversely if any one of the roots of the
det( ) are outside the unit circle they cannot all be cancelled by A and by B once these
are coprime matrices and either ,W or or both, will have roots outside the unit
circle and according to Theorem 4.4 the system is unstable Q.E.D.
Although there is no direct method for obtaining a coprime matrix fraction
representation of the process transfer matrix it is, in general, a simple procedure if the
transfer matrix of the plant is minimal. Stability can be determined by evaluating the
position of the roots of det( Q) which can be done with Jury's criteria.
Theorem 4.5 is very important for the implementation of a design procedure for
multivariable predictive control once it allows a fast initial test of stability of the system.
4.5 Performance Analysis
Crisalle et al. (1989) have proven that for the SISO case if Nu=Ny and 2 =0 then
the predictive controller generates a control that results in the inversion of the plant and
the output matches the reference. A similar result is obtained for the MIMO case as can
be seen next.
THEOREM 4.6 Consider a system with equal number of inputs and outputs, i.e. m = p.
Let A = 0 and Nu, = Ny, for all j, then the resulting closedloop transfer matrix is the
identity, Gy, A H(R + SH)1T I.
Proof Recall equation (3.10), Au = (G' iG + A)'G'c(r y0), in the derivation of the
control law. The condition that A=0 and Nu,=Ny, for all j implies that
(GTOG + A)'GT, = G1. Matrix G is a block lower triangular matrix and its inverse
has a similar structure. Therefore, recalling that K is the matrix of rows corresponding to
the first control move of (G T(G + A)'GT (which is G1 in this case), its elements are
vectors with a nonzero first element followed by zeros, in this case. This particular
format of K has a similar effect of restricting the prediction horizons to Ny = 1 in the
formulation of the control problem, once the elements of K (defined as
k =[k, k2 ky ],j) are k =[k1 0 ... 0],,. Thus KP=KP and KF K,1F
(see Appendix B for details of both proofs). From the definitions of R and S then
S = K,F,
Therefore Q becomes
S= A(I+ z'KlP)AR+ z1KlFBR (4.10)
From Matrix Diophantine equation, for i= 1 and bearing in mind that E1 =I, then
F = z(IAA,) and from equation (3.16) we get PJ= z(BL G,). Substituting in the
expression for Q then
0 = A[I+z'Kz(BL G,)] AR +zKz(I AA,)BR (4.11)
but K, = G1 and
S= AAR + AG'BAR AAR+ G'BR G'AABR .
The coprime factorizations of H imply that BLAR = ALBR. Cancellation of common
terms leads to
sK = G'B,
(4.12)
This way
G, = BR(G'B )G = I
Example 4.6 Consider the process of Example 4.4.
z+1.5 0.5
H z(z1) (Z 1)
0.04 0.05
(z1) ( 1)
A multivariable predictive controller was designed for
following parameters were used: Ny = 1, Nu = 1, 1 = 1, 1, = 1, h
R, S and Tare:
z +1.5z 2.5 5(2z 1) 50(2z 1)
R= z S 3z 3z
R 2(zl) z1 4(2z1) 100(2z1)
Z2 3z 3z
(RA + SB) (z + 2.5)
SB(A2 004
GYW = B(RA + SB) = z
0.04
Z
Q.E.D.
this process and the
= 0, = 0. Matrices
5 50
z z
T= _3 3
4 100
3 3
0
0.5
0.05
Z
from where we see that
Gyr = B(RA + SB)1T=o
It is of interest to establish conditions under which the multivariable predictive
controller produces offsetfree responses to step changes in the set point of any given
output as was proven by Crisalle et al. (1989) for SISO systems.
THEOREM 4.7 Consider a plant H with equal number of inputs and outputs where (i)
H(1) is nonsingular, (ii) the weights (, in the objective function equation are > 0 for all
j and (iii) the resulting closedloop is stable. Then it achieves zero offset response
lim, [r(t) y(t)]= 0
to step changes in the reference r(t).
Proof In the z domain r(z) y(z)= [I G, (z)] r(z). From the Final Value Theorem
lim_ [r(t) y(t)]= lim, A[ Gy(z)] (4.13)
where is a vector of step inputs of amplitude a The transfer matrix G, (z) is
A
Gy, = H(R+SH)' T
From the definitions of R and S (3.19) to (3.21) Gy, becomes
Gy, = H(A(I + zK, )+KH)T,
From Matrix Diophantine equation F = z (I E AA) and
Gy, = H[A(I + z'K,)+Kz'(IAA)H] T
Recall that Kz' = T At z = 1 the delta operator becomes A = 0 and therefore
and substituting in (4.13)(1) = H()(T()H())T()=I
and substituting in (4.13)
Q.E.D.
lim, [r(t) y(t)]= 0
Condition (i) requires the plant to have a nonsingular gain, i.e., there exits the
inverse matrix H'(1). This condition is a necessary requirement for offsetfree
behavior, otherwise, there is no solution u to the steadystate equation r = H(1)u where
r is the final value of the stepfunction in the set point and u is the final value of the
input required to deliver an output that exactly matches the final value of the reference.
Systems with singular gains are functionally uncontrollable (Skogestad, 1996) and hence
perfect tracking for any arbitrary step change is impossible to attain. Finally, condition
(iii) is necessary because ( = 0 for some value of j, then deviations from the set point
for output y, are not penalized in the objective function equation, and consequently, the
controller cannot ensure attainment of perfect tracking.
4.6 Implementation
The mathematical development of this chapter allowed the generation of a Matlab
code for the analysis of stability of unconstrained multivariable predictive control.
Based on Theorem 4.5 a design procedure was created where the critical step is
the generation of a right coprime description of a process model, H, for which no direct
method is known. Maple was used to generate the greatest common divisors for systems
where reduction of MFDs to lowest terms was necessary. In order to avoid the roundoff
errors of floating point operations a rational representation of the coefficients of the
polynomials was done allowing the division algorithm for polynomials to give exact
results. This way the elementary row and column operations, required in the process of
obtaining the greatest common divisors, are executed correctly. The drawback of this
approach resides in the fact that as the convolutions of the polynomials are executed the
rational number representations of the coefficients tend to have increasing numerators
and denominators. This imposes a limitation on the size of MFDs that can be reduced to
lowest terms.
For square systems the 'numerator' matrix possesses (in general) a determinant
and a simple (initial) test of coprimeness is possible once coprime matrices will not have
common roots in their determinant. The presence of common roots, however, does not
imply reducibility of the MFD once poles and zeros can be located in different positions
in the transfer matrix.
Thus given matrices R and S (refer to Chapter 3) and a right coprime MFD of
the process matrix (H=BA') then Q is generated, Jury's criteria is applied to the
determinant of Q to check for stability. In the presence of unstable roots then these
parameters are discarded for another set until only stable roots are present.
Performance is verified by simulating the system in Simulink. The block diagram
structure of Simulink allows the simulation of the closedloop system without the need to
invert R+ SH as would be the case if a numerical integrator were to be used. Although
the inversion of R is still necessary the following algebraic manipulation manages to use
the block diagram properties of Simulink to execute the inversion:
R = (I + z'KP)A = I z'(I AKP)
Let R' = z' (I AKP) then
R=IR'
and R' = (I R')'. This way the inversion of matrix R can be represented in block
diagram as seen in Figure 4.2.
R1(z) = +
R '(z)
Figure 4.2 Inversion of R in Simulink
The Matlab code was modified to generate matrix R'. Simulations done showed very
good results.
4.7 Conclusion
In this chapter the closedloop behavior of a plant with a multivariable
generalized predictive controller is analyzed and a method, based on matrix fraction
description (MFD), is developed to check for stability and determine the closedloop
poles for any choice of tuning parameters. Performance is analyzed by simulating the
closedloop system in Simulink.
Properties such as "zero offset response to step changes in the reference" and
"inversion of the plant" are proven for plants with equal number of inputs and outputs.
CHAPTER 5
THREEPHASE SEPARATOR
Offshore production systems are responsible for the treatment of the petroleum
produced in the sea fields. These systems are subject to intense fluctuations in its feed
due to the multiphase nature of the fluids coming from the reservoir. The feed consists of
a mixture of oil, water and gas with different degrees of dispersion between phases.
Separators are responsible for absorbing fluctuations of the feed as well as for promoting
the separation. Its performance is of vital importance for the quality of both, the oil
exported to refineries and the water discharged to the sea. Economic and environmental
restrictions have led the industry to research means of improving the efficiency of the
equipment. Little has been done, though, in the use of advanced control techniques.
Separators possess three independent SISO PI (proportional and integral) control
loops, one for each phase. Frequently, bad tuning of the controllers is responsible for a
low yield in the separation process. The apparent periodicity of the multiphase flow of
the feed raises hope that it can be predicted at some degree. Although this phenomenon
will not be studied here, future values of set points (reference trajectory) are assumed to
be known.
Therefore this dissertation will, hopefully, show the benefits and lay the formal
grounds for the future development of predictive controllers for three phase separators.
Most separators are horizontal as will be considered in this dissertation.
5.1 Description of the Separator
Separators use gravity as the driving force for separation. For that reason small
droplets tend not to separate and a low efficiency is normally observed.
\
WATER j U1L
Figure 5.1 OilWaterGas Separator
The vessel is divided in two sections separated by a weir; a separation chamber
and an oil chamber. Gravitational force promotes the segregation of phases and a water
rich phase with dispersed oil settles at the bottom of the separation chamber while an oil
rich phase with dispersed water lays atop. An interface between the two phases is
formed. The liquids flow in the direction of the weir and along this path a series of
parallel plates help the liquidliquid separation. The oil rich phase flows over the weir
into the oil chamber. The water rich phase is discharged just before the weir to the water
treatment unit.
An interface level controller manipulates the opening of the water (outlet) valve.
A level controller controls the level of oil in the oil chamber. Pressure inside the vessel is
regulated by a controller acting on the gas valve.
The outputs (measurable states) are pressure, height of water (height of interface)
and height of oil in the oil chamber. The inputs are the opening of the valves.
Parallel plates are used to reduce the length of space required for the separation of
dispersed droplets. These adhere to the plates and coalesce. The plates are inclined
allowing for the fall (in the case of water droplets) or rise (in the case of oil) droplets.
Figure 5.2 Parallel Plates
5.2 Rigorous Model
The model here developed is a simplification of a rigorous model (Nunes, 1994)
which has the following aspects:
The feed is composed of 15 chemical components and one pseudocomponent.
The model possesses a total of 81 variables.
Thermodynamic equilibrium of phases. This is assumed because liquid phase
residence time is approximately 5 minutes while thermodynamic equilibrium
takes approximately 30 seconds to occur. SoaveRedlichKwong (SRK) equation
of state is used.
Droplet distribution of dispersions at the entrance of the separator is known.
Stokes regime is considered for the calculation of terminal velocity of droplets of
water and oil. This is a reasonable assumption since the liquid phases have low
velocities and laminar flow is observed.
Droplets are perfect spheres.
Effects of wall and concentration are negligible in the calculation of terminal
velocities
No coalescence of droplets occurs.
No emulsifying agents present. This implies that the system is a dispersion.
No dispersion of liquid particles in the gas phase.
A parabolic profile for the velocity of water is considered in the formulation.
Droplets are uniformly distributed in space at the entrance of parallel plates.
PARALLEL PLATES
0 o
0 0 U u 0
Figure 5.3 Velocity Profile
5.3 Simplified Model
Simulations done with the rigorous model show that thermodynamic
representation of phases is important in studying the behavior of the gas but has little
effect on the dynamics of liquidliquid separation. Because the main interest here is to
control and improve the oilwater separation, the following assumption is made:
No mass transfer between thermodynamic phases occurs.
The major impact of this assumption is that variations in pressure, temperature or
composition do not cause a mass change of the individual thermodynamicc) phases.
Because fluctuations observed (in practice) are negligible, the assumption is justified and
allows for a simpler model with 7 state variables.
Thermodynamic calculations determine the total amount of each phase present in
the vessel. The degree of dispersion between phases, however, is calculated using
dynamic population balances. To avoid confusion between thermodynamic phases and
dispersions, notation is carefully defined.
VG 'Gas'
VL' 'Oil' ,
VW 'Water' VL
Figure 5.4 Volumes
Following common practice, 'water' will, hereon, refer to the water rich phase
lying at the bottom of the separation chamber as seen in Figure 5.4. Note that it is a
mixture of two thermodynamic phases, namely, the oil particles (dispersed phase) and the
water (continuous phase). Properties associated to the 'water' will have a subscript W
while lower case italic superscripts will refer to the thermodynamic phases; w standing
for the water and f for the oil. Thus V' is the volume of oil in the 'water' and V, is the
volume of water in the 'water'. The total volume of the 'water' is:
vw Wv17+v ;w
Following the same rationale 'oil' will be used to name the mixture of oil
(continuous phase) and water (dispersed phase) present in the separation chamber and the
oil chamber. Distinction, however, must be made between the two since they have
different properties such as volume, concentration of water, etc. Properties associated to
the 'oil' in the oil chamber will have a subscript L while to represent the 'oil' in the
separation chamber a subscript Lsc will be used. Thus VL is the volume of water in the
'oil' in the oil chamber.
Subscript G will be used for the gas phase. Note that the assumption of no
dispersion of droplets of oil or water in the gas phase makes it unnecessary to distinguish
continuous and dispersed phases.
The above definitions allow the following algebraic relations to be established for
volumes and concentrations:
Total volume restriction V = VG + V + VLsc + VL
Volume of 'water' in separation chamber: VW = VW + V
Volume of 'oil' in separation chamber Vcs = VSs + VW
Volume of 'oil' in oil chamber V = L + V
Concentration of dispersed water in 'oil' at separation chamber Cc VSc
VLsc
Concentration of dispersed oil in 'water' at separation chamber C = 
Vw
w V
Concentration of dispersed water in 'oil' at oil chamber C =
VL
5.4 Flow Equations
Figure 5.5 shows the heights and flow rates of each phase. The equations are seen
next:
Figure 5.5 Heights and Flow rates
Qoutw =
CvsW Vd (p pw)+YwhW + YLsc(hT h)
Qweir = 24.88s 2g lengthhw, 0.2(h hw,,,))(hT hw,,)1 5
QoutL
QOUtG
(5.1)
(5.2)
(5.3)
CvLsL dLppL )+yLk
pLO.0693
sVGsGdG(ppG)(p+p)
PG2.832
(5.4)
5.5 Differential Equations
Total height at separation chamber (derived from the total material balance):
dh QinL + Qinw Qweir Qoutw
dt 2 C (D h)h
Height of water (derived from the mass balance of the water phase):
(5.5)
p,0.0693
dh Qin (1 TOG ) QinBSW E(
(5.6)
dt 2C,, (D hw)h
Height of oil at the oil chamber (derived from the mass balance of the oil phase):
dhL Qweir Qut (5
dt 2C; (D h )h
Balance of volume of water dispersed in the 'oil' in the separation chamber:
dV" V"w
Ls QinLBSW(1 E) Qweir LWc (5.8)
dt VLsc
Balance of volume of oil dispersed in the 'water' in the separation chamber
SV= QinTOG(1 E )Qoutw (5.9)
dt VW
Balance of water dispersed in the 'oil' in the oil chamber
dVL = QweirL, s QoutL Vf (5.10)
dt VL VL
Pressure equation (derived from the material balance of gas)
dp 1 (RT)2 (Qin QOutG)
 = p   ("GL + "W OutL Outw) (5.11)
dt VV,VLc VL pMWG
5.6 Conclusion
In this chapter a simplified dynamic model of the separator was created based on
a more rigorous model previously developed. This simplification greatly reduced the
number of states without loss of accuracy for the variables of interest in the development
of a controller for the equipment. This model was programmed in Matlab and
Simulations done in Simulink were used to check it against the rigorous model giving
good results.
hT
hLsc
hL
hw
D
C
C s
C'
Lsc
dL
dw
dG
w
ELsc
eL
QinL
Qinw
QinG
Qweir
5.7 Nomenclature
total height ('water' and 'oil') at separation chamber
height of 'oil' in separation chamber
height of 'oil' in the oil chamber
height of 'water' in separation chamber
diameter of separator
length of separation chamber
length of oil chamber
concentration of water in the 'oil' in the sep. chamber
concentration of oil in the 'water' in the sep. chamber
concentration of water in the 'oil' in the sep. chamber
density of 'oil' in oil chamber
density of 'water' in oil chamber
density of 'gas' in oil chamber
efficiency of removal of water from 'oil' in oil chamber
efficiency of removal of oil from 'water' in oil chamber
efficiency of removal of water from 'oil' in sep. chamber
flow rate of 'oil' into separator
flow rate of 'water' into separator
flow rate of 'gas' into separator
flow rate of 'oil' over weir
[m]
[m]
[m]
[m]
[m]
[m]
[m]
[m3/m3]
[m3/m3]
[m3/m3]
[kgf/cm3]
[kgf/cm3]
[kgf/cm3]
dimensionlesss]
dimensionlesss]
dimensionlesss]
[m3/minutes]
[m3/minutes]
[m3/minutes]
[m3/minutes]
QoutL
Qoutw
QOUtG
V
VG
VW
VLs
VL
VLs
vt
VW
P
SG
T
BSW
TOG
flow rate of 'oil' out of separator
flow rate of 'water' out of separator
flow rate of gas out of separator
Total volume of separator
Volume of the gas
Volume of the 'water'
Volume of the 'oil' in separation chamber
Volume of the 'oil' in oil chamber
Volume of oil in the 'oil' in separation chamber
Volume of dispersed water in 'oil' in separation chamber
Volume of oil in 'oil' in oil chamber
Volume of dispersed water in 'oil' in oil chamber
Volume of oil in the 'water' in oil chamber
Volume of water in the 'water' in oil chamber
Pressure
Constant of gases
opening of the 'oil' valve
opening of the 'water' valve
opening of the 'gas' valve
Temperature
concentration of water in the 'oil' in the feed
concentration of oil in the 'water' in the feed
[m3/minutes]
[m3/minutes]
[m3/minutes]
[m]
[m3]
[m3]
[m3]
[m3]
[m3]
[m3]
[m3]
[m3]
[m3]
[m3]
[kgf/cm3]
[atm.liter/mol.K]
dimensionlesss]
dimensionlesss]
dimensionlesss]
[K]
[m3/m3]
[m3/m3]
CHAPTER 6
DESIGN OF A PREDICTIVE CONTROLLER FOR THE SEPARATOR
In this chapter a multivariable predictive controller is designed for the separator
with the help of the Matlab code previously developed. A linearized model of the
separator is used and the closedloop system is simulated for the analysis of performance.
Many different combination of tuning parameters are tested. For the set of
parameters where stability is verified, performance is analyzed by simulating the closed
loop system. Initially two sets of parameters are chosen, a first one, which results in very
intense oscillations in the height of water and a second set where the parameters are
successfully changed with the purpose of reducing these oscillations. Further on another
set up where the pressure is controlled separately and the two remaining controlled
variables (heights of water and oil) are controlled by a multivariable predictive controller
is proposed and studied.
6.1 Discrete Linear Model of the Separator
To design a predictive controller for the separator a discrete linear model of the
equipment is necessary. The nonlinear model of the separator, developed in Chapter 5,
was identified by evaluating the open loop response of the system to step perturbations in
inputs (opening of valves). The following model was generated where the variables are
now deviation variables.
17 0 2.4
h 206s+1 367s+1 W
L s,
126 169 43
322s+1 330s+1 508s+1
0 0 2.4 sG
2s+l1
This model was discretized with a sampling time of 0.01 minutes or 0.6 seconds to give:
0.000825z1 0.0000650z1
0
1z' 1z1
0.00391z 1 0.00512z 1 0.000846z1
S z1 z1 z1
0 0 0.01197z' sG
0 0
10.995z1
The following list shows the relationship between the variables of the objective
function equation and the variables of the model:
yl output 1................... h height of water
y2 output 2...................... h height of oil
y3 output 3...................... p pressure
ul control action 1................ s opening of water valve
u2 control action 2............... s opening of oil valve
U3 control action 3.............. sp opening of gas valve
The linearized model implies the use of deviation variables (null initial
conditions). Consequently the changes in the reference are deviations around the original
steady state of the nonlinear model. The trajectory of the reference (or set point) is seen
in Table 6.1:
Table 6.1 Variations in magnitude and time of change of set point
Time (minutes) Initial Value Final Value
hw 25 0 0.1
h 50 0 0.1
P 75 0 0.1
6.2 Design of Predictive Controller
Following the development of Chapter 4 a left MFD and a right irreducible MFD
are created for the model.
0.000825z1
BL(z')= 0.00391z1
0
1z1
AL(z1)= 0
0
0.000825z1
BR(z1)= 0.00391z1
0
0 0.0000650z1
0.00512z1 0.000846z1
0 0.0120z'
0
10
0
0
0
10.995z'
0
0.00512z1
0
0.0000651z1
0.000842z'
0.0120z'
1z1 0 0.000396
AR(z1) = 0 1z1 0.000525
0 0 (1 0.995z1)
The determinant of BR has no relevant roots while the determinant of AR has the
following set of relevant roots: {1,1,0.995}. Recall that because the system has same
number of inputs as outputs these roots correspond to the zeros of Gy, and G
respectively.
6.3 Case 1: First Tuning
To tune the controller many different combinations of tuning parameters were
tested and most of them generated unstable systems. Table 6.2 lists a set of parameters
for which stability was verified.
Table 6.2 Tuning parameters for Case 1
hw P
Ny 20 20 20
Nu 8 8 8
2 10 8 5
5 3 7
Matrices R and S are:
1 z1
R(z) = 0
0
0
1z'
0
0
0
1z
1.2 +1.118z 1
S(z')= 0.10730.1004z 1
0.05104 0.04703z'
3.004 +2.797z '
4.953 + 4.612 z'
0.3880.3574 z1
0.1537+0.1433 z'
0.2454 + 0.2288z1
13.98+12.86 z'
and the resulting matrix Q is
K(z') = R(z')A, (z')+ S(z')B, (z') =
11.987z' +0.9882z2 0.01538z 0.01432z2
0.01923z 0.01790z2 11.975z 1 +0.9764z2
0.001555z'1 +0.001433z2 0.001987z'1 +0.001830z2
0.0003962+0.001164z'1 0.0007128z2
0.0005247+0.001752z 1 0.001138z2
1+1.827z 0.8408z2
The determinant of K2(z') is
det(Q(z1))=
z6 +5.7893z5 13.9692z4 +17.9827z 13.0261z2 +5.0343z 0.81101
6
z
whose relevant roots are
p, ={0.9137 0.0772i,0.9813 0.0471i,0.9996 + 0.0069i}.
All roots lie within the unitary circle and the system is stable. Matrix Tis seen next:
T(1,1)= 0.0004 0.0008z0.0012z2 0.0016z3 0.002z4 0.0024z 0.0028z6 0.003 17 0.0035z8
0.00399 0.0043z10 0.0047z1 0.005 l12 0.0055z13 0.0058z14 0.0062z1
0.0066z16 0.007z7 0.0074z18 0.0077z19
T(2,1) = 0.000007 + 0.00002z + 0.00004z2 + 0.00007z3 + 0.0001z4 + 0.00014z + 0.00018z6 + 0.00022z7 + 0.00027z8
+0.0003z9 +0.00036z'0 +0.0004z1 +0.00045z2 +0.000493 + 0.0005414 +0.00058215 + 0.0006z16
+ 0.0007Z7 + 0.0007z18 + 0.00076219
T(3,1)= 0.000054 + 0.0001z + 0.00014z2 + 0.00016z + 0.00018z4 + 0.0002z5 + 0.0002z6 + 0.00021z7 + 0.00021 z
+ 0.00022z9 + 0.00022z10 + 0.00022z1 + 0.00022 + 0.00023z3 + 0.00052314 + 0.00023z5 + 0.00024z16
+0.00024'7 + 0.00024z8 + 0.000252'9
T(1, 2) =0.0011 0.0022z 0.0033z 0.0043z 0.005z4 0.00635 0.007z6 0.0087 0.009z8
0.01z9 0.0lzo1 0.01211 0.01312 0.014z3 0.015214 0.01615
0.016Z16 0.017Z7 0.018Z8 0.019z19
T(2,2)= 0.0019 0.0036 0.00542 0.007z 0.0087z4 0.015 0.012z6 0.013z7 0.015z8
0.016z9 0.018z10 0.019z11 0.02Z1'2 0.013z' 0.013'14 0.014Z'5
0.0271'6 0.029Z'7 0.03z'8 0.0322'9
T(3,2)= 0.0004 +0.00077z + 0.001z2 +0.0013z3 + 0.00144 +0.0015z5 +0.0016z6 +0.0016z7 +0.0016z8
+0.001629 +0.0017z10 +0.0017z1 +0.0017212 +0.0017z1 +0.0018214 +0.001821
+0.001816 +0.001817 +0.00118Z +0.0019219
T(1, 3) = 0.000012 0.000035z 0.000067z2 0.0001 z1 0.00016z4 0.0002z 0.00028z6 0.00034z7 0.00041z8
0.00047z9 0.00054z10 0.0006z1 0.00067z12 0.00073z3 0.0008z14 0.0009z5 0.0009z16
O.OOl17 O.OOz18 O.OOlz19
T(2,3) = 0.00002 0.000056z 0.000 z2 0.00017z3 0.00025z4 0.00034z5 0.00044z6 0.00055z7 0.00065z8
0.0007629 0.00086z10 0.00096z11 0.00 12 0.0012z3 0.0013z14 0.0014z5
0.0015z16 0.001617 0.0017z18 0.00182'9
T(3,3)= 0.0140.026z 0.035z2 0.043z 0.04824 0.05325 0.05526 0.05727 0.058z8
0.05929 0.061z10 0.062z1 0.06412 0.06513 0.066214 0.068Z1
0.069z16 0.07Z17 0.072z18 0.073219
Figure 6.1 shows the results of the simulation of the linear model for the height of
water. Excessive oscillation occurs and this can be explained by the fact that a very small
sampling time was chosen in order to deal with the very fast dynamics of the pressure.
Besides this, during the design of the controller, the greatest weight was placed on the
first control action, Aul (opening of the valve of water), as X1 is 10 and less importance
was given to the error between the reference (set point) and the output (water level), rlyl,
as 01 is 5. This causes a restriction on the opening of the water valve while the level is
free to reach higher values of error.
0.2
0.18
0.16
0.14 
0.12 
hw 0.1
0.08 
0.06
0.04
0.02
0
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.1 Height of Water
The multivariable aspect of the controller can be observed by the apparent
absence of perturbations in the level of water when the height of oil changes (set point
change at 50 minutes) and when the pressure changes (set point change at 75 minutes).
Both these changes will promote an increase in the pressure of the separator, which in
turn would increase the outflow of water. Consequently a decrease in the water level
would be seen. Due to the predictive and multivariable nature of the predictive control
the water valve (see Figure 6.4) closes while the oil valve is closing thus avoiding the
decrease in the level. The end result is a smooth curve of height of water.
0.14
0.12
0.1 
0.08
hL 0.06
0.04
0.02
0
0.02
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.2 Height of Oil
The height of oil oscillates much less. If we recognize that fluctuations in the
height of water directly affect the flow of oil over the weir (which is the feed to the oil
chamber) then we see that most of the oscillations were caused by the behavior of the
water level.
The peak seen at 50 minutes is caused by the change in the set point of level of
oil. The level rapidly reduces and if it were not for the oscillations in the height of water
a final value would soon have been reached. Thus we may conclude that the control of
oil level had a good performance.
The results for the pressure show an excellent tracking of the set point as seen in
Figure 6.3. Because the change in heights of oil and water have very little effect on the
pressure this variable is almost independent of the others. However the opposite is not
true, the pressure has a very strong effect on the heights since the flow rates leaving the
vessel are a direct function of the internal pressure of the separator.
0.12
0.1
0.08
0.06
0 10 20 30 40 50 60 70 80
time(min)
90 100
Figure 6.3 Pressure
The independence of the pressure added to the multivariable and predictive
aspects of the controller create the expectation of a good performance for the control of
pressure. The opening of the valves are seen in Figures 6.4, 6.5 and 6.6.
0.8
0.6
0.4
SW
0.2
0.2
0.4
0.6
0.8
0 10 20 30 40 50 60 70 80
time(min)
Figure 6.4 Opening of the Water Valve
0 10 20 30 40 50 60 70 80
time(min)
Figure 6.5 Opening of the Oil Valve
90 100
90 100
0.05
0
0.05
0.1
0.15
Sp 0.2
0.25
0.3
0.35
0.4 ....
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.6 Opening of the Gas Valve
Thus the control of the oil level and pressure gave good results. The control of
water level, however, did not show a good performance and another set of parameters is
searched next. Looking for a better performance of the predictive controller the weights
are changed. The aim of this new tuning is to reduce the oscillations of the level of water
while maintaining the performance of the other two variables.
6.4 Case 1: Second Tuning
A second set of parameters was tested. Again, a series of different tuning
parameters were tested and from the set of parameters that resulted in a stable closedloop
systems one was chosen. Table 6.3 shows the set of parameters used in this case.
Table 6.3 Tuning parameters for the second tuning of Case 1
hwkh P
Ny 20 20 20
Nu 9 9 9
A 3 5 3
5 1 2
The weight on the control action of the water valve was reduced (from 10 to 3).
The weights on the other control actions were also reduced. The weight on the error
between the level of water and the set point, 4i, was maintained. Thus the relative
importance of this term in the objective function equation is higher and an improvement
in the performance of the control of water level is expected while for the other variables,
oil level and pressure, it should deteriorate. Control actions should be more intense.
Matrices R and Tare:
1z 1 0 0
R(z') = 0 1z1 0
0 0 1 z
3.914 +3.646z' 3.444 + 3.207z' 0.1445 +0.1348z1
S(z')= 0.20470.1913z1 2.757+2.567z' 0.1058+0.09866z1
0.12450.1152z1 0.34210.3168 z' 10.25+9.470 z
and matrix Q is giving by
K(z') = R(z')A, (z)+ S(z1
11.983z' +0.9845z2
= 0.01058z' 0.09854z2
0.0014372' +0.00133022
)BR(z1
0.01764z' 0.01642z2
11.986z' +0.9869z2
0.001752z' +0.001622z2
0.003962+0.001821zl 0.001324z2
0.0005247+0.001567zl 0.0009685z
1+1.8720zl 0.8814z2
de 1) z6 + 5.841z5 14.22z4 +18.47z 13.5022 +5.266z0.8561
det(G(z'))= 6
whose relevant roots are
p={0.9360 0.0725i,0.9855 0.0423i,0.9992+ 0.0106i}.
This lets us conclude that the system is stable. Matrix Tis seen next:
T(1,1)= 0.001354 0.00269z 0.00422 0.0053z3 0.0065824 0.00785z5 0.0091z6 0.010427 0.0116z8
0.012829 0.014z10 0.015z1 0.0165z12 0.0178z13 0.02214 0.02z15
0.0215z16 0.0227z17 0.024z18 0.025z19
T(2,1)= 0.0000137 + 0.00004z + 0.000078z2 + 0.00013z3 + 0.0001924 + 0.000255 + 0.00033z6 + 0.0004z7 + 0.0005z8
+ 0.00058z9 + 0.00067z10 + 0.00076z11 + 0.00085z12 + 0.00094z13 + 0.001214 + 0.001 12'5 + 0.0012z16
+0.0013z17 +0.0014z18 +0.0015z19
T(3, 1) = 0.000095 + 0.00018z + 0.00025z2 + 0.0003z + 0.00035z4 + 0.00039z5 + 0.00042z6 + 0.00045z7 + 0.00047z8
+ 0.00049z9 + 0.0005z10 + 0.0005z" + 0.0005z12 + 0.00056z1 + 0.00058z14 + 0.0006z5 + 0.00062z16
+ 0.0006z7 + 0.00066z18 + 0.00067z19
T(1, 2) = 0.0013 0.0025z 0.00372 0.0049z3 0.006z4 0.0071z5 0.0082z6 0.009z7 0.01z8
0.01129 0.0125zo 0.014z1 0.015Z2 0.0157z3 0.017214 0.018Z5
0.019Z16 0.02z17 0.02118 0.022219
T(2, 2) =0.001 0.002z 0.0032 0.004z3 0.004824 0.005725 0.006526 0.007427 0.00828
0.009z9 0.01z 0.01z 0.012'2 0.013z' 0.013z14 0.014z1
0.015z16 0.01617 0.017z18 0.018219
T(3,2) = 0.00025+ 0.00046z + 0.00064z2 + 0.0008z3 + 0.000924 + 0.00125 + 0.00126 + 0.001227 + 0.0013z8
+0.0013z9 +0.0014z0 +0.001z +0.0015z2 +0.0015z3 +0.001614 +0.0017z5
+ 0.0017z16 + 0.0018z17 + 0.0018z18 + 0.0019z19
T(1, 3) =0.00001 0.00003z 0.000062 0.0001z 0.000144 0.00019z5 0.00025z6 0.0003 l7 0.00037z8
0.00044z9 0.0005z10 0.0005611 0.00063z12 0.00069z1 0.00075z14 0.0008z5 0.00088z16
0.00094Z17 0.OOlzi 0.00l19
T(2,3)= 0.000008 0.000023z 0.000044z2 0.00007z 0.00014 0.00014z5 0.00018z6 0.0002z7 0.00027z8
0.00032z9 0.00037z10 0.0004z1 0.0004612 0.0005z13 0.00055214 0.00059z'
0.00064z16 0.00069z17 0.00073z8 0.0008z'9
T(3,3) =0.0071 0.013z 0.0192 0.023z 0.027z4 0.03 lz 0.034z6 0.036z7 0.038z8
0.04z9 0.042z10 0.044z11 0.04612 0.048z1 0.05z14 0.052z1
0.054Z16 0.05617 0.058z18 0.06z19
By comparison to Figure 6.1 it is clear how this tuning causes less oscillations in
the height of water. The amplitudes are smaller and the attenuation is faster.
0.18
0.16
0.14
0.12 
0.08 
0.06
0.04
0.02
0
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.7 Height of Water
Figure 6.8 shows the results for the height of oil. If we compare it to Figure 6.2
we see that the oscillations have increased. Bearing in mind that with the improvement
of the control of water level the oscillations in the flow rate of oil to the oil chamber
(Qwe,,,) will decrease, then it is easy to conclude that the performance of the control of oil
level has worsened. This is an expected downside to the change made in the parameters.
For the pressure this new tuning has had little effect. The results are very close to
the last tuning (see Figure 6.9). Once again this is caused by the independence of the
pressure with relationship to the other variables.
0.14
0.12
0.1
0.08
0.06
hL 0.04
0.02
0
0.02 
0.04
0.06 I
0 10 20 30 40 50 60 70 80
time(min)
90 100
Figure 6.8 Height of Oil
0 10 20 30 40 50 60 70 80 90 100
time (min)
Figure 6.9 Pressure
1
0.5
0
SW 0.5
1
1.5 ...
0 10 20 30 40 50 60 70 80
time(min)
Figure 6.10 Opening of Water Valve
1
0.8
0.6
0.4
0.2
SL 0
0
0.2
0.4
0.6
0 10 20 30 40 50 60 70 80
time(min)
Figure 6.11 Opening of Oil Valve
90 100
90 100
The control action for the opening of the water valve (see Figure 6.10) is more
aggressive while the opposite has happened to the oil valve and the gas valve. This
behavior was also expected from the modification in parameters.
0.05
0
0.05
0.1
0.15
0.2
Sp
0.25
0.3
0.35
0.4
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.12 Opening of the valve of Gas
It was very difficult to tune the 3x3 predictive control for the separator. The
above tuning does not show a good performance. The reason for this lies in the fact that
the transfer function relating the opening of the gas valve to the pressure has a time
constant very small (T=2) if compared to the other transfer functions. That way a very
small sampling time was required (Ts=0.01minutes), which, for the horizons used
(Ny=20 and Nu=9), caused the loss of the predictiveness of the predictive control once 20
sampling times correspond to 1.2 seconds. The controller became very sensitive to the
tuning parameters and only a few sets of these generated stable closedloop systems.
6.5 Case 2
In this case the control of pressure was left out of the loop and a predictive
controller was designed for a system with two inputs (openings of water and oil valve)
and two outputs (height of water and oil).
17
hw 206s+1 sw
h 126 169 sL
322s+1 330s+1
The absence of the pressure in the transfer matrix allowed the sampling time to be
much higher (0.1 minutes). With this new sampling time the model was discretized
yielding:
0.00825z'
LL7 0.0391z1 0.0512z1 sL
1z' 1 z1
A sensitivity analysis was done to measure the influence of the prediction horizon
on the location of the closedloop poles. Table 6.4 shows the control horizon and weights
adopted.
Table 6.4 Tuning parameters for Case 2
hw_______ h
Nu 9 9
2 1 10
10 5
The prediction horizon and the relevant roots of det(2(z1)) are seen on Table
6.5. It is seen that as the prediction horizon increases the location of the roots of
det(2(z1)) move (slightly) towards the origin. Because the change in the position of the
roots of det(2(z')) are small with prediction horizon, the value of Ny=20 was chosen
for the design of the controller.
Table 6.5 Variations of relevant roots of det(2(z1)) with prediction horizon
Ny Roots of det(Q(z'))
20 0.99860.0138i 0.93620.0726i
30 0.99560.0199i 0.92230.0541i
40 0.99060.0243i 0.92050.0451i
50 0.98390.0260i 0.92040.0432i
60 0.97650.0245i 0.92030.0446i
Matrices R, S and Q are
Sz 112.84+11.91z'
1.4881.393z1
11.888z' +0.8966z2
0.01322z 0.01225z2
25.82+23.92z'1
3.693 + 3.425z'
0.1322z' 0.1225z2
11.981z' +0.9825z2
whose determinant is
det((z )) = z 3.8695z +5.61842 +3.6283z+0.8794
4
and the relevant roots of det((z )) are
and the relevant roots of det(f2(z1)) are
p, = {0.9362 0.0726i,0.9986 0.0138i}
^ [O 1z]
Matrix Tis seen next:
T(1,1)= 0.007410.01405z 0.01998z2 0.02526z 0.02997z4 0.03416z 0.03789z6 0.04122z7 0.0442z8
0.04688z9 0.04956z1' 0.05224z1 0.05492z12 0.0576z3 0.06028z14 0.06296z5 0.06564z16
0.06832z17 0.07z18 0.07368z19
T(2,1)= 0.0000995 + 0.00029z + 0.000566z2 + 0.000916z3 + 0.001336z4 + 0.001817z5 + 0.002352z6 + 0.002934z7 + 0.003558z8
0.004217z9 +0.004876z10 +0.005535z" +0.00619512 + 0.006854z3 + 0.007513z14 + 0.008172z1 + 0.00883 lz16
0.009491z17 +0.01015z8 +0.01081z19
T(1, 2) =0.01725 0.03238z 0.04558Z2 0.05702z3 0.06687z4 0.07531z5 0.08249z6 0.08857z7 0.09369z8
0.0989 0.1023z' 0.1066Z2 0.1109293 0.115214 0.11955 0.1238z16 0.128127 0.1325'8 0.1368219
T(2,2) 0.00228 0.004302z 0.0060882 0.00766z 0.009039z4 0.01024z 0.0113z6 0.0122z7 0.01302z8
0.01372z9 0.01443z1' 0.01513z1 0.0158312 0.01653z3 0.01724z14 0.01794z1 0.01864z16
0.0193517 0.02005z18 0.02075z19
For the PID control of pressure the tuning parameters were:
Table 6.6 Tuning Parameters of PID
Proportional 0.5
Integral 0.1
Derivative 0
The results for this case are very good. The height of water doesn't show the
oscillatory behavior seen in Case 1. From Figure 6.13 the predictive nature of the
controller is clearly seen. The height of water starts rising at 10 sampling times, or 1
minute, before the actual change in the set point occurs. A smooth transition to the final
value follows and a very small overshoot is seen.
0.12
0.1
0.08
hw0.06 
0.04
0.02 
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.13 Height of Water
The peak seen at 50 minutes is due to the change in height of oil. The increase in
oil level elevates the pressure of the vessel and the water valve (Figure 6.16) closes
preventively to counteract this higher pressure.
0.1
0.08
0.06
hL
0.04
0.02
0
0.02 ....
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.14 Height of Oil
A PID controls the pressure and the response is seen in Figure 6.15.
0.18
0.16
0.14
0.12
0.1
P
0.08
0.06
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.15 Pressure
0.1
0
0.1
0.2
SW
0.3
0.4
0.5
0.6
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.16 Opening of Valve of Water
0.4
0.35
0.3
0.25
0.2
L0.15
0.1
0.05
0
0.05
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.17 Opening of the Valve of Oil
0.1
0.05
0.05
0.15
0.2
0.25
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 6.18 Opening of the Valve of Gas
6.6 Conclusion
The design of a predictive controller with three inputs and three outputs for the
separator was difficult due to the difference in time constants between the pressure and
the levels. Initially a small sampling time was required and the search for a stable system
with good performance required many trials. A second case was studied where a
predictive controller with two inputs and two outputs was designed for the control of
levels while a PID controller was used for the pressure. A sensitivity study was done
where the closedloop poles are calculated for different prediction horizons. Tuning the
controller was much easier in this case and a good performance was observed.
CHAPTER 7
SIMULATION AND RESULTS
In this chapter both predictive controllers designed in Chapter 6 are tested against
the nonlinear model of the separator by simulating the closedloop in Simulink of Matlab.
The same series of step changes of set point made in Chapter 6 (see Table 1) is used.
7.1 Results of Case 1
Figure 7.1 shows the change in the height of water. The results resemble the
simulations done for the linearized model in Chapter 6 (see Figure 6.1). The fluctuations
decay very fast and a steady state is reached before the set point of oil height is changed.
0.66
0.64
0.62
0.6
0.58
hw
0.56 
0.54 
0.52 
0.5
0.48
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.1 Height of Water
90 100
At 50 minutes a peak occurs due to the change in height of oil. If we observe the
opening of the water valve in Figure 7.4 we see that the controller initially closes the
valve to counteract the increase in pressure caused by the increase in oil level. This is
done before the change in set point of oil level demonstrating the predictive and
multivariable characteristics of the controller. While this peak is almost unnoticed for the
linear model it is very evident in this case as it is a consequence of the nonlinearities of
the model.
Figure 7.2 shows the response of the oil level. Compared to Figure 6.8 we see
that, here, the decay is slightly faster. Once again the nonlinearities are responsible for
this change in response. In general terms the performance is good.
0.85
0.8 ,
hL
0.75
0.7
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 7.2 Height of Oil
The response of the pressure is very similar to that of the linear model.
9.62
9.6
9.58
9.56
9.54
P
9.52
9.5
9.48
9.46
9.44
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.3 Pressure
1.2
1 
0.8
0.6 l
SW
0.4
0.2 
0
0.2
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.4 Opening of Water Valve
90 100
90 100
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.5 Opening of Oil Valve
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.6 Opening of Gas Valve
1.6
1.4
1.2
1
SL
0.8
0.6
0.4
90 100
0.85
0.8
0.75
Sp
0.7
0.65
90 100
7.2 Results of Case 2
Figure 7.2 shows the response of the water level. Compared to the linear model
(Figure 6.13) the overshoot is smaller while the peak around 50 minutes is very similar.
0.62
0.6 
0.58
0.56
hw
0.54
0.52
0.5
0.48
0 10 20 30 40 50 60 70 80 90 100
time(min)
Figure 7.7 Height of Water
All other simulations show that the results for the nonlinear model are very close
to that of the linear model. It is important to note that the height of the weir (0.9 meters)
is the maximum allowable value for the height of water and that the change in set point
was 0.1 (from 0.5 to 0.6 meters). Thus the new set point corresponds to a change in 20%
in the height of water, which makes it a good test on how the linear model applies to the
separator.
Figures 7.8 through 7.12 show the response of the nonlinear model to the
remaining variables.
0.84
0.82
0.8
0.78
0.76
hO.74
0.72
0.7
0.68
0.66
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.8 Height of Oil
9.65
9.6
9.55 
P
9.5
9.45
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.9 Pressure
90 100
90 100
0.8
0.75
0.7
0.65
0.6
SW
0.55
0.5
0.45
0.4
0.35
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.10 Opening of Water Valve
1
0.95
0.9
0.85
SL 0.8
0.75
0.7
0.65
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.11 Opening of Oil Valve
90 100
90 100
0.72
07
0.68
0.66
0.62 ...
0 10 20 30 40 50 60 70 80
time(min)
Figure 7.12 Opening of Gas Valve
90 100
7.3 Conclusion
The nonlinear simulations gave results close to the linear model of Chapter 6.
Thus we conclude that the linearization was successful and the design of the predictive
controller for the separator using the techniques of linear control is a good approach.
CHAPTER 8
CONCLUSIONS
This dissertation has successfully developed a method based on matrix fraction
descriptions (MFD) for determining the closedloop poles of unconstrained multivariable
predictive control.
The method, proposed here, takes advantage of the polynomial nature of the
matrices of the controller to represent the closedloop transfer matrices as an MFD, a
ratio of "numerator" and "denominator" polynomial matrices. It is shown that if a right
coprime matrix description of the process transfer matrix is done then the closedloop
poles of the system are determined by evaluating the roots of the determinant of the
denominator matrix giving necessary and sufficient conditions for stability. This method
avoids the inversion of transfer matrices which is a numerically difficult task. Instead
only the determinant of the denominator matrix is required and stability is verified with
the help of Jury's criteria by determining if there are any poles that lie outside the unitary
circle. Because in this formulation of predictive control a transfer function model was
used to represent the process, as opposed to overparametrized models such as
step/impulse response models of DMC like methods, the degree of the determinant is
considered minimal. This approach reduces enormously the numerical complexity of the
problem allowing the tuning of multivariable systems with many inputs and outputs.
It is proven that the system has zero offset response to step changes in the
reference, a property known to be valid for the single input single output case. For
systems with equal number of inputs and outputs the "inversion of the plant" is also
proven. In this case it is seen that if the weights on the input are zero then the solution to
the optimization problem is a controller that inverts the plant and the output matches the
reference.
This method becomes an important tool for the analysis of predictive controllers
allowing practitioners and researchers to follow the effect of tuning parameters on the
behavior of the system.
The use of a multivariable predictive controller for an oilwatergas separator is
studied. The nonlinear model of the plant is linearized around a steady state and a three
inputs threeoutputs predictive controller is designed for it. It is seen that the controller
responds positively to changes in the parameters and performance objectives can be
pursued. Another setup is also studied where a controller with two inputs and two
outputs is designed to control the levels of water and oil. The pressure is controlled
separately. Results show agreement with the simulations done for the linear model and it
is concluded that predictive control is a successful control method for oilwatergas
separators.
APPENDIX A
DETAILED DERIVATION OF CONTROL LAW EQUATION
In this appendix intermediate equations used in the derivation of the control law
equation are developed, starting from the model of the plant. Multiple input single output
(MISO) formulation is initially used to show the basic structure of the equations. MIMO
system equations are derived as a combination of these.
A. 1 Derivation of Matrix Diophantine Equation
Consider the Darma model
A, O ... 0 y (t) B, :2 B"' u (t1)
0 422 .. 0 y2(t) B21 B22 B2 p 2 t1)
0 0 Am ym (t)_ Bml m2 Bmp up (t 1)
Analyze, initially, a MISO system of one output, y1(t), and p inputs, u (t),
u2 (t) ... ,u(t). For that the first row of the above matrices can be represented as
follows:
u, (t 1)
u 2(t1)
A4,y (t)=[B 12 .. B] u (A.1)
u,(t 1)
where
4, (q') = I + a q1 + a2,,q +... + a,,, q
