Development of Adaptive Model

Predictive Control Toolbox

Submitted in the partial fulfillment of the requirements for
the degree of

Bachelor of Technology (Hons.) and Master of Technology

by

Vivek Kumar Pandey

Roll no.09D02032
Supervisor:

Prof. Sachin Patwardhan

Department of Chemical Engineering

INDIAN INSTITUTE OF TECHNOLOGY, BOMBAY
2013-2014
i
INDIAN INSTITUTE OF TECHNOLOGY, BOMBAY
Department of Chemical Engineering

Acceptance Certificate

The report for the “Dual Degree Project’’ entitled “Development of Adaptive Model
Predictive Control Toolbox” submitted by Mr.Vivek Kumar Pandey (Roll No-09d02032)
may be accepted for evaluation.

Dept. of Chemical Engineering _____________________

IIT Bombay (Prof. Sachin Patwardhan)

14th July 2014

ii
DECLARATION

I declare that I truly and honestly abide by the honour code that I have signed. I also declare
that I have properly cited the facts and works done by others and that I have not
misinterpreted or falsified any ideas in my submission. I understand the rules and regulations
of IIT Bombay and that any violations of above will result into severe disciplinary action
against me by the Institute.

_______________________

(Signature)

Vivek Kumar Pandey

Roll No- 09D02032

iii
Acknowledgement

I wish to express my heartiest gratitude to Prof. Sachin Patwardhan whose invaluable

guidance and constant motivation has been of immense help in my work. I would also like to
thank him for taking out time to help me understand the important concepts and plan my
work accordingly.

Vivek Kumar Pandey

IIT Bombay

iv
Abstract
Model Predictive Control is one of the best methods for optimal control of dynamical
systems. Nowadays, Model Predictive Control algorithms are being increasingly used in
process and petrochemical industries around the world. In recent times, MPC algorithms are
being explored in other fields like automotive and aerospace industries as well. However,
depending upon the type of systems, industries and need, different industries may use
different formulations. The type of MPC algorithms also depend on the model of the system.
The major factors that are taken into account while implementing a particular type of MPC
algorithm is performance and stability. Since process industries are huge complexes, it
becomes important that controllers should not become unstable. That is to say, along with
performance, it also becomes necessary to ensure the nominal stability. This report attempts
to present a generic formulation of MPC which is nominally stable and which can be used to
derive simpler formulations as and when needed. The domain knowledge of the system can
be used to know which particular formulation will be applicable for that system. The report
presents overview of some of the popular algorithms of MPC. It also outlines how we can
optimize MPC model parameters in the real time and hence update the model continuously.
This is called Adaptive Model Predictive Control. This algorithm is very useful when model
of the plant changes with time.

The generic theoretical formulation of MPC presented in this report is over the Infinite
Horizon. The objective function that has to be optimized for MPC is added with Terminal
state term. The corresponding weighting matrix is found by solving the discrete Lyapunov
Equation.

The report also includes a case study of the Alstom Benchmark Challenge for Gasifier
Control. The report aims to solve the challenge problem partially using MPC. Control of the
linear and non-linear gasifier using the linear first principle model at some operating points
has also been attempted

Effort has also been made to model and control the system.

v
Contents
Acknowledgement iv

Abstract v

List of Figures and Tables vii

Nomenclature ix

1 Introduction 1
1.1 Background…………………………………………… 1
1.2 Motivation for MPC and AMPC……………………… 1
1.3 Introduction to MPC and AMPC……………………... 2
1.4 Literature on MPC and AMPC……………………….. 3
1.5 Objectives…………………………………………….. 4
1.6 Outline of the report………………………………….. 4

2 MPC Toolbox 5
2.1 Modeling and Estimation…………………………….. 5
2.2 Generic MPC Formulation…………………………… 7
2.3 QP Formulation………………………………………. 10
2.4 Approaches to MPM Estimation……………………… 15
2.5 Adaptive Control……………………………………… 18
2.6 Summary……………………………………………… 20

3 Alstom Benchmark Challenge 21

3.1 Introduction …………………………………. …….. 21
3.2 Control System Specifications……………………..… 22
3.3 Results and Discussions …………………………… 24
3.4 Summary ………………………………………… 46

4 Conclusions 47
4.1 Summary…………………………………………….. 47
4.2 Future Work …………………………………………. 47

vi
Appendix & References

List of Figures
1.1. Diagrammatic Representation of MPC Algorithm…………………… 2
2.1. Height in Tank………………………………………………………. 19.
2.2. Variation of Input Voltage …………………………………………. 20
3.1 Schematic Model of Gasifier…………………………………………. 22
3.2 CVGAS vs Time(change in Air flow)………………………..………. 25
3.3 PGAS vs Time ( change in Air Flow)………………….…….………... 25
3.4 TGAS vs Time (change in Air Flow)………………….………...……… 26
3.5 CVGAS vs Time (change in Coal flow)………………………..……….. 26
3.6 PGAS vs Time ( change in Coal Flow)………………….…….………... 27
3.7 TGAS vs Time (change in Coal Flow)………………….………...……. 27
3.8 CVGAS vs Time(change in Steam flow)………………………..……… 28
3.9 PGAS vs Time ( change in Steam Flow)………………….…….………... 28
3.10 TGAS vs Time (change in Steam Flow)………………….………...……. 29
3.11 Output 1 and Output 2 vs Time (Linear Plant 100)………………………..….. 31
3.12 Output 3 and Output 4 vs Time (Linear Plant 100) ……………………………32
3.13 Input 1 vs Time (Linear Plant 100) .………………….…….………........ 32
3.14 Input 2 and Input 3 vs Time (Linear Plant 100)………………….………...…. 33
3.15 Input 4 and Input 5 vs Time (Linear Plant 100)……………………………….. 33
3.16 Output 1 and Output 2 vs Time (Linear Plant 50)…………………………… 34
3.17 Output 3 and Output 4 vs Time (Linear Plant 50) …………………………… 35
3.18 Input 1 vs Time (Linear Plant 50).……………………………………… 35
3.19 Input 2 and Input 3 vs Time (Linear Plant 50)……………………………… 36
3.20 Input 4 and Input 5 vs Time (Linear Plant 50)………………..….. 36
3.21 Output 1 and Output 2 vs Time (Non-Linear Plant 50) ……………… 37
3.22 Output 3 and Output 4 vs Time ( Non-Linear Plant 50)………………….…... 38
3.23 Input 1 vs Time (Non-Linear Plant 50)………………….………..……. 38
3.24 Input 2 and Input 3 vs Time (Non-Linear Plant 50)…………………………. 39
3.25 Input 4 and Input 5 vs Time (Non-Linear Plant 50)…………………………… 39
3.26 Output 1 and Output 2 vs Time (Non-Linear Plant 100)……………………… 40
3.27 Output 3 and Output 4 vs Time(Non-Linear Plant 100)……………………… 41
3.28 Input 1 vs Time (Non-Linear Plant 100)……………………………………… 41
3.29 Input 2 and Input 3 vs Time (Non-Linear Plant 100)………… 42
3.30 Input 4 and Input 5 vs Time (Non-Linear Plant 100) ………………………. 42
3.31 Output 1 and Output 2 vs Time (Setpoint Change 100)……………… 43
3.32 Output 3 and Output 4 vs Time (Setpoint Change 100)…………………… 44

vii
3.33 Input 1 vs Time (Setpoint Change 100)……………. 44
3.34 Input 2 and Input 3 vs Time (Setpoint Change 100)…………. 45
3.35 Input 4 and Input 5 vs Time (Setpoint Change 100)…………. 45

List of Tables
Table 2.1: Interpretation of coupling matrices for different MPM Estimation……… 17

Table 3.1: Constraints on inputs and rate of change of inputs………………………… 23

Table 3.2: Constraints on outputs…………………………………………………… 23

Table 3.3: Step change in the inputs at different instants……………………………. 24

Table 3.4: Transfer function data for step change in Air flow….............................. 29

Table 3.5: Transfer function data for step change in Coal flow …………………… 29

Table 3.6: Transfer function data for step change in Steam flow ………………….. 29

Table 3.7: IAE (Negative Step Disturbance)…………………..…………………… 43

viii
Nomenclature
MPC Model Predictive Control

AMPC Adaptive Model Predictive Control

FSS Feasible System Sets

MPM Model Plant Mismatch

WCHR Char Flow

WAIR Air Flow

WCOL Coal Flow

WSTM Steam Flow

WLS Limestone Flow

CVGAS Fuel Gas Calorific Value

MASS Bed Mass

PGAS Fuel Gas Pressure

TGAS Fuel Gas Temperature

PSINK Sink Pressure

Input 1 Charcoal

Input 2 Air

Input 3 Coal

Input 4 Steam

Input 5 Limestone

Output 1 Calorific Value of the fuel gas

Output 2 Bedmass

Output 3 Pressure of the fuel gas

Output 4 Temperature of the fuel gas

ix
Chapter-1
Introduction
1.1 Background
Model Predictive Control, or MPC, is an advanced method of process control that has been in
use in the process industries such as chemical plants and oil refineries since the 1980s. Most
of the oil companies in the world today use Model Predictive Controllers (Morari and Lee,
1999). Model predictive control refers to a class of algorithms that computes control inputs
by solving a constrained optimization problem over a future horizon (Muske and Rawlings,
1993).

Developed independently by Shell as Dynamic Matrix Control (DMC) in 1979 and Richalet
as Model Algorithmic Control in 1978, these model based control strategies have proved to
be the most widely used in process industries (Muske and Rawlings, 1993, Morari and Lee,
1999, Qin and Badgwell, 2003).

1.2 Motivation for MPC and AMPC

A chemical process is generally characterized by multivariable and coupled interactions. The
system dynamics is highly non-linear and most of the process parameters are constrained
within certain range due to physical limitations of the actuator. Most of the earlier controllers
like PID couldn’t handle multivariable interactions efficiently. This is because a PID
controller typically has one input and one output and thus fails to capture the effects of
multivariable interactions and coupling. Separate Logical based control strategy is required to
co-ordinate among all PID controllers which can lead to chaos and instability. Similarly,
Linear Quadratic Optimal Control gives us optimal results but is unable to handle constraints
on process parameters which we encounter in most of the chemical processes. Hence, an
efficient control algorithm is required which can deal with multivariable interactions, process
non-linearity, constraints on process parameters and at the same time giving optimal results.
The algorithm of Model Predictive Control answers to all these questions. Hence, from
chemical engineering point of view, MPC is very important. The major advantage of MPC
lies in its ability to handle multivariable interactions, process non-linearity, constraints on
process parameters and robustness (Qin and Badgwell, 2003).

