ROBUST MODEL PREDICTIVE CONTROL

UNDER OUTPUT CONSTRAINTS

A.Brdys T.Chang

School of Electronic, Electrical Engineering and Computer, The University of Birmingham, Birmingham B15
2TT, UK; Tel: +44 121 414 4354, Fax: +44 121 414 4291, email: [email protected]

Abstract: Control of a linear time-varying uncertain dynamical system with delayed

inputs is considered. The model parameters, disturbance inputs and model structure errors
are unknown but bounded, and the parameter value can abruptly change. The objective is
to keep the system output within prescribed limits regardless of the uncertainty scenario.
A model predictive type of controller is designed that utilises a set bounded model of the
uncertainty and employs safety zones modifying the original constraints so that the
control input feasibility can be guaranteed. The controller is applied to quality control in a
benchmark Drinking Water Distribution System (DWDS), and its performance is
validated by simulation. Copyright © 2002 IFAC

Keywords: transport delay, uncertain dynamic systems, predictive control, bounding

method, robust control, constraints.

1. INTRODUCTION that uses a set bounded state estimation and an extra

constraint added to the optimisation problem in order
Model Predictive Control (MPC) is an effective tool to improve the computational efficiency. In the paper
to deal with time-delay and hard constraints in a of Gossner (1997) and Chisci (2001) satisfied results
systematical way. The issues of optimisation have been obtained about constraints fulfilment,
feasibility, stability and performance for linear stability and state convergence under bounded
systems have been investigated to a sufficient disturbance.
maturity in the literatures (Mayne, 2000). The In practical the input constraint was included in the
situation gets highly complicated when the output optimisation problem explicitly, and can be solved as
constraints are present and there are model-reality hard constraints. In this paper, as in (Brdys, et al,
differences. Several properties of MPC have to be 1998), a method of fulfilling output constraint by
rechecked, such as the performance, stability and introducing so called safety zones into the original
constraint fulfilling. An overview of robustness in constraints under the uncertainty of model,
MPC was given in (Beporad and Morari, 1999). The measurement and actuator errors and modelling
feasible controller action may violate the constraints errors is proposed. For the modified set of constraints
when applied to the physical system. Much of the a feasible region in a control space is shrank
work on the MPC that was carried out in rent years comparing to the original one. Hence, with suitably
addressed this issue. A linear state feedback based chosen safety zones the model feasible control input
MPC was investigated in (Kothare, 1996; also satisfies constraints of a physical system if the
Kouvaritakis, 1999; Scokaert, 1998) employing a uncertainty radius is not too large. Given proposed
concept of invariant set. A feedback MPC was control action and an estimate of a system state at a
introduced in (Scokaert, 1998) in order to reduce present time instant the system output response is
conservatism of the controller, hence allowing predicted over a prediction horizon based on the
handling larger uncertainty radius. Manipulating the model of system dynamics and uncertainty. The
reference trajectory was efficiently used in (Gilbert, safety zones are validated and redesigned at each step
1995; Bemporad and Mosca, 1998; Angeli, et al, of MPC, if necessary, by applying a robust output
2001) in a case of unknown disturbances. In a recent prediction based on a set-membership bounded
papers (Bemporad and Garulli, 2000) a linear system uncertainty model (Milanese, et al, 1996; Brdys,
with unmeasured state and bounded disturbance 1999).
input and under input and state constraints was
investigated. Feasible MPC algorithm was developed
The paper is organised as follows. The problem is
formulated in section 2. The robust predictive θ (t ) ∈ Θ(t ) (3)
Θ(t ) ={θ ∈ R n : θ l (t ) ≤ θ ≤ θ u (t )}
controller is performed in section 3, and its design is ∆

described in section 4. The controller is applied to

quality control in a benchmark Drinking Water n = q + iu − il + 1
Distribution System (DWDS), and its performance is θ l (t ) = [b1l (t )

bql (t ) a il (t )

a il (t )]
l u
validated by simulation. The results are presented in
section 5. The quality control in DWDS is at present θ u (t ) = [b1u (t )

bqu (t ) a iu (t )
l

a iu (t )]
u

extremely important application field still waiting for ε 0 (t ) ∈ Ε 0 (t ) (4)

practically sound solution (Polycarpou, et al, 2001).
Ε 0 (t ) ={ε 0 ∈ R : ε 0l (t ) ≤ ε 0 ≤ ε 0u (t )}
∆

2. PROBLEM FORMULATION The envelopes known at t 0 covers the period over

