Title: Design and analysis of multivariable predictive control applied to an oil-water-gas separator
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00100821/00001
 Material Information
Title: Design and analysis of multivariable predictive control applied to an oil-water-gas separator a polynomial approach
Physical Description: Book
Language: English
Creator: Nunes, Giovani Cavalcanti, 1961-
Publisher: University of Florida
Place of Publication: Gainesville Fla
Gainesville, Fla
Publication Date: 2001
Copyright Date: 2001
 Subjects
Subject: Predictive control   ( lcsh )
Stability   ( lcsh )
Chemical process control   ( lcsh )
Chemical Engineering thesis, Ph. D   ( lcsh )
Dissertations, Academic -- Chemical Engineering -- UF   ( lcsh )
Genre: government publication (state, provincial, terriorial, dependent)   ( marcgt )
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )
 Notes
Summary: ABSTRACT: This dissertation uses a polynomial-operator technique to study stability and performance of unconstrained multivariable predictive control. A simple and direct way to determine stability of the closed loop system is developed. It is shown that to guarantee stability of the closed loop 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 closed loop. 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 single-output 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.
Summary: ABSTRACT (cont.): The use of a multivariable predictive control for an oil-water-gas 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 oil-water-gas 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.
Summary: KEYWORDS: predictive control, dynamic systems, multivariable control, optimization, optimal control
Thesis: Thesis (Ph. D.)--University of Florida, 2001.
Bibliography: Includes bibliographical references (p. 115-117).
System Details: System requirements: World Wide Web browser and PDF reader.
System Details: Mode of access: World Wide Web.
Statement of Responsibility: by Giovani Cavalcanti Nunes.
General Note: Title from first page of PDF file.
General Note: Document formatted into pages; contains viii, 118 p.; also contains graphics.
General Note: Vita.
 Record Information
Bibliographic ID: UF00100821
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: oclc - 49261928
alephbibnum - 002763563
notis - ANP1585

Downloads

This item has the following downloads:

Phd ( PDF )


Full Text










DESIGN AND ANALYSIS OF MULTIVARIABLE PREDICTIVE CONTROL
APPLIED TO AN OIL-WATER-GAS 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 CLOSED-LOOP BEHAVIOR ................................................... 23
4.1 Stability Analysis of Closed-loop Multivariable Systems ........................................23
4.2 Closed-loop 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 THREE-PHASE 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. CLOSED-LOOP TRANSFER MATRICES.................... .......... ............. 110
D 1 Disturbances in the Closed-loop........................................................... .............. 110
D.2 Closed-loop 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 OIL-WATER-GAS SEPARATOR: A POLYNOMIAL APPROACH

By

Giovani Cavalcanti Nunes

August 2001

Chairman: Tim Anderson
Major Department: Chemical Engineering

This dissertation uses a polynomial-operator technique to study stability and

performance of unconstrained multivariable predictive control. A simple and direct way

to determine stability of the closed-loop system is developed.

It is shown that to guarantee stability of the closed-loop 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 closed-loop.

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 single-output 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 oil-water-gas 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 oil-water-gas 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 high-order 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 finite-impulse response models and finite-step

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)







t-3 t-2 t-1 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 single-output single-input process is

Ny2 Nu
J= (i)(r(t+i) -(t+i It))2 + i)Au(t+i-1)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 non-minimum phase plants. It uses transfer function model

for the plant leading to lower-order representations.

Different transfer function models are used by GPC such as Controller Auto-

Regressive Integrated Moving Average (CARIMA), Controller Auto-Regressive Moving