When the process is highly non-linear or the system is time varying, then using a linearized
process model with constant parameters to control the plant may not give satisfactory result
over large deviation from operating points. However, it is very difficult to develop a reliable
dynamic non-linear deterministic model for a chemical process. Chemical processes are high
dimensional and distributed parameter systems and uncertainties in the parameters of models
developed from the first principle can be high. Even if such a model developed, using it on-
line may prove difficult due to high computational time (Qin and Badgwell, 2003).

1|Page
Another approach which can deal with the process non-linearity efficiently is that we can
estimate parameters of the process model every time we get a new measurement. The
controller parameters are adjusted in real time based on the input –output data collected. This
method is helpful in controlling a time varying or initially unknown system (Tanaskovic et al,
2013). One of the most efficient control strategies is the Adaptive Model Predictive Control.
The major advantage of Adaptive MPC lies in the idea that model structure still remains the
same and the model parameters adapt with time as the new measurement is available
(Adetola et al, 2009).

1.3 Introduction to MPC and AMPC

Model Predictive Controllers require a good model for future predictions. The entire essence
of MPC lies on how good the model is. This is because MPC relies on forecasting the output
over a moving window and getting the manipulated inputs by solving a constrained non-
linear optimization problem (Morari and Lee, 1999). In most cases, MPC acts as a
supervising controller to co-ordinate among different PID controllers.

The figure 2.1 below shows a pictorial representation of MPC algorithm.

Figure 1.1: Diagrammatic Representation of the MPC Algorithm

2|Page
As we can see in the figure 2.1, the manipulated input is predicted over a moving window.
The future inputs are then used to get the predicted future output trajectory using prediction
equations. The control problem is then solved as a non-linear constrained optimization
problem. The first move of the optimized manipulated inputs is then applied to the plant and
the entire optimization problem is again solved at the next instant (Muske and Rawlings,
1993).

This class of algorithms are also referred to as Receding Horizon control or Moving horizon
control (Muske and Rawlings, 1993).

A typical mathematical formulation of MPC is composed of three different parts. Firstly,

system identification is done where we identify a model which can be used for future
predictions. Secondly, state estimation which includes designing a suitable observer. Thirdly,
an optimization problem needs to be formulated over a moving horizon which is solved using
different non-linear programming techniques to get the manipulated input profile.

Modeling of the system can be done either using first principles, empirically or a combination
of both. Depending on how the modeling is done, we will end up getting a white box model,
grey box model or black box model.

The modeling for the systems in this report is done empirically using Time Series Analysis.
The modeling is described in detail in the subsequent sections

AMPC Algorithm works similar to the MPC algorithm. In AMPC, the parameters of the
models are estimated in the real time. That’s to say, online model identification at each time
instant is performed. Different algorithms have appeared in literatures over the past few years
for parameter adaptation.

One of the most popular algorithms used for online model identification is the concept of the
Recursive Least Square Estimation (Astrom and Wittenmark, 1994). However, constrained
formulation of AMPC is not as straightforward as MPC. In fact, some of the adaptive control
strategies including Recursive Least Square Estimation give poor results when constraints are
involved (Tanaskovic et al, 2013).

1.4 Literatures on MPC and AMPC

Since the use of MPC extended over different fields, different algorithms have appeared in
the literatures for different applications. The earliest works of the MPC were called DMC,
GPC and MAC.

In one of the important works on Model Predictive Control with Linear Models by Muske
and Rawlings (1993), a new model predictive algorithm was developed which predicted the
trajectory over an infinite horizon. This guaranteed nominal stability independent of the
tuning parameters. The work introduced a terminal state constraint and a terminal weighting
which is found by solving a discrete Lyapunov Equation.

3|Page
The field of Adaptive Control has extensively been studied for a long time and there are well
established theories but as far as the constrained adaptive control is concerned, little has been
done in this area. Hence, practical application of adaptive control is limited as most of the
systems are constrained. In one the work on AMPC for constrained linear systems by
Tanaskovic et al (2013), a novel adaptive feedback control technique for uncertain linear
systems is proposed which is able to handle constraints on input and output as well as
measurement noise. The algorithm uses a novel Set Membership identification algorithm to
on-line refine the Feasible System Sets. The computed FSS is then used by the controller for
tracking. The solution of problem requires only Linear and Quadratic programming
techniques. While much have been done in the field of linear AMPC, the non-linear AMPC
has received very little attention. In one of the works for AMPC for constrained non-linear
systems by Adetola et al (2009), the non-linear MPC is combined with adaptation mechanism
for the uncertainty set. The algorithm uses a Min-max and Lipschitz based method approach
to provide robustness for MPC Controllers. In this approach, the parameters are estimated
only when improved estimates of the parameters are available.

1.5 Objectives
In the literature, multiple formulations of MPC and AMPC are available. In particular, there
are several ways of handling the model plant mismatch. Each one of them has their
advantages and disadvantages, and, for practical implementation, a user has to select one that
suits best for the problem under consideration. This part of the work is aimed at development
of a unified MPC formulation, which can be reduced to any of the available formulation by
making appropriate choices of coupling matrices in the generalized MPC formulation. To
achieve efficient on-line implementation, the resulting constrained optimization problem has
been cast as QP problem. It is desired to test the proposed toolbox using the Alstom coal
gasifier as one of the test examples. Development of MPC relevant dynamic model is in
progress.

1.6 Outline of the Report

The rest of the report is organized as follows. Chapter 2 gives an overview of the generalized
MPC formulation. This chapter describes some of the important work that has been done in
the field of MPC and AMPC and a generic algorithm for MPC. The generic algorithm is used
to develop a MPC Toolbox. This chapter also includes generalized QP formulations for MPC
problem over infinite horizon. Chapter 3 contains a case study of the Alstom Benchmark
Challenge II on Gasifier Control. In particular, modeling and simulations aspects of the
system are discussed. Chapter 4 presents conclusion reached so far and scope of the future
work.

4|Page
Chapter-2
MPC Toolbox
This chapter describes the generic algorithm of MPC that can be used to derive all simpler
forms. All the possible cases of MPC formulation has been taken care of. Towards the end of
the chapter, the generic formulation has been mapped to simpler formulations.

2.1 Modeling and Estimation

Let’s consider that our system to be controlled can be represented by the following linear
state space model.

x(k  1)  x(k )  u (k )  w(k )

(2.1)
y (k )  Cx(k )  v(k )

x(k )  R n represents the state variables, u(k )  R m represents the manipulated inputs ,
y(k )  R r represents the measurement variables, w R n represents the state noise and v  R r
represents the measurement noise .

The noises w R n and v  R r are assumed to be zero mean white noises with co-variances
as shown

Q  E[( w( k ) w( k ) T ]
Pwv  E[( w( k )v ( k ) T ] (2.2)
R  E[( v( k )v( k ) T ]

The model described above can be developed either by linearization of a first principle model
or state realization of time series model developed from the input-output data.

The above model can be used to develop a predictor and current state estimator form of state
estimator. These state estimators can either be derived using pole placement design
(Luenberger Observer) or Kalman Filtering approach. In fact, if we develop model from the
data, what we end up getting is the predictor form directly.

Predictor

xˆ (k / k  1)  xˆ (k  1 / k  2)  u (k  1)  Lp e(k  1)
(2.3)
e(k )  y(k )  Cxˆ (k / k  1)

5|Page
Current State estimator

xˆ (k / k  1)  xˆ (k  1 / k  1)  u (k  1)
e(k )  y (k )  Cxˆ (k / k  1) (2.4)
xˆ (k / k )  xˆ (k / k  1)  Lc e(k )

Estimation of plant model Mismatch

In order to account for the model plant mismatch (MPM), the state space model developed
earlier can be augmented with extra artificial states.

State Augmentation

The state space model as shown earlier in Eq. 2.1 is augmented with extra artificial states. It
must be ensured that extra states should not exceed the number of measurements or else the
additional states will not be observable.

The augmented state space model looks like

x(k  1)  x(k )  u u (k )    (k )  w(k )

 (k  1)   (k )  w (k )
(2.5)
 (k  1)   (k )  w (k )
y (k )  Cx(k )  C (k )  v(k )

Here   R s and   R t are the input and output disturbances vectors while vectors w  R s
and w  Rt are the zero mean white noise sequences with covariance Q and Q
respectively.

The artificial states are expected to account for any unmeasured disturbances and the MPM
and hence help in eliminating the offset.

The above state space model can also be written in the augmented form as

xa (k  1)   a xa (k )  uau (k )  wa (k )
(2.6)
ya (k )  Ca xa (k )  v(k )

The individual matrices look like

 x(k )    0 u 


xa (k )    (k )  a  0 I 0 ua  0 (2.7)
 (k )  0 0 I  0


6|Page
Q 0 0  w(k ) 
   
R1a  E[ wa (k ) wa (k ) ]  0
T
Q 0 wa   w (k )
0 0 Q   w (k ) 
  
R2a  E[v(k )v(k )T ]  R

 Pwv 
 E[ wa (k )v(k ) ]   0 
 
T
R12a (2.8)
Ca  C [0] C
 0 

2.2 The Generic MPC formulation

Let us assume that at any instant k , we define our manipulated input for the next p instants
as follows:

U p  {u(k / k ), u(k  1/ k ),.......u(k  p  1/ k )} (2.9)

For carrying out the future predictions, we use the following generic form of predictor

zˆ (k  j  1)  zˆ (k  j )  uu (k  j / k )   ˆ (k  j / k )
ˆ (k  j  1 / k )  ˆ (k  j / k ) (2.10)
ˆ (k  j  1 / k )  ˆ (k  j / k )
yˆ (k  j / k )  Czˆ (k  j )  Cˆ (k  j )

In practice, we restrict the future degree of freedom to q manipulated moves.

U q  {u(k / k ), u(k  1/ k ),.......u(k  q  1/ k )} (2.11)

Using the input blocking concept, we span the q manipulated moves as following

u (k  j )  u (k / k ) : j  m0  1.......m1  1
u (k  j )  u (k  m1 / k ) : j  m1  1.......m2  1
...........................................................
(2.12)
u (k  j )  u (k  mi / k ) : j  mi  1.......mi 1  1
............................................................
u (k  j )  u (k  mq 1 / k ) : j  mq 1  1.......mq  1

where m j ’s are selected such that

m0  0  m1  m2  .....  mq 1  mq  p

The degree of freedom that we have here i.e. q is called the control horizon.

7|Page
Using the model equation for prediction above, the outputs is predicted over the infinite
horizon. In the infinite horizon formulation, the input value is made zero after some time.

Let us assume that the desired future set point trajectory at any instant k is given by

{ yr (k  j / k ) : j  1,2...... p} . (2.13)

Model Predictive Control in Infinite Horizon Formulation

One of the uncertainties that we have while using the finite horizon formulation is that its
stability is not guaranteed. That is to say for some value of tuning parameters, our solutions
may diverge. However, this problem can be overcome by utilizing the infinite horizon
formulation. Gilbert et al have proved the stability of infinite horizon MPC. In practice, we
can’t keep predicting our output over the infinite horizon. Muske and Rawlings (1993) have
developed a novel way of incorporating the infinite horizon formulation over a finite horizon
by using the solution of discrete Lyapunov Equation.

The infinite horizon problem formulation as described in Muske and Rawlings (1993) has
been followed in this report.

Given the set point trajectory, the model predictive control problem at the sampling instant k
is defined as a constrained optimization problem whereby the future manipulated input moves
in the set U q are determined by minimizing the objective function J (k ) defined as follows

[U q ]opt  Min
U q J (k )

J (k )   (k  p / k )T w (k  p / k )

p 1
   (k  j / k )T we (k  j / k )
j 1

q 1
  u (k  m j / k )T wu u (k  m j / k )
j 0

q 1
 [u (k  m j / k )  u s ]T wu [u (k  m j / k )  u s ] (2.14)
j 0

 (k  j / k )  yr (k  j / k )  yˆ (k  j / k ) : j  1,2,....... p (2.15)

represents predicted control error and

u (k  m j / k )  u (k  m j / k )  u (k  m j  1 / k ) : j  1,2,......q  1
(2.16)
u (k  m0 / k )  u (k / k )  u (k  1)

represents changes in the manipulated input moves

8|Page
Here ( we , w ) represents positive definite error weighting matrices, wu represents the
positive semi-definite the input move weighting term and wu represents positive semi-
definite the input weighting matrix. The us (k ) represents a reference manipulated input at the
k ’ th instant. The us (k ) is found by minimizing output tracking error

Min
x s ,u s   ( yt (k )  Cxs (k )  C (k ))T ( yt (k )  Cxs (k )  C (k )) (2.17)

subject to the following constraints

u min  u s  u max
x  (2.18)
I   ...... u  s   0
u s 

Here yt (k ) is the output target and  is a positive definite weighting matrix for output
tracking error.

The above minimization problem can be formulated as QP where

H  uT I   C T CI  u

T

  
(2.19)
F (k )  2 uT I   C T  yt  CI      C
T 1

subject to: As u s  Bs (2.20)

 I  u 
As    ; Bs   max 
 I   umin 

When the poles of  are inside unit circle, the terminal state weighting matrix can be found
by solving the discrete Lyapunov Equation as

w  C T weC  T w (2.21)

The MPC minimization problem is subject to the following constraints

Model Prediction Equations

zˆ(k  1)  zˆ(k )  u u (k )   ˆ (k )
(2.22)
yˆ (k )  Czˆ(k )  Cˆ (k )

9|Page
Manipulated input constraints

u L  u (k  m j / k )  u H : j  0,1,2,.........q  1
(2.23)
u L  u (k  m j / k )  u H : j  0,1,2,.........q  1

Output Constraints

yL (k  j )  yˆ (k  j / k )  yH (k  j ) (2.24)

The resulting optimization problem can be solved by using any nonlinear programming
method. The controller is implemented in the moving horizon framework. Thus, after solving
the optimization problem, only the first move uopt (k / k ) is implemented on the plant, i.e.

u(k )  uopt (k / k ) (2.25)

2.3 QP Formulation of MPC

Initially, let’s assume that the control horizon (q) and the prediction horizon ( p) are equal.
Defining the future input vector U p (k ) , the target steady state input vector U s (k ) and
predicted output trajectory Yˆ (k ) as follows


U p (k )  u (k / k )T u (k  1 / k )T ........ u (k  p  1 / k ) 
T T


U s (k )  u s (k )T u s (k )T ........ u s (k ) 
T T
(2.26)


Yˆ (k )  yˆ (k  1 / k )T yˆ (k  2 / k )T ........
T
yˆ (k  p  1 / k ) zˆ (k  p)T 
T

The prediction model can now be expressed as a single vector equation

Yˆ (k )  S x zˆ(k )  SuU p (k )  S  ˆ (k )  Sˆ (k ) (2.27)

Here the block matrices are as shown below

 C   Cu [0] .... [0] [0]

 C 2   C  Cu .... 
   u 
S x   ....  Su   ....... ......  (2.28)
 p 1   p 2 
C  C u C u 
  p    p 1u  p  2 u u 

C   C 
C   C (  I n ) 
   
S    ....  S   ....  (2.29)
   p2 p 3 
C  C (    .......  I n ) 
 0  ( p 1   p  2  .......  I n ) 
 

10 | P a g e
The matrix S u is often referred to as Dynamic Matrix of the system.

In order to incorporate the input blocking constraints, the matrix S u needs to be modified.
Consider matrices {I mi : i  1,2,......, q} each of the dimensions (mi  mi 1  1)m  m and defined
as follows

I m 
I 
I mi   m (2.30)
 ... 
 
 I m  ( mi  mi1 1) mm

Let’s also define the transformation matrix pq as

pq  blockdiag I m1  I m1 .... I m1 

( pm pm)
(2.31)

Here I m represents m m identity matrix. Now, let us define the future input vector over the
control horizon (q) .


U q (k )  u(k / k )T u ( k  1 / k )T ........ u (k  q  1 / k ) 
T T
(2.32)

Using the input blocking constraints, the future input moves over the control horizon can be
defined over the prediction horizon as shown

U p (k )  pqU q (k ) (2.33)


U p (k )  u(k / k )T ...... u(k / k )T .......... u(k  mq 1 / k )........ u(k  mq 1 / k ) 
T T
(2.34)

The vector form of the output trajectory is now modified to be consistent with the input
moves over the control horizon. The modified equation looks like

Yˆ (k )  S x zˆ(k )  SuU q (k )  S  ˆ (k )  Sˆ (k ) (2.35)

The above equation uses the modified form of the S u matrix defines as follows

 Cu [0] .... [0] [0]

 C  Cu .... 
 u 
S u   ....... ......   pq (2.36)
 p2 
C u C u 
  p 1u  p  2 u u 

11 | P a g e
Defining the future setpoint trajectory vector R(k )


R ( k )  y r ( k  1 / k )T y r ( k  2 / k )T ...... yr (k  p  1 / k )T xs (k ) 
T T
(2.37)

Here xs (k ) is the target steady state state vector

Defining the predicted error vector E (k ) at any instant k as follows

E (k )  R(k )  Yˆ (k ) (2.38)

The MPC objective function J (k ) can now be re-formulated as follows


 E (k ) WE E (k )  U q (k ) WU U q (k ) 
T T

J (k )    (2.39)
 [U q (k )  U s (k )] WU [U q (k )  U s (k )]

T


Here the weighting matrix are described as follows

WE  blockdiag we we ........... w ( pn pn)

WU  blockdiag wU  wU ........... wU 
( qmqm)
(2.40)
WU  blockdiag wU  wU ........... wU 
( qmqm)

Here U q (k ) is defined as follows

U q (k )  U q (k )  0u(k  1) (2.41)

 Im [0] [0] [0]  Im 

 I Im [0] [0] [0]
 m 0    (2.42)
   ... 
   
 [0]  Im I m  ( qmqm) [0] ( qmqm)

The conventional MPC problem can now be represented as following


 E (k ) WE E (k )  U q (k ) WU U q (k ) 
T T

J (k )    (2.42)
 [U q (k )  U s (k )] WU [U q (k )  U s (k )]

T


subject to the following constraints

Yˆ (k )  S x zˆ(k )  SuU q (k )  S  ˆ (k )  Sˆ (k ) (2.43)

YL  Yˆ (k )  YH
u u L  U q (k )  u u H (2.44)
u u L  U q (k )  u u H

12 | P a g e
Here the matrix u is defined as

I m 
I 
u   m  (2.45)
....
 
 I m  ( qmm )

which contains I m stacked times and


YL  yL (k  1)T yL (k  2)T ........ yL (k  p  1)T zˆL (k  p)T 
T

(2.46)
 y 
T
H ( k  1) yH (k  2)T yH (k  p  1)T zˆH (k  p)T
T
YH ........

The above MPC objective function can be reduced to the conventional Quadratic Objective
function as follows:

1
J (k )  U q (k )T HU q (k )  U q (k )T F (k ) (2.47)
2

Here the matrices

H  2( S uT WE S u   T WU   WU )

F (k )  2 ( S uT WE ) (k )  ( T WU 0 )u (k  1)  WU u u s (k )  (2.48)

 (k )  R(k )  S zˆ(k )  S ˆ (k )  S ˆ (k )
x   

Now, the Quadratic Programming formulation of MPC can be defined as

1
min U q (k )T HU q (k )  U q (k )T F (k ) (2.49)
2
U q (k )

subject to: AU (k )  B

Here

 I qm   u u H 
 I    u u L 
 qm   
    u u H  0 u (k  1) 
A  B  (2.50)
     u u L  0 u (k  1) 
 Su   YH  S x zˆ (k )  S  ˆ (k )  Sˆ (k ) 
   ˆ 
  S u   YL  S x zˆ (k )  S   (k )  Sˆ (k )

13 | P a g e
Future Setpoint Trajectory
The set point trajectory can be generated using two approaches.

Using a reference system of the form

xr (k  j  1/ k )   r xr (k  j / k )  I   r r (k )  y(k )
yr (k  j  1 / k )  y(k )  Crf xr (k  j  1 / k ) (2.51)

for: j  0,1,...... p  1
 r  diag 1  2 ......  r  (2.52)
0  1  1 for i  0,1,......r are tuning parameters
Here
Initial condition: xr (k / k )  0
r (k )  R r represents the setpoint vector
Typically : 0.8  1  0.99 (2.53)

The coefficient matrices in the setpoint trajectory equation should be selected such that the
steady state gain of the reference system should be equal to Identity matrix i.e.

Crf ( I   r ) 1 r  I (2.54)

Typically, the transfer function of matrix of the reference system is diagonal with unity gain
first order low pass filter.

The reference system trajectory can be selected alternatively as follows

Let’s define the signal rf (k ) as

rf (k )   r rf (k  1)  I   r r (k )
r f (0)  0
yr (k  j / k )  rf (k ) for j  1,2,..... p (2.55)
Typically: 0.8   r  0.99

Matlab implementation of the above formulation is included in the Appendix A

In the first method, the filtering of the trajectory is done inside the controller and the
controller is aware of the global setpoint while in the second case, a filtered trajectory is
already fed to the controller and the controller maintains the setpoint fed to it a particular
instant, i.e. controller is unaware of the global setpoint.

14 | P a g e
2.4 Approaches to MPM Estimation
Here, in this section, it is shown that formulation discussed above can describe all the simpler
formulations of MPC available in the literature.

The choice of artificial states and the coupling matrices can be done in different ways.

State Augmentation Approaches:

 Output bias formulation (or Dynamic Matrix Control (DMC)):

In this approach, the elements of vector  can be viewed as a bias in r measured outputs.
The drifting disturbances are assumed to be affecting only the measured outputs as additive
biases. This can be achieved by specifying the coupling matrices as follows

  0 ; C  I (2.56)

 Input bias formulation:

In this approach, the elements of the vector  can be viewed as bias in r manipulated inputs.
When the number of manipulated inputs equal the number of measured outputs ( r  m ), then
we can choose

  u ; C  0 (2.57)

 Disturbance bias formulation:

When the state space model is derived from the first principle, it is possible to choose

  d (2.58)

provided number of disturbance variable (d  r ) . If d  r , then r linearly independent

column of d can be chosen as  .

The augmented state space model can then be used to develop the observers.

Prediction Estimator

xˆ a (k / k  1)   a xˆ a (k  1 / k  2)  uau (k  1)  L paea (k  1)
(2.59)
ea (k )  y(k )  Ca xˆ a (k / k  1)

Current State estimator

xˆ a (k / k  1)   a xˆ a (k  1 / k  2)  uau (k  1)
ea (k )  y (k )  Ca xˆ a (k / k  1) (2.60)
xˆ a (k / k )  xˆ a (k / k  1)  Lca ea (k )

15 | P a g e
In all the three cases described above, the artificial states  (k ) and  (k ) are estimated along
with the original state x(k ) using the observer.

In another approach, the artificial states  (k ) and  (k ) are estimated without designing an
augmented state estimator.

 Innovation Bias Approaches

In these formulations, the filtered innovation and model residuals are used for estimation of
artificial states  (k ) and  (k ) depending upon whether we are using a predictor based
approach or current state estimator based approach.

The vectors  (k ) and  (k ) contain information about the plant model mismatch as well as
the measurement noise. In order to eliminate the high frequency measurement noise and to
shape the controller’s response to the model plant mismatch, the innovation and residuals are
filtered to estimate the  (k ) and  (k ) .

Prediction Estimator based MPM Estimation

  L p ; C  I (2.61)

The artificial states  (k ) and  (k ) are estimated by the following low-pass filter equation

ˆ (k )    ˆ (k  1)  [ I    ]e(k )
ˆ (k )  ˆ (k ) (2.62)
ˆ (0)  ˆ (0)  0

   diag[1 2 ... r ]

Here 0   i  1 and i  1,2,....., r are the tuning parameters.

Current State Estimator based MPM estimation

  Lc ; C  I (2.63)

The artificial states  (k ) and  (k ) are estimated by the following low-pass filter equation

ˆ (k )    ˆ (k  1)  [ I    ]e(k )
ˆ (k )  ˆ (k  1)  [ I   ][ y (k )  Cxˆ (k / k )] (2.64)
ˆ (0)  ˆ (0)  0

   diag[ 1 2 ... r] (2.65)

Here 0   i  1 and i  1,2,....., r are the tuning parameters.

16 | P a g e
Another way of innovation bias is to use the augmentation and estimate  (k ) by designing
an augmented state observer and then add it as a bias in innovation.

 Innovation Bias with Augmentation

The following formulation explains it clearly

z (k  1)  z (k )  u u (k )  L[  (k )  e(k )]
 (k  1)   (k )  w (k ) (2.66)
y (k )  Cz (k )   (k )  e(k )

From above it is clear that in order to filter  (k ) , we are using a state estimator and  (k ) is
then added as a bias to the innovation e(k ) .

 DMC Approach

This is perhaps one of the oldest approaches to estimate plant model mismatch. In this
approach, an open loop observer is used and the entire plant model mismatch is added to the
output as follows

xˆ (k  1)  xˆ (k )  u u (k )
(2.67)
yˆ (k )  Cxˆ (k )  ˆ (k )

The proposed formulation reduces to the DMC formulation if we let

C  I ;   0 (2.68)

Thus, we see above that there are around 11 possible cases for estimation of plant model
mismatch.

Following Table 2.1 lists all the possibilities discussed above

Table 2.1 Interpretation of coupling matrices for different approach of MPM Estimation

Estimation of MPM  C
Input Bias u 0
State Augmentation Output Bias 0 I 
Disturbance Bias d 0
Predictor Form Lp I 
Innovation Bias
Filter Form Lc I 
DMC ( L =0) 0 I 

17 | P a g e
In order to use the state augmentation for estimation of the plant model mismatch, it is
important to consider the cross-covariance between the state and measurement noise. In order
to develop a state estimator for such a system, a generalized formulation of Kalman filter
which takes into account the cross-covariance is needed.

Using this generalized formulation of Kalman Filter (Astrom and Wittenmark, 1994), the
state estimate has been done for state augmentation approaches.

2.5 Adaptive Control

Adaptive Control is a very important control algorithm for dynamical systems where the
system parameters are function of time. The algorithm estimates the model parameters at
different time instants so that model prediction is still good enough. There are different
algorithms for Adaptive Control. One of the very popular algorithms that can be applied on
Linear Time Series models is the Recursive Least Square estimation.

Recursive Least Square Estimation

Recursive computational methods are generally used when the measurements are obtained
sequentially. Rather than solving the conventional least square problem again for all the
measurements, this algorithm helps us to solve the problem just by including the new
measurement and using the past results. As described in Astrom and Wittenmark (1994), the
algorithm works in the following way.

Let  be the true value of the model parameters.

 = a1  cn 
T
an b11 b21  b1n b2n c1 

(2.69)

Let  (N ) denote the least square measurement based on N measurements. Then the
algorithm for the recursive estimation using least square is given by
  
 ( N  1)   ( N )  K ( N )[ yN 1   T ( N  1)( ( N )] (2.70)

K ( N )  P( N ) ( N  1)[1   T ( N  1) P( N ) ( N  1)]1 (2.71)

P( N  1)  [ I  K ( N ) T ( N  1)]P( N ) (2.72)

Here, K (N ) is termed as weighting factor and P is proportional to the variance of the

estimates.

y N 1 is the ( N  1)th measurement

The least square estimate can be interpreted as Kalman Filter for the process

18 | P a g e
 (k  1)   (k ) (2.73)

y(k )   T (k ) (k )  e(k ) (2.74)

 T (k )   y(k  1)   y (k  n) u(k  1)  u ( k  n ) (2.75)

The concept of Recursive Least Square Estimation is widely used in the field of Adaptive
Model Predictive Control to estimate the parameters of the Time Series Models.

2.5.1 Results and Discussions

The Quaruple tank system is first modelled by second order ARMAX model. This initial
model is then used for Adaptive Control. The following plots show the variations of outputs
and inputs with time.

Figure 2.1: Height in Tank (cm)

19 | P a g e
Figure 2.2: Input (V)

2.6 Summary
In this chapter, a generic formulation for MPC with MPM compensation has been developed.
With a suitable choice of model parameters, this formulation can be reduced to any of the
available forms of MPM compensation discussed in the MPC literature. The chapter also
shows the results of Adaptive Control using the recursive least square algorithm.

20 | P a g e
Chapter-3
Alstom Benchmark Challenge
This chapter describes the case study of Alstom Benchmark Challenge II on gasifier control
as mentioned in the Dixon and Pike (2006). The chapter briefly describes the process of
gasification, control system specification and performance tests that need to be carried out.
Towards the end, the results for the modeling of the gasifier are shown.

3.1 Introduction
The origin of the Alstom Gasifier control problem dates back to the late 90’s when efforts
were made in UK to analyse the feasibility of a small scale power generation plant which is
based on Air Blown Gasification Cycle (Dixon and Pike, 2006). The analysis also included
the modeling and the control aspect of the plant. The coal gasifier was one of the major
components for the power plant. The gasifier was found to be a multivariable system with
four outputs and five inputs and exhibiting highly coupled and stiff dynamics.

Gasification Process

The gasifier can be considered as a reactor where coal is gasified with steam. Coal is mixed
with limestone which absorbs sulphur from the coal and then injected into the gasifier at high
pressure. The air and injected steam fluidize the solids in the gasifier and reacts with the
carbon and volatiles from the coal, thus producing a low calorific value gas. The remaining
char is then removed (Dixon et al, 2000). Gasifier offers a clean and efficient way of
producing energy from high sulphur oil, refinery tars, petroleum coke etc (Dixon and Pike,
2006). The gas so obtained as a result of gasification is then used to drive turbines to produce
energy.

Figure 3.1 shows a schematic of a typical gasifier as mentioned in Dixon et al, 2000.

Alstom Benchmark Challenge I

The first problem on the gasifier control was based on the three linear first principal models
that were developed at different load points (Dixon et al, 2000). The linear models were
derived by linearizing the non-linear gasifier model. The challenge included designing
control algorithms that can be used for regulatory problems complying with the inputs and
outputs constraint as well as specifications for the disturbance tests.

Alstom Benchmark Challenge II

The second gasifier control problem was introduced with the first principal non-linear model
of the gasifier in Matlab/Simulink. The non-linear model takes into account all the physical
process that takes place during gasification, for eg, drying process, desulphurization process,
pyrolysis process, gasification process, mass and energy balance. The non-linear model has
been validated against a wide range of industrial data (Dixon and Pike, 2006).

21 | P a g e
Figure 3.1: Schematic model of gasifier

The five controllable inputs of the plant are coal flow, charcoal flow, air flow, steam flow and
limestone flow. In closed loop, the limestone flow will be 10 percent of the coal flow as it is
used to absorb the sulphur in the coal. The four controlled outputs are the calorific value of
the fuel gas, fuel gas temperature, fuel gas pressure and bed mass. In addition, there are two
disturbances: PSINK which represents pressure changes induced as position of the gas
turbine inlet of the fuel gas changes and the change in the coal quality.

3.2 Control System Specification

The challenge requires that the given non-linear model of the gasifier should be treated as the
plant (Dixon and Pike, 2006). The problem needs to be solved by first modeling the system
from scratch and them implementing the control algorithm. However, while implementing
the controller, certain constraints on the inputs and outputs has to be satisfied.

Table 3.1 and Table 3.2 lists the constraints on the inputs and outputs respectively

22 | P a g e
Table 3.1: Constraints on inputs and rate of change of inputs

Inputs (kg/s) Max (kg/s) Min (kg/s) Rate, kg/s

WCHR 10 0 0.2
WAIR 20 0 1.0
WCOL 6.0 0 1.0
WSTM 3.5 0 0.2

Table 3.2: Constraints on outputs

Outputs Constraints (in absolute value)

CVGAS (KJ/kg) 10 (KJ/kg)
MASS (bar) 0.1 (bar)
PGAS (kg) 500 (kg)
TGAS (oC) 1 (oC)

Controller Performance:

The controller has to be subjected to different tests: The tests included the pressure
disturbance tests, load change tests and change in coal quality tests. This is to ensure that
controller is robust to these changes.

Pressure Disturbance Tests:

The pressure disturbance tests were the following (Dixon and Pike, 2006).

 Firstly, giving a step disturbance of -0.2 bar to the system at t=30s and running the
simulation for 5 minutes.
 Secondly a sine wave disturbance of Amplitude 0.2 bar and frequency of 0.04 Hz.
The simulation needs to be run for 5 minutes. The frequency of 0.04 Hz corresponds
to low frequency movements of the inlet valve in reaction to changes in grid
frequency.

The pressure test has to be performed for the 100%, 50% and 0% load operating point.

Load Change Tests (Dixon and Pike, 2006)

The load change test is important to ensure that controller works well for different operating
points of the plant. It may be possible that due to change in demand, we may not be operating
our plant at the same operating point which is used to design the controller. The test requires
that change in load must be restricted to 5% per minute.

Coal Quality Tests (Dixon and Pike, 2006)

This test is required to make sure that controller is robust to changes in coal quality. Hence,
the control algorithm should also take care of the effect of an incremental change in the coal
quality upto 18%.

23 | P a g e
Benchmark PI Controller (Dixon and Pike, 2006).

The paper also presents the results obtained by controlling the gasifier using a set of PI
controllers. These results have to be used as benchmark while developing the candidate
controller.

3.3 Results and Discussions

Modelling of the Gasifier

The gasifier has an inherent PI control to control Bedmass by manipulating Charcoal flow.
This is because Bedmass is marginally stable. In the closed loop, the limestone flow is a fixed
ratio of coal flow which is 0.1 as it is used to absorb sulphur. Hence, out of 5 manipulating
inputs, only 3 degrees of freedom are left. The model identification is thus reduced to three
outputs (CVGAS, PGAS and TGAS) and three inputs (WAIR, WCOAL and WSTM). The
modeling of the gasifier is first done using step response model. This is because identification
of ARX or ARMAX, either the model order has to be increased or Time delay information is
required. The non-linear plant is run in the open loop by supplying input in the form of
rectangular pulse. The band of the pulse input is in the order of 10 percent of absolute values
of inputs at the operating points. This is because plant is highly non-linear and linear
approximation will be valid only near the operating points. A large input band may cause the
plant to drift away from the operating region and hence the linear model so developed may
lead to large model-plant mismatch and hence may not give satisfactory result while
implementing MPC. The data thus obtained is then fitted to a suitable SISO transfer function
model. The 9 FOPTD or SOPTD models obtained is arranged to form a transfer function
matrix. The transfer function matrix so obtained is then converted to a state space model.

The open loop input rectangular pulse as applied is shown in Table 3.1. The deviational form
of inputs and outputs are used for the system identification. The sampling time is 30 sec.

Table 3.3 Step change in the inputs at different instants

Time (sec) Air (kg/s) Coal(kg/s) Steam(kg/s)

0 0 0 0
20000 0.87 0.42 0.135
50000 -1.74 -0.84 -0.27
80000 0.87 0.42 0.135

Using step response data, 9 SISO transfer functions are identified. A transfer function matrix
is formulated using the SISO models. Transfer function matrix is then converted to Minimal
State Space realization. The state space model so obtained is of order 10.

The following figure shows the model fit for different SISO models.

24 | P a g e
Figure 3.2: Calorific value of the fuel gas corresponding to rectangular pulse in Air flow

Figure 3.3: Pressure of fuel gas corresponding to the rectangular pulse in Air flow

25 | P a g e
Figure 3.4: Temperature of the fuel gas corresponding to rectangular pulse in Air flow

Figure 3.5: Calorific Value of the fuel gas corresponding to rectangular pulse in Coal flow

26 | P a g e
Figure 3.6: Pressure of the fuel gas corresponding to rectangular pulse in Coal flow

Figure 3.7: Temperature of the fuel gas corresponding to rectangular pulse in coal flow

27 | P a g e
Figure 3.8: Calorific Value of the fuel gas corresponding to rectangular pulse in steam flow

Figure 3.9: Pressure of the fuel gas corresponding to rectangular pulse in steam flow

28 | P a g e
Figure 3.10: Temperature of the fuel gas corresponding to rectangular pulse in steam flow

The following table lists the gains, zeroes, poles and time delay of the transfer functions

Table 3.4 Transfer function data (Rectangular Pulse in Air)

Outputs Kp Tz T p1 Tp 2 Td
CVGAS -44121 2965.8 1004.3 0 8.5683
PGAS 13309 877.82 1105.2 0 28.173
TGAS 20.75 131.68 928.15 0 10.251

Table 3.5 Transfer function data (Rectangular Pulse in Coal)

Outputs Kp Tz T p1 Tp 2 Td
CVGAS 59834 1450.7 727.97 0 0.10
PGAS 1042.3 1977.7 709.98 0 0
TGAS -17.157 219.92 720.89 0 8.82

Table 3.6 Transfer function data (Rectangular Pulse in Steam)

Outputs Kp Tz T p1 Tp 2 Td
CVGAS -47497 -3897.2 936.67 24.83 0
PGAS 7649.7 1780.6 854.25 0 0
TGAS -46.667 -11.44 934.3 0 76.77

29 | P a g e
The transfer function model so identified has the following structure.

K p (1  Tz s )
y ( s) 
(1  Tp1s)(1  Tp 2 s )

The continuous state space model parameters obtained after the transfer function is converted
to the State Realization are shown below. This state space model is discretized using a
sampling interval of 30 sec.

dx/dt  Ax  Bu
y  Cx  Du

  0.0009957 0 0 0 0 0 0 0 0 0 
 0  0.0009048 0 0 0 0 0 0 0 0 
 
 0 0  0.001077 0 0 0 0 0 0 0 
 
 0 0 0  0.001374 0 0 0 0 0 0 
 0 0 0 0  0.001408 0 0 0 0 0 
A 
 0 0 0 0 0  0.001387 0 0 0 0 
 0 0 0 0 0 0  0.04134  0.005503 0 0 
 
 0 0 0 0 0 0 0.007813 0 0 0 
 0 0 0 0 0 0 0 0  0.001171 0 
 
 0 0 0 0 0 0 0 0 0  0.00107

 8 0 0 
 2 0 0 
 
0.125 0 0 
 
 0 8 0 
 0 2 0   1.303  105 1.192  105 0 
B   
 0 0.125 0  D   1.057  105 2.903 2173
 0 0 128   2.944  5.234 0.21 
  
 0 0 0 
 0 0 2 
 
 0 0 0.25

10.7 0 0  10.2 0 0 62.17  2.042 0 0 

C   0 1.239 0 0  1.311 0 0 0 3.206 0 
 
 0 0 0.1535 0 0  0.1323 0 0 0  0.2007

30 | P a g e
Control of the Gasifier
The Alstom Challenge on gasifier control included linear models developed by linearizing the
first principle non-linear model at 0, 50 and 100 percent operating points. Firstly, the gasifier
is controlled using the linear model and linear plant. Later, linear model is used to control the
non-linear plant. It is not necessary to include the PI controller in the linear plant. Hence, the
linear plant and linear model control doesn’t include any PI control. The PI controller is also
removed for the non-linear plant linear model simulation.

Linear Model and Linear Plant MPC (100% Operating Point)

The following figures shows the results obtained for Linear Plant and Linear Model MPC
(100 % operating point). It can be seen that although the negative step disturbance is rejected,
there is violation of constraints in Pressure of the fuel gas. The Calorific value and
temperature of the fuel shows oscillatory behavior. From figures, we can conclude that in
absence of PI controller, the Bedmass drifts away from the operating point. Moreover, the
settling time of the outputs differ drastically. While the pressure of fuel gas settles in less than
15 minutes, the Temperature of the fuel gas takes approximately 2.5 hours to settle. In case of
inputs, most of them stay within the constraint, but the coal mass flow hits the upper
constraint.

Figure 3.11: Output1 and Output 2 vs Time

31 | P a g e
Figure 3.12: Output 3 and Output 4 vs Time

Figure 3.13: Input 1 vs Time (Closed loop)

32 | P a g e
Figure 3.14: Input 2 and Input 3 vs Time (Closed loop)

Figure 3.15: Input 4 and Input 5 vs Time (Closed loop)

33 | P a g e
Linear Model and Linear Plant MPC (50%)

From the following figures, it can be inferred that in case of 50% operating point, output
constraints are violated in a number of cases and hence it is necessary to solve a constrained
optimization problem specifying the output constraints as well.

Moreover, the settling time of Bedmass and Temperature is very large. However, tuning
parameters may be adjusted to give a small settling time. The response of Bedmass and
temperature is oscillatory as well.

Figure 3.16: Output1 and Output 2 vs Time

34 | P a g e
Figure 3.17: Output 3 and Output 4 vs Time

Figure 3.18: Input 1 vs Time (Closed loop)

35 | P a g e
Figure 3.19: Input 2 and Input 3 vs Time (Closed loop)

Figure 3.20: Input 4 and Input 5 vs Time (Closed loop)

36 | P a g e
Linear Model and Non-Linear Plant MPC (50%)

The MPC algorithm is applied for rejecting the negative step disturbance for the case of Non-
Linear plant and linear model as well. The following figures shows the results obtained at the
50% operating point. It can be inferred from the figure that when the negative step
disturbance is applied, Calorific value, Pressure and Temperature of the fuel gas violates the
constraints. However, Bedmass although stays within constraints but takes a very long time
(approx.) to settle. The manipulated variables stay within constraints.

Figure 3.21: Output 1and Output 2 vs Time

37 | P a g e
Figure 3.22: Output 3 and Output 4 vs Time

Figure 3.23: Input 1 vs Time (Closed loop)

38 | P a g e
Figure 3.24: Input 2 and Input 3 vs Time (Closed loop)

Figure 3.25: Input 4 and Input 5 vs Time (closed loop)

39 | P a g e
Non-linear plant and Linear Model MPC (100 %)

The following figures show the results for MPC implementation at 100 % operating point.
Compared to 50 % operating point, in this case the constraint violation occurs only in the
Pressure of the fuel gas. However, the constraint is violated only in the immediate
neighborhood of the time when negative change in pressure disturbance is applied. The
bedmass in the open loop is controlled by a PI controller and hence bedmass oscillates around
the setpoint before settling. Most of the manipulated variables stay within constraint.
However, coal and limestone flow hits the uppers constraints.

Figure 3.26: Output 1 and Output 2 vs Time

40 | P a g e
Figure 3.27: Output 3 and Output 4 vs Time

Figure 3.28: Input 1 vs Time (Closed loop)

41 | P a g e
Figure 3.29: Input 2 and Input 3 vs Time (Closed loop)

Figure 3.30: Input 4 and Input 5 vs Time (Closed loop)

42 | P a g e
Table 3.7 shows the values of IAE for the three cases considered. It can be seen that IAE for
the 50 % operating point is considerably high which can also be inferred from the plots. IAE
for non-linear plant at 100% is more than linear plant at 100 % which is expected due to
model plant mismatch.

Table 3.7: Integral Absolute Error (Negative step disturbance)

Output Linear (100%) Linear (50%) Non-linear (100%) Non-linear (50%)

CVGAS(J/Kg) 2.9*10^8 1.7102*10^9 5.52*10^8 1.13*10^12
Bedmass (kg) 1.66*10^4 8.08*10^7 7.48*10^4 3.10*10^7
PGAS(N/m^2) 4.25*10^8 1.2674*10^9 3.95*10^9 2.08*10^10
TGAS(K) 10.44 240.7604 59.02 4.08*10^4

Step Change in Setpoint (Linear Model and Non-Linear Plant MPC (100%))

Apart from the control of the gasifier in case of disturbance, an attempt is also made to
control the gasifier by providing a setpoint change in few outputs. The setpoint change is
given only for CVGAS, PGAS and TGAS. The corresponding magnitude of change is 0.2
KJ/kg, 0.02 bars and 0.2 K respectively. Bedmass value is kept constant. The CVGAS shows
inverse response behavior. The output initially shows deviation in opposite direction before
finally settling at the setpoint. Similar inverse response behavior is there in case of TGAS

Figure 3.31: Output 1 and Output 2 vs Time

43 | P a g e
Figure 3.32: Output 3 and Output 4 vs Time

Figure 3.33: Input 1 vs Time (Closed Loop)

44 | P a g e
Figure 3.34: Input 2 and Input 3 vs Time (Closed Loop)

Figure 3.35: Input 4 and Input 5 vs Time (Closed Loop)

45 | P a g e
3.4 Summary
In this chapter, a brief summary of the Alstom Coal gasifier control problem has been
presented. A MIMO output error model is developed by analyzing the positive and negative
step response data. This preliminary model has given an idea about approximate process
gains and time constants. This information will later be used to design PRBS inputs to excite
the plant and develop MPC relevant model that captures the effect of manipulated inputs as
well as unmeasured disturbances. The chapter also covers partial solution to Alstom gasifier
challenge problem. The candidate MPC algorithm using linear model shows satisfactory
result for some operating points. The negative step disturbance is regulated by the MPC.
This result will later be extended by using the identified model for Adaptive Control.

46 | P a g e
Chapter-4
Conclusion
4.1 Summary
In the literature, based on the approach for handling MPM, multiple formulations of MPC are
available. This part of the work was aimed at the development of a unified MPC formulation,
which can be reduced to any of the available formulation by making appropriate choices for
the coupling matrices in the generalized MPC formulation. To achieve efficient on-line
implementation, the resulting constrained optimization has been cast as a QP problem.
Development of a Matlab Toolbox based on the unified formulation is still posing some
difficulties in the development.

It is desired to test the proposed Toolbox using Alstom coal gasifier as one of the test
examples. A preliminary MIMO model with OE structure has been identified by carrying out
open loop simulation studies on the Alstom gasifier simulator. The model identification
exercise has revealed important features of the dynamics such as presence of Right and Left
Half Plane zeroes and non-linearity in the dynamics.

These simulation studies on Alstom Gasifier highlights some of the important things to be
kept in mind while designing input excitation for the system identification and hence the
controller for the system.

From the modeling exercise, it is inferred that modeling is one of the most important step in
designing the controller. Different modeling technique can be used depending upon the
information. While modeling a system from scratch, OE+ARMA model can be a good
option. This model can form a basis for developing more complex model such as ARMAX. It
can also be seen that while developing OE model, the step changes must not be very large or
else model may give poor result away from the operating point.

The control results shows satisfactory trend for some operating points but needs to be
improved to take into account all the operating points and rejects all kinds of measured or
unmeasured disturbances.

4.2 Future Work

The main aim of this work is to develop a Tool-Box for Adaptive MPC. Later, MPC will be
combined with adaptive algorithm.

Once a suitable linear time series model is identified, Adaptive Control will be used for
disturbance rejection. For modeling the Alstom Problem, some non-linear system
identification can be used which takes into account the multivariable interactions. Later, the
other tests mentioned in the gasifier challenge will be completed.

47 | P a g e
Appendix A

Generic Codes for MPC

1. mpc_structure- This function creates a structure named MPC

% This function returns a mpc structure which is the fed to mpc_controller
% function for getting the control input

function
mpc=mpc_structure(mpc,n_op,n_ip,setpoint_filter,rk,rk_f,phy_r,phy,gamma_u,L
,C_mat,yk,uk,xkhat_p,xkhat_c,Pkhat_aug,betakhat,etakhat,phy_beta,phy_eta,Q_
beta,Q_eta,p,q,MPM_type,state_aug_type,Estimator_type,R_mat,we,wu,wdelu,uk_
l,uk_h,yk_l,yk_h,xk_l,xk_h,deluk_l,deluk_h,uk_guess)
% function
mpc=mpc_structure(mpc,rk_f,phy,gamma_u,gamma_beta,C_mat,C_eta,yk,xkhat,we,w
u,wdelu,uk_l,uk_h,yk_l,yk_h,delu_l,delu_h)

mpc=struct('n_op',n_op,'n_ip',n_ip,'setpoint_filter',setpoint_filter,'rk',r
k,'rk_f',rk_f,'phy_r',phy_r,'phy',phy,'gamma_u',gamma_u,'L',L,'C_mat',C_mat
,'yk',yk,'uk',uk,'xkhat_p',xkhat_p,'xkhat_c',xkhat_c,'Pkhat_aug',Pkhat_aug,
'betakhat',betakhat,'etakhat',etakhat,'phy_beta',phy_beta,'phy_eta',phy_eta
,'Q_beta',Q_beta,'Q_eta',Q_eta,'p',p,'q',q,'MPM_type',MPM_type,'state_aug_t
ype',state_aug_type,'Estimator_type',Estimator_type,'R_mat',R_mat,'we',we,'
wu',wu,'wdelu',wdelu,'uk_l',uk_l,'uk_h',uk_h,'yk_l',yk_l,'yk_h',yk_h,'xk_l'
,xk_l,'xk_h',xk_h,'deluk_l',deluk_l,'deluk_h',deluk_h,'uk_guess',uk_guess);

2. mpc_controller- This function calls the individual subroutines.

% This function returns the value of contraol input uk which is obtained by
% solving Linear MPC problem using Quadratic programming. The function
% contains several nested functions which do the individual work as
% explained before them and inside
function mpc=mpc_controller(mpc)

% This function returns the estimates of model-plant mismatch variables and

% constants like betakhat, etakhat, gamma_beta and C_eta;
mpc=mpc_mod_plant_mismatch(mpc);

% This function returns the setpoint trajectory. Note that the trajectory
% can be filtered inside or outside. For filtering inside case, it is
% filtered in the function
mpc=mpc_setpoint_trajectory(mpc);

% This function calculates the target state which is required for solving
% the infinite horizon MPC problem.
mpc=mpc_steady_state_target(mpc);

% This function facilitates in building up the block matrices which are

% required in the formulation of quadratic programming problem.
mpc=mpc_block_mat_quadprog(mpc);

% This function returns the handles of quadprog like Hessian, F vector,

% A_mat, b_vec etc in terms of the block matrices defined in the previous

48 | P a g e
% funtions
mpc=mpc_quadprog_handle(mpc);

options=optimset('Algorithm','interior-point','TolX',1e-15,'TolCon',1e-
15,'MaxFunEvals',10000);

% This function is the inbuilt matlab function for solving QP problem. The
% output is the sequence of desicion variables whose length depend on
% number of control horizons and number of inputs.
% u=quadprog(mpc.H_mat,mpc.F_vec,mpc.A_mat,mpc.b_vec,[],[],mpc.uk_l_mat,mpc
% .uk_h_mat,mpc.uk_guess',options);
%
u=quadprog(mpc.H_mat,mpc.F_vec,mpc.A_mat,mpc.b_vec,[],[],[],[],mpc.uk_guess
',options);
u=quadprog(mpc.H_mat,mpc.F_vec,[],[],[],[],[],[]);
% u
% pause
% Here we reconstruct the guess value for uk that is fed to the quadprog
n_ip=mpc.n_ip;
mpc.uk_guess=[mpc.uk_guess(n_ip+1:end) u(1:n_ip)'];

% This is the actual control input that is sent to the plant and the state
% estimator. From the sequence of u, we choose only the first set of imput.
uk=u(1:n_ip);
mpc.uk=uk;

3. mpc_mod_plant_mismatch- This function estimates the model plant

mismatch using different approaches like State Augementation and
Innovation Bias.
% This function returns the estimates of model-plant mismatch variables and
% constants like betakhat, etakhat, gamma_beta and C_eta;
function mpc=mpc_mod_plant_mismatch(mpc)

% Extraction of the fields from the mpc structure

yk=mpc.yk;
uk=mpc.uk;
xkhat_p=mpc.xkhat_p;
xkhat_c=mpc.xkhat_c;
betakhat=mpc.betakhat;
etakhat=mpc.etakhat;
phy=mpc.phy;
gamma_u=mpc.gamma_u;
L=mpc.L;
C_mat=mpc.C_mat;
R_mat=mpc.R_mat;
phy_beta=mpc.phy_beta;
phy_eta=mpc.phy_eta;
n_op=mpc.n_op;
MPM_type=mpc.MPM_type;
state_aug_type=mpc.state_aug_type;
Estimator_type=mpc.Estimator_type;
Q_beta=mpc.Q_beta;
Q_eta=mpc.Q_eta;

Q_mat=L*R_mat*L';
Pwv_mat=L*R_mat;

49 | P a g e
xkhat=xkhat_p;

xkhat_aug=[xkhat;betakhat;etakhat];
% Pkhat_aug=eye(size(phy_aug)); % whether to put this part as handle or
some big covariance here, ask sir

ek=yk-C_mat*xkhat_p;
if strcmp(MPM_type,'Innovation_bias')
betakhat=phy_beta*betakhat+(eye(n_op)-phy_beta)*ek;

if strcmp(Estimator_type,'Predictor')
etakhat=betakhat;

elseif strcmp(Estimator_type,'Current_state')
etakhat=phy_eta*etakhat+(eye(n_op)-phy_eta)*(yk-C_mat*xkhat_c);
end
gamma_beta=L;
C_eta=eye(n_op);

elseif strcmp(MPM_type,'state_augmentation')
if strcmp(state_aug_type,'input_bias')
gamma_beta=gamma_u;
C_eta=zeros(n_op);
elseif strcmp(state_aug_type,'output_bias')
gamma_beta=zeros(length(phy),n_op);
C_eta=eye(n_op);
elseif strcmp(state_aug_type,'disturbance_bias')
gamma_beta=gamma_d;
C_eta=zeros(n_op);
elseif strcmp(state_aug_type,'innovation_bias')
gamma_beta=L;
C_eta=eye(n_op);
end

phy_aug=daug(phy,eye(n_op),eye(n_op));
phy_aug(1:length(phy),length(phy)+1:length(phy)+n_op)=gamma_beta;
% phy_aug
% pause
gamma_u_aug=vertcat(gamma_u,zeros(n_op),zeros(n_op));
C_mat_aug=horzcat(C_mat,zeros(n_op),C_eta);
R1a_mat=daug(Q_mat,Q_beta,Q_eta);
R12a_mat=vertcat(Pwv_mat,zeros(n_op),zeros(n_op));
R2a_mat=R_mat;
Pkhat=mpc.Pkhat;
Pkhat_aug=Pkhat;

% Using the Kalman filter formulation as given in astrom and

wittenmark.
% Ask sir, how to incorporate the predictor here with wk and vk not
% independent

if strcmp(Estimator_type,'Current_state')

Kfk_mat=Pkhat_aug*C_mat_aug'*inv(C_mat_aug*Pkhat_aug*C_mat_aug'+R2a_mat);
Kvk_mat=R12a_mat*inv(C_mat_aug*Pkhat_aug*C_mat_aug'+R2a_mat);

Kk_mat=phy_aug*Kfk_mat+Kvk_mat;

50 | P a g e
xkhat_aug=phy_aug*xkhat_aug+gamma_u_aug*uk+Kk_mat*(yk-
C_mat_aug*xkhat_aug);
Pkhat_aug=Pkhat_aug-
Pkhat_aug*C_mat_aug'*(inv(C_mat_aug*Pkhat_aug*C_mat_aug'+R2a_mat))*C_mat_au
g*Pkhat_aug;

xkhat_aug=xkhat_aug+Kfk_mat*(yk-C_mat_aug*xkhat_aug);
Pkhat_aug=phy_aug*Pkhat_aug*phy_aug'+R1a_mat-
Kk_mat*(C_mat_aug*Pkhat_aug*C_mat_aug'+R2a_mat)*Kk_mat';

betakhat=xkhat_aug(length(phy)+1:length(phy)+n_op);
etakhat=xkhat_aug(length(phy)+n_op+1:length(phy)+2*n_op);
elseif strcmp(Estimator_type,'Predictor')

Kk_mat=(phy_aug*Pk_aug*C_mat_aug'+R12a_mat)*inv(R2a_mat+C_mat_aug*Pkhat_aug
*C_mat_aug');
xkhat_aug=phy_aug*xkhat_aug+gamma_u_aug*uk+Kk_mat*ek;
Pkhat_aug=phy_aug*Pkhat_aug*phy_aug'+R1a_mat-
Kk_mat*C_mat_aug*Pkhat_aug*phy_aug';
end

if strcmp(state_aug_type,'innovation_bias')
betakhat=betakhat+ek;
etakhat=etakhat+ek;
end
elseif strcmp(MPM_type,'DMC')
gamma_beta=0*L;
C_eta=eye(n_op);
if strcmp(Estimator_type,'Predictor')
etakhat=ek;

elseif strcmp(Estimator_type,'Current_state')
etakhat=yk-C_mat*xkhat_c;
end
end
mpc.gamma_beta=gamma_beta;
mpc.C_eta=C_eta;
mpc.betakhat=betakhat;
mpc.etakhat=etakhat;

4. mpc_setpoint_trajectory- This function is used to generate setpoint

trajectory
% This function returns the setpoint trajectory. Note that the trajectory
% can be filtered inside or outside. For filtering inside case, it is
% filtered in the function
function mpc=mpc_setpoint_trajectory(mpc)
phy_r=mpc.phy_r;
yk=mpc.yk;
rk=mpc.rk;
n_op=mpc.n_op;
% xr=mpc.xr;
xr=mpc.rk;
rk_f=mpc.rk_f;

51 | P a g e
% rk_f
% pause
setpoint_filter=mpc.setpoint_filter;

if strcmp(setpoint_filter,'filtering_inside')
xr=phy_r*xr+(eye(n_op)-phy_r)*(rk-yk);
yr=yk+xr;
rk_fs=yr;
elseif strcmp(setpoint_filter,'filtering_outside')
rk_fs=rk_f;

end

mpc.rk_fs=rk_fs;

5. mpc_steady_state_target- This function is used to calculate the terminal

target state.
% This function calculates the target state which is required for solving
% the infinite horizon MPC problem.
function mpc=mpc_steady_state_target(mpc)

% Extracting the relevant fields from the mpc structure

rk_fs=mpc.rk_fs;
C_mat=mpc.C_mat;
phy=mpc.phy;
gamma_u=mpc.gamma_u;
uk_l=mpc.uk_l;
uk_h=mpc.uk_h;
n_ip=mpc.n_ip;

% Invoking the unconstrained analytical solution

us=inv(C_mat*(inv(eye(length(phy))-phy)*gamma_u))*rk_fs;
xs=inv(eye(length(phy))-phy)*gamma_u*us;

% If us violates the constraints on uk, then the constrained optimization

% problem is invoked. A Quadratic programming problem is solved that tends
% to minimise the difference bewteen setpoint and the resulting model
% output.
for i=1:n_ip
if (us(i)<uk_l(i))||(us(i)>uk_h(i))
optimization='constrained';
else
optimization='';
end
end
us
pause

if strcmp(optimization,'constrained')
% Extracting the relevant fields needed for constrianed optimization
% problem
etakhat=mpc.etakhat;
C_eta=mpc.C_eta;

52 | P a g e
omega=mpc.omega;
n_ip=mpc.n_ip;
uk_guess=mpc.uk_guess(1:n_ip);

% Formulating the constraints on inputs

A_mat=[-eye(2);eye(2)];
b_vec=[-uk_l;uk_h];

% Constructing the Hessian and F

H_mat_steady_state=2*(gamma_u'*(inv(eye(length(phy))-
phy))'*C_mat'*omega*C_mat*(inv(eye(length(phy))-phy))*gamma_u);
F_vec_steady_state=-2*((rk_fs'-
(C_eta*etakhat)')*C_mat*(inv(eye(length(phy))-phy))*gamma_u);

% Calculating us and xs using the quadprog

us=quadprog(H_mat_steady_state,F_vec_steady_state,A_mat,b_vec,[],[],[],[],u
k_guess);
xs=inv(eye(length(phy))-phy)*gamma_u*us;
end

% The xs and us so obtained are included in the structure

mpc.xs=xs;
mpc.us=us;

6. mpc_block_mat_quadprog- This function builds the block matrices used

to calculate the handle of Quadprog.
% This function facilitates in building up the block matrices which are
% required in the formulation of quadratic programming problem.
function mpc=mpc_block_mat_quadprog(mpc)

% Extracting the parameters from the mpc structure

rk_fs=mpc.rk_fs;

xs=mpc.xs;
phy=mpc.phy;
gamma_u=mpc.gamma_u;
gamma_beta=mpc.gamma_beta;
C_mat=mpc.C_mat;
C_eta=mpc.C_eta;
p=mpc.p;
q=mpc.q;
n_ip=mpc.n_ip;
n_op=mpc.n_op;
we=mpc.we;
wu=mpc.wu;
wdelu=mpc.wdelu;
yk_l=mpc.yk_l;
yk_h=mpc.yk_h;
xk_l=mpc.xk_l;
xk_h=mpc.xk_h;

% Defining the empty block matrices which will later be filled

Rk=[];
Sx=[];
Su=[];

53 | P a g e
S_beta=[];
S_eta=[];
We_mat=[];
Wu_mat=[];
Wdelu_mat=[];
psi_mat=[];
psi_0_mat=[];
psi_u_mat=[];
psi_mat2=[];
Im_horz=[];
Im_vert=[];
psi_pq=[];
yk_lvec=[];
yk_hvec=[];
% This matrix fills the last row of C_eta to make it compatible to terminal
% state and hence infinite horizon formulation
C_eta_star=eye(length(phy),n_op);

% This matrix helps in the formation of some matrices

Im=eye(n_ip);
w_inf=dlyap(phy',C_mat'*we*C_mat);
% This loop returns the value of Su(later multiplied by
psi_pq),Sx,S_beta,S_eta,We_mat
for i=1:p

Suc=[];
for k=1:p
if(k<i)
Suc=[Suc;zeros(n_op,n_ip)];
elseif(k>=i)&&(k<p)
Suc=[Suc;C_mat*((phy)^(k-1))*gamma_u];
else
Suc=[Suc;(phy)^(k-1)*gamma_u];
end
end
Su=[Su Suc];

if(i<p)
yk_lvec=[yk_lvec;yk_l];
yk_hvec=[yk_hvec;yk_h];
Rk=[Rk;rk_fs];
Sx=[Sx;C_mat*(phy^i)];
S_eta=[S_eta;C_eta];
We_mat=blkdiag(We_mat,we);

else
yk_lvec=[yk_lvec;xk_l];
yk_hvec=[yk_hvec;xk_h];
Rk=[Rk;xs];
Sx=[Sx;phy^i];
S_eta=[S_eta;C_eta_star];
We_mat=blkdiag(We_mat,w_inf);
end

S_beta_row=zeros(size(phy));
for j=1:i
S_beta_row=S_beta_row+phy^(j-1);
end

54 | P a g e
if(i<p)
S_beta=[S_beta;C_mat*S_beta_row*gamma_beta];
else
S_beta=[S_beta;S_beta_row*gamma_beta];
end
end

for j=1:q-1
psi_mat2=blkdiag(psi_mat2,-Im);
Im_vert=vertcat(Im_vert,zeros(size(Im)));
end
psi_mat2=[psi_mat2 Im_vert];

for i=1:q
Wu_mat=blkdiag(Wu_mat,wu);
Wdelu_mat=blkdiag(Wdelu_mat,wdelu);
Im_horz=horzcat(Im_horz,zeros(size(Im)));

if(i<2)
psi_0_mat=vertcat(psi_0_mat,Im);
else
psi_0_mat=vertcat(psi_0_mat,zeros(size(Im)));
end

psi_mat=blkdiag(psi_mat,Im);
I_qm=eye(q*n_ip);
psi_u_mat=vertcat(psi_u_mat,Im);

Im_i=[];
for j=1:p/q % Note here that p/q should be an integer.
Im_i=[Im_i;Im];
end
psi_pq=blkdiag(psi_pq,Im_i);
end
psi_mat2=[Im_horz;psi_mat2];
psi_mat=psi_mat+psi_mat2;
Su=Su*psi_pq;

% size(I_qm)
% size(psi_u_mat)
% pause
% size(yk_lvec)
% pause

% Storing the values of all the block matrices in the mpc structure.
mpc.Rk=Rk;
mpc.Sx_mat=Sx;
mpc.Su_mat=Su;
mpc.S_beta_mat=S_beta;
mpc.S_eta_mat=S_eta;
mpc.We_mat=We_mat;
mpc.Wu_mat=Wu_mat;
mpc.Wdelu_mat=Wdelu_mat;
mpc.psi_mat=psi_mat;
mpc.psi_0_mat=psi_0_mat;
mpc.psi_u_mat=psi_u_mat;
mpc.I_qm_mat=I_qm;
mpc.yk_lvec=yk_lvec;
mpc.yk_hvec=yk_hvec;

55 | P a g e
7. mpc_quadprog_handle- This function calculates the Hessian, F and
constraints for the quadprog.
% This function returns the handles of quadprog like Hessian, F vector,
% A_mat, b_vec etc in terms of the block matrices defined in the previous
% funtions
function mpc=mpc_quadprog_handle(mpc)
Rk=mpc.Rk;
Sx_mat=mpc.Sx_mat;
Su_mat=mpc.Su_mat;
S_beta_mat=mpc.S_beta_mat;
S_eta_mat=mpc.S_eta_mat;
psi_mat=mpc.psi_mat;
psi_0_mat=mpc.psi_0_mat;
psi_u_mat=mpc.psi_u_mat;
I_qm_mat=mpc.I_qm_mat;
We_mat=mpc.We_mat;
Wu_mat=mpc.Wu_mat;
Wdelu_mat=mpc.Wdelu_mat;
deluk_l=mpc.deluk_l;
deluk_h=mpc.deluk_h;
xkhat=mpc.xkhat_p;
betakhat=mpc.betakhat;
etakhat=mpc.etakhat;
uk_prev=mpc.uk;
us=mpc.us;
uk_h=mpc.uk_h;
uk_l=mpc.uk_l;
yk_lvec=mpc.yk_lvec;
yk_hvec=mpc.yk_hvec;

% yk_l=mpc.yk_l;
% yk_h=mpc.yk_h;

% yk_lvec=[];
% yk_hvec=[];
% for i=1:p
% yk_lvec=[yk_lvec;yk_l];
% yk_hvec=[yk_hvec;yk_h];
% end
% Rk
% pause

% uk_guess=mpc.uk_guess;
zetak=(Rk-Sx_mat*xkhat-S_beta_mat*betakhat-S_eta_mat*etakhat);
H_mat=2*(Su_mat'*We_mat*Su_mat+psi_mat'*Wdelu_mat*psi_mat+Wu_mat);
F_vec=-
2*((Su_mat'*We_mat)*zetak+(psi_mat'*Wdelu_mat*psi_0_mat)*uk_prev+Wu_mat*psi
_u_mat*us);

% A_mat=vertcat(I_qm_mat,-I_qm_mat,psi_mat,-psi_mat,Su_mat,-Su_mat);
% b_vec=vertcat(psi_u_mat*uk_h,-
psi_u_mat*uk_l,psi_u_mat*deluk_h+psi_0_mat*uk_prev,-psi_u_mat*deluk_l-
psi_0_mat*uk_prev,yk_lvec-Sx_mat*xkhat-S_beta_mat*betakhat-

56 | P a g e
S_eta_mat*etakhat,Sx_mat*xkhat+S_beta_mat*betakhat+S_eta_mat*etakhat-
yk_hvec);
A_eins=vertcat(I_qm_mat,-I_qm_mat,psi_mat,-psi_mat);
b_eins=vertcat(psi_u_mat*uk_h,-
psi_u_mat*uk_l,psi_u_mat*deluk_h+psi_0_mat*uk_prev,-psi_u_mat*deluk_l-
psi_0_mat*uk_prev);
A_mat=A_eins;
b_vec=b_eins;
% A_mat=vertcat(I_qm_mat,-I_qm_mat,psi_mat,-psi_mat);
% b_vec=vertcat(psi_u_mat*uk_h,-
psi_u_mat*uk_l,psi_u_mat*deluk_h+psi_0_mat*uk_prev,-psi_u_mat*deluk_l-
psi_0_mat*uk_prev);

uk_l_vec=psi_u_mat*uk_l;
uk_h_vec=psi_u_mat*uk_h;

mpc.H_mat=H_mat;
mpc.F_vec=F_vec;
mpc.A_mat=A_mat;
mpc.b_vec=b_vec;
mpc.uk_l_vec=uk_l_vec;
mpc.uk_h_vec=uk_h_vec;

57 | P a g e
References

Adetola V Darryl D, Guay M. Adaptive Model Predictive Control for constrained non-linear
systems. Systems and Control Letters. 2009; 59:320-326

Astrom KJ and WittenMark B. Computer Controlled Systems, Prentice Hall, India, 1994.
Chapter-13, Identification.

Bemporad A, Trecate GF, Mignone D, Morari M, Torrisi FD. Model Predictive Control-Ideas
for the next generation. Automatisierungstechnik

Dixon R and Pike AW. Alstom Benchmark Challenge on Gasifier Control. IEE Proceedings-
Control Theory Applications. 2006; 153

Dixon R, Pike AW and Donne MS. Alstom Benchmark Challenge II on Gasifier Control.
Proceedings of the Institution of Mechanical Engineers, Part 1: Journal of Systems and
Control Engineering. 2000; 214:389

Johansson KH. The Quadruple Tank Process. A Multivariable Laboratory process with
Adjustable Zero. IEEE Transcations on Control Systems Technology.2000; 8.

Kadam JV , Schlegel M, Marquardt W, Tousain RL, Hessem DH, Berg J, Bosgra OH. A Two
Level Strategy of Integrated Dynamic Optimization and Control of Industrial Processes- a
case study. European Symposium on Computer Aided Process Engineering; 2002, 511-516

Kadam JV, Marquardt W, Schlegel M, Backx T, Bosgra OH, Brouwer PJ, Dunnebier G,
Hessem Dv, Tiagunov A, Wolf Sd. Proceedings Foundation of Computer –Aided Process
Operations. 2003

Morari M and Lee JH. Model Predictive Control: past, present and future. Computer and
Chemical Engineering. 1999; 23:667-682

Muske KR and Rawlings JB. Model Predictive Control with Linear Models. AIChE Journal.
1993; 39.

Patwardhan SC, Lecture Notes on Advanced Process Control. IIT Bombay.

Qin SJ, Badgwell TA. A survey of industrial Model Predictive Control Technology. Control
Engineering Practice. 2013; 11: 733-764

Tanaskovic M, Fagiano L, Smith R, Goulart P, Morari M. Adaptive Model Predictive Control

for Constrained Linear Systems. European Control Conference (ECC). 2013

58 | P a g e