2.1. Time-varying Bounded Model and control horizon and the uncertainty radius must be
Control Objective sufficiently small so that the control objectives can
be robustly achieved. Although there is a modelling
The discrete time-varying SISO system is described error, there is no one set of true parameter values
by the input-output model: matching the model and plant for all inputs. The
model structure reflects well real plant dynamics so
q i = iu
y (t ) = ∑ bi (t ) y (t − i ) + ∑ a i (t )u (t − i ) + ε 0 (t ) (1) that the parameter and modelling error bounds are
i =1 i = il reasonable.
The control aims at keeping the plant output y p (⋅)
where y (t ) , u (t ) , bi (t ) , a i (t ) and ε 0 (t ) ∈ R , t ∈ Z ,
within the output constraints described by the lower-
y (t ) is the output, u (t ) is the plant input, bi (t ) upper bounds:
and a i (t ) are time-varying model parameters, ε 0 (t )
is the composite error of modelling and disturbance y min (t ) ≤ y p (t ) ≤ y max (t ) (5)

over the control horizon t ∈ [t 0 , t 0 + Tc ] , and there are

input, the integer i l and i u define range of delays in
the control input, and Z denotes set of integers.
Let constraints on the control input to be satisfied:

θ (t ) = [θ 1 (t ) θ n (t ) ] u min (t ) ≤ u (t ) ≤ u max (t )
(6)
= [b1 (t ) bq (t ) a i (t ) a i (t ) ] | u (t ) − u (t − 1) |≤ ∆u max
l u

n = q + iu − il + 1 Let us denote by Ω the set of all feasible control

ϕ (t − 1) = [ y (t − 1) y (t − q) u (t − i l ) u (t − i u )]
T

trajectories, that is:

∆
Then equation (1) can be compactly written as: Ω ={u (⋅) : u (t ) ∈ R, c(u (t ), y p (u (⋅))(t )) ≤ 0,
y (t ) = θ (t )ϕ (t − 1) + ε 0 (t ) (2) for t ∈ [t 0 , t 0 + Tc ]} (7)

The θ (⋅) and ε 0 (⋅) constitute the model uncertainty. where c(⋅) is the input and output constraint
conditions represented by (5) and (6), and
Let y [pt , t ] (⋅) denote the plant response to input
0
y p (u (⋅))(t ) is the plant response to the input u (⋅) .
u[ t ,t ] (⋅) over time period [t 0 , t ] . Clearly there exist
0 As the control problem is under constraints and the
the trajectories of θ (⋅) and ε 0 (⋅) so that with these inputs are delayed the Model Predictive Control
scenarios of the model uncertainty the model (MPC) will be applied to design the controller.
response equals to the plant response, that is Hence, it is assumed that Tc < Tm .
y [ t ,t ] (⋅) = y [pt ,t ] (⋅) . Different inputs require different
0 0

scenarios of θ (⋅) and ε 0 (⋅) in order to match the 2.2 Parameter Piece-Wise Constant Bounded Model
plant responses. It is assumed that the control input is
valued on a compact set so that the trajectories of Parameter θ (t ) can abruptly change their values at
θ (⋅) and ε 0 (⋅) can be bounded above and below time instant t j ∈ [t0 , t0 + Tm ] , j = 1,2, 

Np and
over a time interval [t 0 , t 0 + Tm ] by the bounded t N = t0 + Tm . Within the time interval [t j −1 , t j ] the
p

envelopes θ (⋅) , θ (⋅) and ε (⋅) , ε (⋅)

u,0 l ,0 u ,0
0
l ,0
0 changes can be assumed slow and also, only certain
respectively, where t 0 is the initial time and Tm is parameters are active, that is their values are nonzero.
the modelling horizon. The tightest envelopes are not The instants t j are assumed known. Under these
known. It is assumed that at t 0 a priori envelopes assumptions the whole model horizon Tm can be
θ (⋅) , θ (⋅) and ε (⋅) , ε (⋅) sufficiently well
u l u
0
l
0
partitioned into N p time slots, defined as:
approximate the unknown envelopes so that no
S j ={t ∈ Z : t j −1 ≤ t < t j }, j = 1,2,
∆
bounding modifications will be considered during 

, Np (8)
control design. Hence the following holds:
 where H m is the control sequence dimension, H p is
Let I j define set of active parameter indices over
the prediction horizon, H m ≤ H p , and Q is positive-
S j . Hence, θ j = {θ i }i j j ∈I j
is the active parameter definite matrix. The first part of the index is to
sub-vector over S j . Hence, θ (t ) = θ j
for t ∈ S j , minimise the total control energy and at the same
time attempts to avoid abrupt control value change.
and θ is treated as constant. The corresponding
j
For the water quality control by chlorine injection,
bounds on θ j are calculated as follows: this part is to minimise the total chlorine dose

( )
applied. The terminal output of the finite predictive
θ min
j
= {min
t∈S
θ il (t ) }i ∈I j j j
horizon is penalised by the second part of the index,
= {max (θ (t ) )}
j
which tends to keep the terminal output close to the
θ j
max
u
ij i j∈I j reference value y r by applying a tuning knob µ .
t∈S j

When apply the nominal model at k=1 all the past

Hence, the parameter bounding orthotope Θ j over outputs and inputs required are taken from the
S j can be calculated as: measurements. When k increases the unavailable
measurements are replaced by the model responses.
The delays imply that at least i l < H m . However,
{ }
∆
Θj =θ j ∈R : θ min ≤ θ j ≤ θ max ,t ∈S j
dim( I j ) j j
(9)
preferably i u < H m . Hence, availability of
sufficiently tight uncertainty bounds needed for
The uncertainty in ε 0 (⋅) is tackled similarly.
achieving the control objectives and the delay range
Finally, the model with piecewise constant parameter may be faced.
is given as: Notice that the deterministic nominal model allows
to quickly generate control sequence that is optimal
Model ℜ(⋅): y (t ) = θ (t )ϕ (t − 1) + ε 0 (t ) (10) for the selected scenario of the uncertainty. However,
θ (t ) = θ , ε 0 (t ) = ε , for t ∈ S j
j
0
j
(11) it is suboptimal, if feasible, for the real plant. The
optimality robustness can be improved by
θ ∈Θ j j
ε ∈Ej j
(12)
0
formulating the optimisation problem as the minmax
one. This however is not further pursued in the paper
It is understood that the vector ϕ (⋅) changes its as solving constrained minmax problem would
structure following changes of the parameter vector. greatly increase the computational burden even for
Model ℜ(⋅) of (10)-(12) is called the parameter linear-quadratic problem at hand. The feasibility
piece-wise constant bounding model that will be problem however, needs to be addressed.
applied for MPC design. Piecewise constant model is Minimisation of the performance index giving by
preferred for practical reason of the computational (14) over the constraints described by (7) where the
and model identification efficiency. nominal model is used to replace unknown plant
mapping y p (⋅) yields the solution that may not be
3. CONTROLLER STRUCTURE feasible when applied to the plant. In order to
guarantee the feasibility the nominal model based
3.1 MPC Based on Nominal Model and Modified constraints are modified by introducing so called
Constraints safety zones. The modified constraints define
narrower set for the control actions and this is the
The nominal model N(⋅) is defined by the nominal price to be paid for model-reality differences. The
scenario of uncertainty. Let us take some kind of safety zones can be viewed as the feasibility margins
centre of Θ j as the nominal parameter value for that accommodate the uncertainty. It will be
technically proceeded as follows.
t ∈ S j , denoted as θˆ j . For example:
The output constraints over the prediction horizon
described by the upper and lower limits are:
θˆ j = (θ min
j
+ θ max
j
)/2 (13)
[
Y min = y min (t + 1) 

y min (t + H p ) ]T
(16)
= [y )]
Hence, the nominal parameter trajectories satisfy
(t + 1) (t + H p
max max max T
Y 

y (17)
θˆ(t ) = θˆ j for t ∈ S j . It is assumed without any loss
of generality that εˆ 0j = 0 . The optimisation problem The limits are modified by the safety zones σ l and
to be solved by MPC at time instant t can be σ u as:
formulated as:
Ysmin = Y min + σ l Ysmax = Y max − σ u
J (Uˆ ) = Uˆ T QUˆ + µ[ y (t + H p | t ) − y r ] 2 (14)
σ l = [σ 1l σ Hl ]T σ u = [σ 1u σ Hu ] T
[ ]



p p

Yˆ = y (t + 1 | t ) y (t + H p | t )
T

[σ σ ]∈ ∑


l u

Uˆ = [u (t | t ) u (t + H m − 1 | t )]
T


∆
∑ ={[σ σ u ] : σ il ≥ 0, σ iu ≥ 0, for i = 1
l 

H p and
y (t + k | t ) = θˆ(t + k )ϕ (t + k − 1) , for k = 1, ,Hp
(t + i ) + σ < y (t + i ) − σ }


min l max u
y i i (18)
(15)
The modified model based constraints read:
 | u c (t ) − u (t ) |≤ ε e
Ys min
≤ Yˆ ≤ Ysmax (19)
where y m (t ) is the plant output measurement, u c (t )
The condition (18) is needed in order to make sure
is the controller output or the actuator input and
that the modified constraints define nonempty set.
The input constraints are treated as hard constraints
ε m and ε e are the error bounds.
in the optimisation formulation. The actuator error The plant output prediction at t over [t , t + H p ] is
ε e will also be considered in implementing a control performed based on a priori information contained in
command. Hence, the input constraints are modified the past inputs and output measurements, future
in order to cater for the actuator error as: inputs, input-output model equations and uncertainty
bounds. This a priori information has been described
U min ≤ Uˆ ≤ U max (20) in a form of equalities and inequalities constraining
outputs over [t , t + H p ] . Any output trajectory
− ∆U max ≤ ∆Uˆ ≤ ∆U max
U min = [u min (t ) + ε e u min (t + H m − 1) + ε e ]

T satisfying these constraints can be the plant response.
The robust output prediction provides the intervals:
U max = [u max (t ) − ε e u max (t + H m − 1) − ε e ]

T

Ypl = [ y lp (t + 1 | t ) y lp (t + H p | t )]T
∆U max = [∆u min (t ) − 2ε e ∆u min (t + H m − 1) − 2ε e ]



T

Ypu = [ y up (t + 1 | t ) 

y up (t + H p | t )]T
Finally, the MPC optimisation task at t reads:
over [t , t + H p ] bounding the plant output values
arg min ( J (Uˆ )) over [t , t + H p ] . Hence,
Uˆ (σ l , σ u ) = Uˆ

subject to Uˆ ∈ Ω s
(21)
y lp (t + k | t ) ≤ y p (t + k ) ≤ y up (t + k | t ) (23)
∆
 N (•)
Ω s =Uˆ ∈ R H m : c s (Uˆ , Yˆ ) ≤ 0, Yˆ = Uˆ → Yˆ  for k = 1, ,Hp
î 


where c s (⋅) is the mapping describing the modified Based on model ℜ(⋅) of (10)-(12), the plant output
input and output constraints (20) and (19). After the bounding constraints at t+k can be summarised as:
optimisation problem has been solved, only the first
m = 1, , k 

value of Û is to be applied as the current controller

y (t + m | t ) − θ (t + m)ϕ (t + m − 1) ∈ Ε(t + m)
output at time instant t :
θ (t + m ) ∈ Θ (t + m )
u c (t ) = u (t | t ) | y m (t + k − m) − y p (t + k − m) |≤ ε m , for k − m ≤ 0,
| u (t + k − m) − u (t + k − m) |≤ ε e
c

At the next time instant t + 1 , the predictive horizon

moves forward and the minimisation process of
where E (t ) = E j , Θ(t ) = Θ j for t ∈ S j . Let P k
solving Û is repeated and thus is a receding horizon
predictive control scenario. denote the set of all y (t + k | t ) satisfying the above
If Uˆ (σ l , σ u ) is feasible for the plant then u (t | t ) is constraints. Hence, the k − th step robust output
applied. The feasibility is assessed by performing at t prediction at time instant t can be defined as:

min [y (t + k | t )]
a robust prediction of the plant response to
Uˆ (σ l , σ u ) . If the feasibility of the proposed control y lp (t + k | t ) = (24)
subject to y (t + k | t ) ∈ P k
can not be guaranteed then the safety zones σ , σ l u
max[y (t + k | t )]
are redesigned and new control actions are generated y up (t + k | t ) = (25)
as before. This repeats till suitable safety zones are subject to y (t + k | t ) ∈ P k
determined.
If the predicted plant output satisfy:

3.2 Assessment of Feasibility by Robust Prediction Ypl ≥ Y min and Ypu ≤ Y max (26)

Set bounded model of uncertainty that is used in this

then clearly, the assessed control sequences Û is
paper enable us to calculate upper and lower
guaranteed to be feasible. In other words it is
envelopes bounding real plant response to a specific
robustly feasible.
input. Comparing these envelopes against the bounds
defining the plant output constraints allows to assess
the input feasibility. An algorithm for the envelope 3.3 Operation of Robust MPC Controller
calculating shall now be presented.
The errors in the plant output measurements and in In general, determining suitable safety zones requires
executions of the control inputs (actuator error) are a number of iterations to be performed. The
bounded as: algorithm for iterative calculation of the safety zone
will be described in the next section. The overall
| y m (t ) − y p (t ) |≤ ε m (22)
robust MPC controller is of iterative type and it exists such value Λ* of the multiplier Λ that the
operates as follows: safety zones can be calculated by solving the
following problem (Fletcher, 1987):
(i) Let [σ l σ u ] = 0 , solve Û using (21);
(ii) Calculate Y pl , Y pu using (24),(25) respectively; arg min φ (σ l , σ u , Ψ , Λ* )
[σ l σ u]= (29)
If (26) is satisfied go to (iv) [σ l σ u ]
Else go to (iii) subject to [σ l σ u ] ∈ ∑
(iii) Redesign [σ l σ u ] , and calculate Û based on
this safety-zone design, then go to (ii) If the optimum value of the penalty function is equal
to zero then the wanted safety zones are found. If it is
(iv) Let u c (t ) = u (t | t ) nonzero then it means that the MPC controller is not
able to produce control that feasibility can be
robustly assessed. In this situation, assuming good
4. ROBUST PREDICTIVE CONTROLLER controllability of the plant, the uncertainty radius
DESIGN must be reduced in order to regain robustly feasible
operation of the MPC controller. Noticing that
In the previous section a structure of robust MPC solution for (29) is not unique, finding the smallest
controller has been proposed. An implementation of safety zones remains an open problem. Following
this structure needs dedicated algorithms for solving (Powell, 1969; Hestenes, 1969) an algorithm for
variety of problems. Firstly, calculating Uˆ (σ l , σ u ) simultaneous solving (29) and finding Λ* shall be
requires solving constrained linear quadratic derived. A key assumption needed is second order
Fr chet differentiability of the mapping C (⋅,⋅) . This

optimisation problem. A number of efficient solvers

exist to perform this task. Secondly, performing the mapping is a composition of the mappings describing
robust plant output prediction requires solving the the nominal model based MPC generation (21),
nonlinear and non-convex optimisation problems robust plant output prediction (24-25) and f (⋅) (28).
(24) and (25). An approximated solving approach is Clearly, the later one is smooth. The former ones are
proposed by piece-wise linearisation, with the defined by a solution of a constrained optimisation
linearisation error included in the modelling error.
problem parameterised by σ l , σ u and by Uˆ (σ l , σ u )
The final problem to be solved is a linear mixed-
integer programming (MIP) problem (Brdys, 1999), respectively. As the parameters enter the constraints
and it can be solved by using a standard solver. the elegant sufficient conditions for differentiability
Thirdly, in order to calculate suitable safety zones a of the solution do not exist (Hager, 1979; Clarke,
penalty function related to the constraints (26) is 1983). However, for broad class of problems the
employed. Let us define: differentiability holds (Findeisen, et al, 1980). Hence,
an existence of all derivatives needed is assumed in
C (σ l , σ u ) = [ f (V1 ) 

f (V 2 H )] T
p
this sequel. The algorithm I can now be stated as:
(27) Algorithm I:
where
V = [V1 

[
V 2 H ]T = (Y min − Y pl ) T
p
(Y pu − Y max ) T ]
T
(i) Initially set k = 0, Λ = Λ(1) , Ψ = Ψ (1) ,
and || C ( 0 ) || ∞ = ∞ ;
 x if x > ε+ (ii) Find the minimizer of (29) x(Ψ , Λ ) and denote
 C = C (x) ;
f ( x) = ( x + ε ) / 4ε +
+ 2
if −ε+ ≤ x ≤ ε+ (28)
 0 if x < −ε +
1 1
î (iii) If || C || ∞ > || C (k ) || ∞ , ∀i :| C i |> || C (k ) || ∞ set
4 4
ψ i = 10ψ i , and go to (ii);
where ε + is small positive number. Notice, that if (iv) Set k = k + 1, Λ( k ) = Λ, Ψ ( k ) = Ψ , C ( k ) = C ;
C (σ l , σ u ) = 0 holds then (26) holds as well. (v) Set Λ = Λ( k ) − Ψ ( k ) C ( k ) and go to (ii) until the
Moreover, the constraint C (σ l , σ u ) = 0 can get constraint is fulfilled with desired accuracy,
arbitrarily close to (26) by setting ε + sufficiently where x = [σ l , σ u ] .
small. The multiplier (shifted) type of penalty
function associated with (27) is defined as (Fletcher, The algorithm can ensure global convergence with
1987): the convergent rate of 0.25 by online tuning in step
(iii).
1 If Λ* is known, (29) can be attempted by Newton
φ (σ l , σ u , Ψ , Λ ) = [C (σ l , σ u )] T Ψ[C (σ l , σ u )] methods. Applying the Newton algorithm we obtain
2
that:
− ΛC (σ l , σ u )
∇φ ( x) = AΨC ( x) − AΛ* where A = ∇C (x )
where Ψ = diag 2 H ψ i , Λ = [λ1 

λ 2 H ] , λi are
∇ 2φ ( x ) = ∇ 2 C ( x)[ΨC ( x) − Λ* ] + AΨA T = W x
p p

the multipliers. Under rather mild conditions there

For large values of ψ i , the following approximation the node 5 mainly controls node 16. There are
holds (Fletcher, 1987): interactions between the two control loops. Hence it
is a 2-input 2-output system. The chlorine
[ A T W x−1 A] −1 ≈ Ψ (30) concentrations at the monitored nodes should be
maintained within 0.20[mg/l]-0.30[mg/l] limits.
Hence, assuming the inverse of A exist, Hence, y max = 0.30[mg / l ] and y min = 0.20[mg / l ] .
The 2-input 2-output controller for the booster
W x−1 ≈ ( A T ) −1 Ψ −1 A −1 stations is to be designed.

and Newton method of solving (29) yields: 5 7

201 201 200

x ( k +1) = x ( k ) − ( A T ) −1 [C ( x ( k ) ) − Ψ −1 Λ* ] (31)
6 18 13 14

It is a special property of our problem that A is 12

6 13 25 14 16

square matrix so that existence of A −1 is not unusual.

8
4 15 26

As Λ* does not depend on Ψ then Ψ −1 Λ* ≈ 0 for 27 15

large values of ψ i . Finally, a greatly simplified 5

11 24

algorithm is obtained as:

7 8 19
23
9
x ( k +1)
=x (k )
− ( A ) C(x )
T −1 (k)
(32) 17
4 3
16

22
Notice that the iteration (32) does not require Λ* any
longer. It is surprising at the first glance. However, 10
9 21 12
the formula (32) has a very appealing form. Namely, 3

new safety zones are calculated by a correction of the 18 20

present ones using the extend of constraint violation. 10 19 11

As calculating ( A T ) −1 may be still computationally 100

101 101 1 1 2 2
17

demanding further simplification is proposed that

consist in replacing ( A T ) −1 by constant scalar gain. Fig. 1 A Drinking Water Distribution Network
The resulting algorithm becomes of standard
relaxation type: The pump schedule, flows in the pipes and tank
levels are obtained from the upper level controller
Algorithm II: optimising in an integrated manner the quality and
quantity (Brdys, et al, 2001), and the demand
(i) Set x = [σ l σ u] = 0 ; prediction is assumed known precisely. The same as
for the quantity control a 24-hour control horizon is
(ii) Solve Û using (21), if C ( x) = 0 is satisfied
considered. The sampling period is 5 minutes
then go to (iv); yielding totally 288 steps for the whole control
(iii) Using x ( k +1) = x ( k ) − νC ( x ( k ) ) find new safety- horizon.
zone [σ l σ u ] , go to (ii); The water network was implemented using
EPANET2.0(Rossman, 2000) in order to simulate the
(iv) Set u c (t ) = u (t | t ) . plant responses. The chlorine transportation from the
inject node to the monitored node without tanks in
where v is the relaxation gain and its possible choice the DWDS can be described by a input-output model
is v = max(diag[∇C (0)]) . presented in (Zierolf, 1998):

5. SIMULATION RESULT y (t ) = ∑ α p (t )u (t − d p (t )) (33)

p ∈P ( t )

A drinking water distribution case-study system was where P(t ) is the set of the paths from the injection
investigated. The network structure is illustrated in
Fig. 1. There are 16 nodes, 27 pipes and 3 tanks in node to the monitored node, α p (t ) is the chlorine
the system. The water is pumped from the source decay factor related to path p and d p (t ) is the
(node 100 and node 200) and is also supplied by the detention time in path p . It describes a time-varying
tanks (node 17,18,19). Node 16 and node 8 are
system with time-varying delayed inputs. In order to
selected as monitored nodes as they are the most
make the model applicable for control and
remote nodes from the source. The chlorine
concentrations at these nodes are the two plant estimation, the delayed time d p (t ) is discretized as
outputs y1 and y 2 . In order to maintain the chlorine in (Polycarpou, 2001). Finally, a time-varying model
with constant delays is obtained and extended to a
level at the monitored nodes, two booster stations are
MIMO form including tanks:
installed at the nodes 5 and 10. The chlorine
concentrations at the injection nodes are the inputs 2 2 i = 24

u 1 and u 2 . The node 10 mainly controls node 8 and y n (t ) = ∑ b ni (t ) y n (t − i ) + ∑ ∑ a m , n , i (t )u m (t − i )

i =1 m =1 i=6
where n = 1,2 , am.n.i describes an impact of the chlorine at the booster station nodes. The simulation
results have illustrated good performance of the
injection input m that is delayed by i steps on the controller.
output n. Notice that the time-varying delays are
replaced by 19 constant delays from range of 6-24,
yielding 76 parameters. As described in section 2 the ACKNOWLEDGEMENT
whole control horizon is partitioned into eight time This work was partly supported by the Polish State
slots generating a piece-wise constant parameter Committee for Scientific Research under grant No.8
model, and only certain parameters are active within T11A 022 16.
certain time slots. Fig.2 shows the envelopes
bounding the parameters a1,1, 7 over a whole horizon
of 288 steps. Notice that the parameters are active,
that is they are nonzero, only over certain time REFERENCES
periods. The centre values of these parameter
envelopes are taken as the parameter values in the Angeli, D., A. Casavola and E. Mosca (2001). On
nominal model. With these parameter estimation feasible set-membership state estimators in
results, the modelling error was ±4% of the plant constrained command governor control.
output value. Automatica, 37, 151-156.
In the following presentation of the simulation Bemporad, A., and E. Mosca (1998). Fulfilling Hard
results, the chlorine concentration was scaled by the Constraints in Uncertain Linear Systems by
factor of 0.25mg/l, so the upper and lower output Reference Managing. Automatica, 34, 451-461.
limit (0.2mg/l and 0.3mg/l) are converted into 0.8 Bemporad, A. and M. Morari (1999). Robust Model
and 1.2 respectively, and the output reference is Predictive Control: A Survey. In: Robustness in
y r = 1.0 . The measurement error and actuator error Identification and Control (A. Garulli, A. Tesi
were ±2.5% of the measurement value and and A. Vicino (Ed)), 207-226. Springer, London.
controller output value respectively, and the Bemporad, A., and A. Garulli (2000). Output-
corresponding error bounds of ε m , ε e were obtained. feedback predictive control of constrained linear
systems via set-membership state estimation.
The MPC described in section 2 was designed using Int. J. Control, 73, 655-665.
algorithm II to calculate the safety zones, where Brdys, M.A., J.T. Dud and P. Tatjewski (1998).
H m = H p = 36 . The controller starts with zero Improving optimality in multilayer control
safety zones. Its operation over a whole horizon and systems by tighter constraint control and
the output constraint violation are illustrated in Fig. supervision. Proc. of the 8th IFAC/IFORS/
3. The violation is about 5% above the output limit IMACS/IFIP Symposium on Large Scale
around steps t=140 and t=230. The operation of the Systems: Theoty & Applications.Vol.1, 468-473,
controller over the same time period but with the Paris.
modified output constrains by safety zones is shown Brdys, M.A. (1999). Robust Estimation of Variables
in Fig.5, hence achieving the feasibility. The safety and Parameters in Dynamic Networks. Proc. of
zones generated at step t=204 are illustrated in Fig. the 14th IFAC World Congress, Beijing,
4. The relaxation gain used was υ = 1.0 . P.R.China, 5-9 July, 1999.
Brdys, M.A., T. Chang, and K. Duzinkiewicz (2001).
Intelligent Model Predictive Control of Chlorine
6. CONCLUSION Residuals in Water Distribution Systems. Proc.
of World Water & Environmental Resources
A robust MPC controller has been developed for
Congress, 20-24 May, Orlando, Florida.
keeping an output of a linear time varying systems
Chisci, L., J.A. Rossiter and G. Zappa (2001).
under uncertainty within prescribed limits with
Systems with persistent disturbances: predictive
delayed inputs. The uncertainty in: time varying
control with restricted constraints. Automatica,
parameters, measurement errors, actuator errors and
37, 1019-1028.
modelling error has been modelled applying set-
Clarke, F.H. (1983). Optimization and nonsmooth
bounded models. The safety zones have been
analysis. Wiley, New York.
introduced to modify the model-based output
Findeisen, W., F.N. Bailey, M. Brdys, K.
constraints so that control input feasibility can be
Malinowski, P. Tatjewski, A. Wozniak (1980).
robustly achieved. The safety zones have been
Control and Coordination in Hierarchical
iteratively designed at each generation time instant of
Systems. John Wiley & Sons, Chichester-New
the MPC based on the envelopes bounding the
York.
predicted plant output responses. Algorithms for
Fletcher, R. (1987). Practical Methods of
generating the safety zones based on the constraint
Optimisation. John Wiley & Sons, Chichester-
violations extend over the MPC prediction horizon
New York.
have been derived and their convergence has been
Gibert, E., I. Kolmannowski (1995). Discrte-time
analysed. An efficient simple relaxation scheme has
reference governors for systems with state and
been designed to reduce the computational burden.
control constraints and disturbance inputs. Proc.
The controller has been applied to a drinking water
Of the 34th IEEE Conference on Decision and
distribution system to maintain chlorine
Control, 1189-1194.
concentration at a monitored demand node within
prescribed limits by controlling injection of the
Gossner, J.R., B. Kouvaritakis and J.A. Rossiter Output of Node 8
(1997). Stable generalized predictive control 1.3

1.2
with constraints and bounded disturbances.
1.1
Automatica, 33, 551-568.
1
Hager, W.W. (1979). Lipschitz continuity for
0.9
constrained processes. SIAM J. Control and 0.8
Optimisation, 17, no.3, 321-338. 0.7
0 50 100 150 200 250 300
Kothare, M.V., V. Balakrishnan and M. Morari
(1996). Robust constrained model predictive Chlorine Injection at Node 1
2.2
control using linear matrix inequalities.
Automatica, 32(10), 1361-1379. 2.1

Kouvaritakis, B., Rossiter, J.A., and Schuurmans, J. 2

(1999). Efficient robust predictive control.

1.9
Proc. Of the American Control Conference,
4283-4287, San Diego, California. 1.8
0 50 100 150 200 250 300
Hestenes, M.R. (1969). Multiplier and gradient Fig.3 Controller Performance with Zero Safety Zones
methods, J. Opt. Theo. Applns., 4 303-320. (Dashed Line: Upper Limit
Mayne, D.Q., J.B. Rawlings, C.V. Rao and P.O.M. Dash-Dot Line: Lower Limit)
Scokaert (2000). Constrained model predictive
control: Stability and optimality. Automatica,
36, 789-814. 1.3

Milanese, M., J. Norton, H. Piet-Lahanier and E.

Walter (Eds) (1996). Bounding Approaches to 1.2

System Identification. Penum, New York.

Polycarpou, M.M., J.G. Uber, Z. Wang, F. Shang, 1.1
M.A. Brdys (2001). Feedback control of water
quality. IEEE Control Systems Magazine.
1
(Accepted for Publication).
Powell, M.J.D. (1969). A method for nonlinear
constraints in minimization problems, in 0.9

Optimization (Ed. R. Fletcher). Academic Press,

London. 0.8

Rossman, L.A. (2000). EPANET 2.0 for Windows.

Water Suply and Water Resourses Division, 0.7
0 5 10 15 20 25 30 35
National Risk Management Research
Laboratory, Cincinnati, OH 45268. Fig.4 An Example of Safety Zone Design
Scokaert, P.O.M., D.Q. Mayne (1998). Min-Max (Solid Line: Modified Output Constraints
feedback model predictive control for Dashed Line: Robust Output Prediction Envelope
constrained linear systems. IEEE Transactions Dash-Dot Line: Output Constraints)
on Automatic Control, 43(8), 1136-1142.
Zierolf, M.L., M.M. Polycarpou, J.G. Uber (1998).
Output of Node 8
Development and auto-calibration of an input- 1.3
output model of chlorine transport in drinking 1.2
water distribution systems. IEEE Transactions 1.1

on Control Systems Techonology, 6(4), 543-553. 1

0.9

0.8
0.14
0.7
0 50 100 150 200 250 300

0.12 Chlorine Injection at Node 1

2

1.9
0.1

1.8

0.08 1.7

1.6
0.06
1.5
0 50 100 150 200 250 300

0.04 Fig.5 Constraints Fulfilment After Safety Zone

Introduction
0.02 (Dashed Line: Upper Limit
Dash-Dot Line: Lower Limit)
0
0 50 100 150 200 250 300

Fig. 2 Piece-wise Constant Linear Model: Parameter

Example