Average (CARMA), etc. In this dissertation a Deterministic Auto-Regressive 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 input-output pairs can be considered as a single-input single-output (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 multiple-inputs and

multiple-outputs (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

open-loop 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 mono-variable 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

non-minimum 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 Predictor-Based Self-Tuning

Control. Extended Prediction Self-Adaptive Control (EPSAC) by De Keyser et al.

(1985) proposes a constant control signal starting from the present moment while using a

sub-optimal predictor instead of solving a Diophantine equation. Later on the input was

replaced by the increment in the control signal to guarantee a zero steady-state 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. State-space 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 dual-mode

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 end-point

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

state-space 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 on-line solution to

open-loop 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 Hamilton-Jacobi-Bellman (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+i-1) (3.1)
=1 =1

Nu2 Nup
+1 2 (i)AU2 (t +i-1)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 vector-matrix format.

This will simplify the notation and help visualization of the solutions.

Claim 3.1. Equation (3.1) can be written in compact vector-matrix 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 G-AuT 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(t-1). The predicted values of the output can be generated

using Matrix Diophantine equation,

E (q-')AA(q-')+ q-' (q-) = I, (A.9)

where E (q-1) and 1F(q-1) 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(t-1) (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 + i-1)+ 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+i-1)=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+ q-yK 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 shift-operator while

matrix T(q) is in forward shift-operator.

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 CLOSED-LOOP BEHAVIOR

In this chapter the closed-loop system is analyzed for stability and performance.

For that a practical method to determine the closed-loop 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 Closed-loop Multivariable Systems

The issue of ensuring stability of a closed-loop 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

closed-loop is still not available for multivariable systems. Although unconstrained

systems present analytical solutions to the problem of minimizing the objective function,

obtaining the closed-loop 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 closed-loop 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 closed-loop 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 closed-loop system, in the

sequel, the complex variable z is substituted for the forward shift operator q. This

substitution is allowed by the analogy between shift-operator 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 Closed-loop Transfer Function Matrices

Figure 4.1 shows the block diagram of the multivariable predictive control applied

to a plant H(z).


Figure 4.1 Closed-loop System









Let H(z) be the process plant. Using block-diagram algebra operations, the

dynamics of the system shown in Figure 4.1 can be captured by the closed-loop 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 = NC-1 and D = DC-1 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 Smith-McMillan form. Every rational transfer matrix G(z) can be reduced

to its Smith-McMillan 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 Smith-McMillan

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 Smith-McMillan form of H(z) is given by


0
M(z) z(z 1) 2.
z+2.5
0
z-1

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 Smith-McMillan form as


M(z) =diag (z-) z 2.5
z(z-1)' z-1









The Smith-McMillan 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 closed-loop

transfer matrices that will determine instability. The Smith-McMillan 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 non-square systems, however,

the concept of zeros is unclear.

The closed-loop 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 closed-loop 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 closed-loop

transfer matrices.









THEOREM 4.3 Consider the closed-loop 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)A-1 (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(z-1) is always polynomial when expressed in backward shift-operators.

In forward shift-operators, 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

shift-operators, 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 shift-operator 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 shift-operator version of

G(z-1).

DEFINITION 4.2 The conjugate of a polynomial M(z-1) = mo + mz- +... + mz-'l in

backward shift-operator is defined as MA(z) =M(z-')zl =mz nAz" + mzz-' +...+mn .

The degree of the polynomial in backward shift-operator is the degree of its conjugate

polynomial. 0

DEFINITION 4.3 The conjugate of a polynomial matrix M(z-1) 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(z-1) )

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(z-1) possesses a conjugate matrix

Q(z) = Q(z-')z"n where no excess zs are present.

When transforming backward shift-operators MFDs into forward shift-operators

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- +...+bz-nB
A(z-') a +az-1 +...+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 closed-loop 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)z-nH (z(z)z-nK)-1 = A(z)z- nHIzn-1(z)

Multiply matrices z"""z to get Z mz-n 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,(z-1) = B(z-1)-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)2-1(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(z-1)) det(zn") det(,ln )

Note that for non-square systems the determinant of det(Q(z)) does not exist. For the

analysis of stability and performance the above development shows that the backward

shift-operator polynomial matrices contain all the relevant information.

DEFINITION 4.6 Let N(z-')D- (z-1) be a coprime MFD of G(z-1). The finite roots

of det(D-1 (z-')) are the relevant poles of G(z-1). 0

Example 4.5 Consider the process of Example 4.4

1 -'+ z-21.5 0.5z-1

H(z-) = -z-1 l-z-
0.04z-1 0.05z-1
1 z-z 1 z-1






36

A^1-) -z-1 0 z [ 0]
0 1-z-1 0 :

z-1 + Z-21.5 0.5z-1
S 0.04z-1 0.05z-1

-1+)2.5z1 0 -2
-2z-' 1 _0 1 _

then

A(z) =A(z 1)zH I-z-1 0 z2 ] = [Z(z-) o)]
0 --z-' 0 z 0 (z-1)

(z=B(z-) '+z-21.5 0.5z-1z 2 0Fz+1.5 0.5
(z 0.04z-1 0.05z- 0 z 0.04z 0.05

Q(1) (Z-1)Z 1+2.5z-1 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.03z-2 + 0.075z-3)(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 shift-operator to obtain the determinants


det(Q(z-1)) =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 closed-loop 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 closed-loop 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 = BA-1 be a coprime

matrix fraction description of H. The closed-loop system will be stable if and only if the

roots of det( Q) are within the unit circle.









Proof The closed-loop 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 closed-loop 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, = G-1. 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 G-1 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+ z-1KlFBR (4.10)

From Matrix Diophantine equation, for i= 1 and bearing in mind that E1 =I, then

F = z(I-AA,) 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 +z-Kz(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(z-1) (Z 1)
0.04 0.05
(z-1) ( 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(z-l) z-1 4(2z-1) 100(2z-1)
Z2 3z 3z


(RA + SB) (z + 2.5)
SB(A-2 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 offset-free 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 closed-loop 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 + z-K, )+KH)-T,

From Matrix Diophantine equation F = z (I E AA) and

Gy, = H[A(I + z-'K,)+Kz'(I-AA)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 offset-free

behavior, otherwise, there is no solution u to the steady-state equation r = H(1)u where

r is the final value of the step-function 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 round-off

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 closed-loop 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=I-R'

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 closed-loop 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 closed-loop

poles for any choice of tuning parameters. Performance is analyzed by simulating the

closed-loop 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
THREE-PHASE 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 Oil-Water-Gas 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 liquid-liquid 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 pseudo-component.

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. Soave-Redlich-Kwong (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 liquid-liquid separation. Because the main interest here is to

control and improve the oil-water 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 dLp-pL )+yLk
pLO.0693

sVGsGdG(p-pG)(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 V-V,-VL-c- VL p-MWG




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 closed-loop 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.000825z-1 0.0000650z-1
0
1-z-' 1-z-1
0.00391z 1 0.00512z 1 0.000846z-1
S- z-1 -z-1 -z-1

0 0 0.01197z' sG
0 0
1-0.995z-1

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.00391z-1
0


1-z-1
AL(z-1)= 0
0

-0.000825z-1
BR(z-1)= -0.00391z-1
0


0 0.0000650z-1
-0.00512z-1 0.000846z-1
0 -0.0120z-'


0

1-0
0


0
0
1-0.995z-'


0
-0.00512z-1
0


-0.0000651z-1
-0.000842z-'
-0.0120z-'


1-z-1 0 -0.000396
AR(z-1) = 0 1-z-1 -0.000525
0 0 -(1- 0.995z-1)

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- z-1
R(z-) = 0
0


0
1-z-'
0


0
0
1-z


-1.2 +1.118z 1
S(z-')= 0.1073-0.1004z 1
0.05104 0.04703z-'


-3.004 +2.797z '
-4.953 + 4.612 z-'
0.388-0.3574 z-1


-0.1537+0.1433 z-'
0.2454 + 0.2288z-1
-13.98+12.86 z-'













and the resulting matrix Q is


K(z-') = R(z-')A, (z-')+ S(z-')B, (z-') =

1-1.987z' +0.9882z-2 0.01538z --0.01432z-2
0.01923z- -0.01790z-2 1-1.975z 1 +0.9764z-2
-0.001555z'1 +0.001433z-2 -0.001987z'1 +0.001830z-2


-0.0003962+0.001164z'1 -0.0007128z-2
-0.0005247+0.001752z 1 -0.001138z-2
-1+1.827z- -0.8408z-2


The determinant of K2(z-') is


det(Q(z-1))=


-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.0008z-0.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.014-0.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), rl-yl,

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 closed-loop

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:

1-z 1 0 0
R(z-') = 0 1-z-1 0
0 0 1- z-

-3.914 +3.646z-' -3.444 + 3.207z-' -0.1445 +0.1348z-1
S(z-')= 0.2047-0.1913z-1 -2.757+2.567z-' -0.1058+0.09866z-1
0.1245-0.1152z-1 0.3421-0.3168 z-' -10.25+9.470 z-












and matrix Q is giving by

K(z-') = R(z-')A, (z-)+ S(z-1
1-1.983z-' +0.9845z-2
= 0.01058z-' -0.09854z-2
-0.0014372-' +0.0013302-2


)BR(z-1
0.01764z-' -0.01642z-2
1-1.986z-' +0.9869z-2
-0.001752z-' +0.001622z-2


-0.003962+0.001821z-l -0.001324z-2
-0.0005247+0.001567z-l -0.0009685z
-1+1.8720z-l -0.8814z-2


de 1) -z6 + 5.841z5 -14.22z4 +18.47z -13.5022 +5.266z-0.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 closed-loop 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.0512z-1 sL
1-z-' 1- z-1

A sensitivity analysis was done to measure the influence of the prediction horizon

on the location of the closed-loop 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(z-1)) are seen on Table

6.5. It is seen that as the prediction horizon increases the location of the roots of

det(2(z-1)) 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(z-1)) 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.488-1.393z-1

1-1.888z-' +0.8966z-2
0.01322z- -0.01225z-2


-25.82+23.92z'1
-3.693 + 3.425z-'

0.1322z-' -0.1225z2
1-1.981z-' +0.9825z-2


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(z-1)) are


p, = {0.9362 0.0726i,0.9986 0.0138i}


^ [O 1-z-]












Matrix Tis seen next:


T(1,1)= -0.00741-0.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 closed-loop 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 closed-loop 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 closed-loop 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 closed-loop 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 closed-loop

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 over-parametrized 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 oil-water-gas separator is

studied. The nonlinear model of the plant is linearized around a steady state and a three-

inputs three-outputs 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 oil-water-gas

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 (t-1)
0 422 .. 0 y2(t) B21 B22 B2 p 2 t-1)

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(t-1)
A4,y (t)=[B 12 .. B] u (A.1)

u,(t -1)

where

4, (q-') = I + a q1 + a2,,q +... + a,,, q




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs