NONLINEAR PREDICTIVE CONTROL OF A BENCHMARK CSTR

H. Chen, A. Kremling and F. Allgöwer∗

Institut für Systemdynamik und Regelungstechnik, Universität Stuttgart
Pfaffenwaldring 9, 70550 Stuttgart, Germany

Keywords: Benchmark problem, reactor control, non- • A complete set of performance objectives including
linear control, nonlinear model predictive control an uncertainty description is given.
A more detailed discussion on the reasons and implica-
Abstract tions of the first two points can be found in [5].
Furthermore, we propose a solution to the control of
The intention of this paper is twofold: First, we want this benchmark problem based on on-line optimization
to propose a benchmark problem for nonlinear control in a receding horizon manner. Such control schemes are
system design. This problem involves control of a highly usually termed nonlinear model predictive controllers [9].
nonlinear chemical reactor that exhibits interesting prop- The paper is structured as follows: In Section 2 a com-
erties, like a change of steady state gain at the main op- plete description of the reactor is given including the
erating point. Secondly, we give a “reference” solution mathematical model with parameters. Section 3 briefly
to this benchmark problem based on a nonlinear model describes the specific operating point and the control
predictive control scheme. Despite the difficulty of the problem. We suggest a multi-input and a single-input
problem this controller achieves stability and good per- control problem. Finally, we propose a “reference” so-
formance to setpoint changes and disturbances in a ro- lution to those control problems in Section 4 where we
bust way. explain the structure and parameters of the controller
used.
1 Introduction
2 Description of the CSTR [4]
For economic and chemical engineering reasons, it is of-
ten desirable to operate chemical reactors such that the The reactor under consideration is a continuous stirred
production of a wanted product is maximized. For this, a tank reactor with a cooling jacket in which cyclopentenol
reactor is considered at an operating point where optimal is produced from cyclopentadiene by acid-catalyzed elec-
yield with respect to a desired product is achieved. How- trophylic hydration in aqueous solution. This reactor was
ever, operation at this point can considerably complicate first described in [4] and is adopted unchanged in this
the design of control system, and is then the motivation paper. The reaction scheme and parameters are derived
to propose a benchmark problem for nonlinear control by theoretical modeling based on physical properties de-
system design, which is based on a specific continuous scribed in the literature for a real process. Details on the
stirred tank reactor (CSTR) that is described in [1, 4]. derivation of the chemical parameters and the physico-
The benchmark problem that will be described in Sec- chemical background can also be found in [4].
tions 2 and 3 is characterized by a number of interesting Fig. 1 shows a schematic diagram of the reactor. The
features: main reaction is given by the transformation of cyclopen-
tadiene (substance A) to the product cyclopentenol (sub-
• The steady state gain changes its sign at the oper-
stance B). The initial reactant cyclopentadiene also re-
ating point. Thus, linear controllers (with integral
acts in an unwanted parallel reaction to the by-product
action) will not be able to stabilize this reactor and
dicyclopentadiene (substance D). Furthermore, cyclopen-
accomplish satisfying performance [6].
tanediol (substance C) is formed in an unwanted con-
• The zero dynamics changes its stability property at secutive reaction from the product cyclopentenol. This,
this operating point. Therefore, the qualitative be- so-called van der Vusse reaction, is described by the fol-
havior of the CSTR differs considerably for different lowing reaction scheme:
setpoints and disturbances. k1 k2
A −→ B −→ C
(1)
• The problem has a “real world” background. 3
2A −→
k
D .
∗ author to whom all correspondence should be addressed:

phone: +49-711-6856193, Fax: +49-711-6856371, The flow V̇ fed to the reactor contains only cyclopenta-
email: allgower@rus.uni-stuttgart.de diene (substance A) with concentration cA0 and temper-

Proc. 3rd European Control Conference ECC’95, Rome/Italy 3247

 erance of 0.02 mol
l ). The yield Φ of product B is defined
V, A
as the ratio between product concentration cB and con-
centration of initial reactant cA0 in the feed1
cB |S
Φ= , (3)
cA0
mK and is a measure for the effectiveness of the production.
VR
This optimal operating point is found by optimization
of the steady state yield with respect to the design vari-
ables steady state feed flow VV̇R |S , steady state heat re-
moval Q̇K |S and feed temperature ϑ0 |S . It is described
by the following values:
QK
V, A, B, C, D
mol mol
cA0 |S = 5.10 cA |S = 2.14
l l
Figure 1: Schematic representation of the CSTR mol
ϑ0 |S = 104.9 oC cB |S = 1.09
l
 (4)
ature ϑ0 . The heat removal Q̇K can be withdrawn from V̇ 
the coolant by an external heat exchanger.  = 14.19 h−1 ϑ|S = 114.2 C o
VR 
S
The dynamics of the reactor can be described by the  kJ
following nonlinear differential equations that are derived Q̇K  = −1113.5 ϑK |S = 112.9 oC .

from component balances for substances A and B and S h
from energy balances for the reactor and cooling jacket: The operating point considered here is different from
the one suggested in [4]: The yield is increased by
V̇
ċA = (cA0 − cA ) − k1 (ϑ)cA − k3 (ϑ)cA 2 (2a) about 20%; the selectivity2 with respect to initial reac-
VR tant A is also increased by more than 50%; and in ad-
V̇ dition, the energy consumption is considerably reduced.
ċB = − cB + k1 (ϑ) cA − k2 (ϑ)cB (2b)
VR Thus, this operating point of optimal yield is very desir-
V̇ 1  able for economic reasons.
ϑ̇ = (ϑ0 − ϑ) − k1 (ϑ)cA ∆HRAB By use of control we want to guarantee the production
VR ρCp
 of cyclopentenol with desired purity cB |s despite varia-
+k2 (ϑ)cB ∆HRBC + k3 (ϑ)cA 2 ∆HRAD tions in the feed temperature ϑ0 . We distinguish two
kw AR control problems: A multi-input problem and a (more
+ (ϑK − ϑ) (2c) demanding) single-input problem.
ρCp VR
1  
ϑ̇K = Q̇K + kw AR (ϑ − ϑK ) , (2d) 3.1 Multi-input control problem
mK CP K
cA ≥ 0, cB ≥ 0 . For the multi-input control problem we suggest that the
flow rate normalized by the reactor volume VV̇R and the
The concentrations of substances A and B are cA and cB
respectively. The temperature in the reactor is denoted heat removal Q̇K are used as manipulated variables. We
by ϑ, while the temperature in the cooling jacket is given assume that the product concentration cB and the reac-
by ϑK . The reaction velocities ki are assumed to depend tor temperature ϑ can be measured. However, only the
on the temperature via the Arrhenius law concentration cB is of interest to be controlled. We want
  to be able to produce substance B with concentration cB
Ei
ki (ϑ) = ki0 ·exp , i = 1, 2, 3 . (2e) in the following range (setpoint):
ϑ/o C + 273.15
mol mol
Values for the physical and chemical parameters in 0.8 ≤ cB |s ≤ 1.09 .
l l
equations (2a)–(2e) are given in Table 1. It should
be noted that most parameters are only known within The feed of the reactor is assumed to come from an
bounds. upstream unit. Therefore, the feed temperature ϑ0 is
assumed to vary between
3 Control problem at the point of opti- 100 oC ≤ ϑ0 ≤ 115 oC
mal yield
and is considered as unmeasurable disturbance.
The reactor is operated at a point where optimal yield 1 By |S we denote the steady state value of a variable.
with respect to the product B is achieved (within a tol- 2 The selectivity S is defined as S = c cB
−c
.
A0 A

3248
 Table 1: Physico-chemical parameters for the CSTR

Name of parameter Symbol Value of parameter

collision factor for reaction k1 k10 (1.287 ± 0.04) · 1012 h−1

collision factor for reaction k2 k20 (1.287 ± 0.04) · 1012 h−1
1
collision factor for reaction k3 k30 (9.043 ± 0.27) · 109 mol A·h
activation energy for reaction k1 E1 −9758.3K
activation energy for reaction k2 E2 −9758.3K
activation energy for reaction k3 E3 −8560K
kJ
enthalpies of reaction k1 ∆HRAB (4.2 ± 2.36) mol A
kJ
enthalpies of reaction k2 ∆HRBC −(11.0 ± 1.92) mol B
kJ
enthalpies of reaction k3 ∆HRAD −(41.85 ± 1.41) mol A
density ρ (0.9342 ± 4.0 · 10−4 ) kg
l
kJ
heat capacity Cp (3.01 ± 0.04) kg·K
kJ
heat transfer coefficient for cooling jacket kw (4032 ± 120) h·m 2 ·K

surface of cooling jacket AR 0.215 m2

reactor volume VR 0.01 m3
coolant mass mK 5.0 kg
kJ
heat capacity of coolant CP K (2.0 ± 0.05) kg·K

considering the two parameter sets given in Table 2. Pa-

Table 2: Two extreme cases for parameter uncertainty
rameters not listed in Table 2 are assumed to have their
nominal value. Those sets are related to two extreme
Parameter Case 1 Case 2 cases chosen by physical considerations, that represent a
kind of worst-case deviation from the nominal values.
k10 [h−1 ] 1.327 · 1012 1.247 · 1012
k20 [h−1 ] 1.327 · 1012 1.247 · 1012 3.2 Single-input control problem
k30 [ mol1A·h ] 8.773 · 109 9.313 · 109 The second control problem of interest involves derivation
kJ
∆HRAB [ mol of a single-input controller. Here, only the flow rate VV̇R
B] 6.56 1.84
kJ is available as manipulated variable. Heat removal Q̇K
∆HRBC [ mol A] −9.08 −12.92
is held constant at the steady state value given by (4).
kJ
∆HRAD [ mol A ] −40.44 −43.26 We consider the same disturbances and constraints in the
manipulated variable VV̇R as in the multi-input case. Only
the range of setpoint values is reduced to
The controller has to compensate the effects of changes
in the setpoint value cB |s and of disturbance in ϑ0 simul- mol mol
0.95 ≤ cB |s ≤ 1.09 .
taneously. The manipulated variables are constrained by l l
Again the same test signals and uncertainty sets can be
V̇
3 h−1 ≤ ≤ 35 h−1 chosen for testing the controllers.
VR
kJ kJ
−9000 ≤ Q̇K ≤ 0 . 4 Nonlinear predictive control of the
h h
In order to test the performance of the controller, we CSTR
suggest step changes in the setpoint and feed tempera-
4.1 Model predictive control
ture from their maximal value to the minimal value and
back. The maximal steady state offset should not exceed Model predictive control (MPC) has become an attrac-
0.02 mol
l (control tolerance). Robustness can be tested by tive feedback strategy, especially for linear or nonlin-

3249
ear plants subject to input and state constraints. Dur- lead to poor performance or even instability. Thus, a
ing the last decade, many formulations have been devel- state estimator is needed in the case of output feedback.
oped for linear and nonlinear, stable and unstable, and Here, we use a continuous-discrete extended Kalman fil-
non-minimum-phase plants (e.g. [2, 8, 9]). This method ter (EKF) [3], which is an extension of the linear Kalman
has been successfully applied also in the process indus- filter to the more general case of nonlinear system with
try (e.g. [10]). In MPC the controller predicts the behav- discrete output measurements
ior of a plant over a prediction horizon using the plant
model and measurement, and determines a manipulated ẋ(t) = f (x(t)) + w(t) (6a)
variable sequence that optimizes some open-loop perfor- yk = g (x(tk )) + v k , (6b)
mance objective over the prediction horizon. This ma-
nipulated variable sequence is implemented until the next where w(t) and {v k } are random variables that represent
measurement becomes available. Then, the optimization noise. For details see for example [3].
problem is solved again. For a more detailed description of NMPC we refer for
The general formulation for nonlinear model predictive example to [2, 7, 9].
control (NMPC) may be stated (by some abuse of nota-
tion) as 4.2 Simulation results

min J uNc , x(t), y(t)

(5a) In this section we present a “reference” solution to the
u Nc benchmark problem described in Sections 2 and 3. The
objective function for NMPC is taken as the integral
subject to
square error of the controlled variable from its setpoint:
= f x(t), uNc ,


ẋ(t) t ∈ (tk , tk+Np ] (5b) Z tk+Np

2
J = cB (t) − cB |s (t) dt . (7)
= g x(t), uNc


y(t) (5c) t k

and It is very time consuming to solve the on-line optimiza-

tion problem posed. The complexity depends mainly on
Nc


h1 x(t), u = 0 (5d) the number of independent variables in uNc . In order to
Nc reduce the number of independent variables, a technique


h2 x(t), u ≥ 0 (5e)
called “blocking” is used [7]. The idea is not to allow the
with initial condition manipulated variables to vary at every future sampling
time but to require them to be constant over several sam-
x(tk ) = xk , (5f) pling periods. In addition, we also set the control hori-
zon Nc smaller than the prediction horizon Np , and let
where x is the state vector, y is the vector of controlled
the manipulated variables be constant after the control
variables, Np and Nc are the length of the prediction
horizon. In physical terms, a manipulated variable se-
and control horizon respectively corresponding to sam-
quence of only Nc steps cannot make the system follow
pling period ∆T , satisfying Np ≥ Nc . The scalar function
the setpoint exactly over all Np steps, when Nc < Np [7].
J is the performance objective, equations (5b) and (5c)
Therefore, only the setpoint change in the control horizon
represent the plant dynamics and output function, equa-
is considered for each optimization.
tions (5d) and (5e) describe equality and inequality con-
With the EKF not only the four states (cA , cB , ϑ, ϑK )
straints which have to be satisfied in practical application
but also the unmeasurable disturbance (feed tempera-
(for example, actuator saturation). The sequence
ture ϑ0 ) are estimated. In order to estimate the distur-
uNc = [uk , uk+1 , · · · , uk+Nc −1 ]T (5g) bance it is assumed that the disturbance is constant but
unknown, thus, satisfying the differential equation
is the discrete manipulated variable vector, with uk+i =
ϑ̇0 (t) = 0 . (8)
u(tk + i∆T ), i = 0, 1, · · · , Nc − 1. The manipulated vari-
ables are considered in a discretized way, because it is in Here, the reactor temperature ϑ is used as secondary mea-
general unrealistic to solve optimization problem (5) an- surement to enhance the estimation performance.
alytically over the prediction horizon Np with continuous
manipulated variables.
4.2.1 Controller tuning
In this formulation it is assumed that the full state
measurement is available at time tk , i.e., xk in initial con- The NMPC controller with EKF has a number of tun-
dition (5f) is known. However, in most applications, not ing parameters (among others the prediction horizon Np ,
all of the states can be measured directly. In this case the the control horizon Nc , and the parameters for the EKF).
initial condition can merely be set equal to the states pre- There are no general criteria to achieve stability and ro-
dicted at the previous sampling time. Of course, model- bustness for nonlinear model predictive control. In par-
plant mismatch and unmeasurable disturbances will then ticular, for NMPC with EKF no separation principle is

3250
known to hold. Hence, we choose the parameters for It is shown in Fig. 3 that the EKF quickly recovers from
the EKF in such a way that it has a good estimation the initial estimation error and gives very good estimates
performance. It is known that for linear non-minimum-
4
phase systems the closed-loop system can become unsta-
ble if the control action is too aggressive (i.e. the predic- 3.5

Concentration
Multi-input
tion horizon is too “short” or the control horizon is too

cA [ mol ]
3



ĉA


l
“long”) [7]. This happens also in the nonlinear case, espe- 2.5 Single-input
cially when the operating point is at the point of optimal 2 H
Y
0 HcA 500 1000 1500 2000 2500
yield as in the case of the CSTR considered. Here, we Time [s]

choose the prediction horizon to be Np = 200, the con- 120

Coolant Tempera-
 K
ϑ̂
trol horizon to be Nc = 3 with a sampling period of 20 115 

seconds. And we also make the manipulated variables to 110
I
@
@ϑ 

ture [o C]
Single-input
maintain constant over two sampling periods.
 Multi-input
K
105

100
0 500 1000 1500 2000 2500
4.2.2 Setpoint tracking Time [s]
110
H
Y
H ϑ̂

Temperature
Fig. 2 shows the closed-loop response of the CSTR to step 108 0 Multi-input, Single-input

changes in the setpoint cB |s from maximum to minimum ϑ0

and back to maximum value without disturbance. During

106 



[o C]
104

Feed
1.2
cB , cB |s [ mol ]

102
Single-input 0 500 1000 1500 2000 2500
l

H
Concentration

Time [s]
1 j
H
 Multi-input
Figure 3: State and disturbance estimation for setpoint
0.8 tracking ( , − · − true and estimated variable in the
0 500 1000 1500 2000 2500
single-input case; −−, · · · true and estimated variable in
40
Time [s] the multi-input case).
variable V̇ [ 1 ]
VR h

30
of the states and the disturbance for both the multi-input
Manipulated

20
and single-input case.
10
In the Sections 4.2.3 and 4.2.4 we assume that there
0
0 500 1000 1500 2000 2500
is no estimation error in the initial states, but the initial
Time [s]
0
estimated disturbance ϑ̂0 takes the nominal value in (4).
variable Q̇K [ kJ ]
h

-2000
Manipulated

-4000 4.2.3 Disturbance rejection

-6000

-8000 Because the operating point is at optimal yield, it is very

-10000
0 500 1000 1500 2000 2500
difficult to maintain this operating point especially when
Time [s]
the feed temperature ϑ0 is decreased to 100o C. Fig. 4
Figure 2: Setpoint tracking of the CSTR controlled with shows the closed-loop response of the CSTR to step dis-
NMPC ( single-input; −− multi-input; − · − setpoint turbances in the feed temperature. At time t = 0s
for the single-input case; · · · setpoint for the multi-input the feed temperature changes abruptly from the “steady
case). state” value ϑ0 = 104.9oC to 115oC and then at time t =
1500s down to 100oC, while the setpoint remains con-
the first 400 seconds the controller deals with the initial stant at the optimal value. The system is stable and only
estimation error. The true and estimated initial states a steady state offset of −0.008 moll (single-input), which
and disturbance (xk , ϑ0k and x̂k , ϑ̂0k ) are is smaller than the required control tolerance, can be ob-
  served. In principle, the NMPC controller is designed
xk
= [2.14 1.09 114.2 112.9 104.9]T , for zero offset, but the maximally achievable concentra-
ϑ0k
  tion cB in steady state is of course lower than the optimal
x̂k value of cB = 1.09 mol
= [2.50 1.09 114.2 114.0 110.0]T . l when the feed temperature is re-
ϑ̂0k
duced to ϑ0 = 100 C, ( VV̇R and Q̇K being free within their
o

Thus, a rather large initial error was assumed. constraints and the other parameters being not changed).
It can be seen from Fig. 2 that both in the multi-input It is clear that the manipulated variable sequence is not
and single-input case only three steps are needed to bring totally stationary at the last optimization step, because
the concentration of product B into the required control the controller always tries to overcome the offset. Alto-
tolerance (±0.02 mol
l ) of the setpoint without large over- gether, the controller achieves a very good disturbance
shoot. attenuation.

3251
 1.12 1.2
cB , cB |s [ mol ]

cB , cB |s [ mol ]
Setpoint Single-input


 
l

l


Concentration

Concentration
1.1
1
1.08

1.06 0.8
 Multi-input

1.04 0 500 1000 1500 2000 2500 3000 3500 4000

0 500 1000 1500 2000 2500 3000 3500 4000
Time [s]
Time [s]
25 40

variable V̇ [ 1 ]
variable V̇ [ 1 ]

VR h
VR h

20 30

Manipulated
Manipulated

15 20

10 10

5 0
0 500 1000 1500 2000 2500 3000 3500 4000 0 500 1000 1500 2000 2500 3000 3500 4000
Time [s] Time [s]
0 0

variable Q̇K [ kJ ]
variable Q̇K [ kJ ]

h
h

-2000
-2000

Manipulated
Manipulated

-4000
-4000
-6000
-6000 -8000

-8000 -10000
0 500 1000 1500 2000 2500 3000 3500 4000 0 500 1000 1500 2000 2500 3000 3500 4000
Time [s]
Time [s]

Figure 4: Disturbance attenuation of the CSTR con- Figure 5: Robustness to parameter uncertainties
trolled with NMPC ( single-input; −− multi-input). ( single-input; −− multi-input).

4.2.4 model predictive control scheme. Despite the difficulty of

Robustness to parameter uncertainties
the problem very satisfying control performance in both
The physico-chemical parameters of the benchmark prob- setpoint tracking and disturbance attenuation is robustly
lem are only known within bounds. To study the ro- achieved with the controller proposed.
bustness of the NMPC controller with EKF, two extreme
cases of model-plant mismatch as described in Section 3 References
will be discussed.
Fig. 5 shows how the NMPC controller with EKF [1] S. Engell and K.-U. Klatt. Nonlinear control of a non-
tries to make the system track setpoint changes at time minimum-phase CSTR. In Amer. Contr. Conf., pages
t = 400s and t = 2300s, against model-plant mismatch 2941–2945, Los Angeles, CA, 1993.
and disturbances: There is model-plant mismatch in the [2] C.E. Garcia, D.M. Prett, and M. Morari. Model Predic-
first 1800 seconds according to Case 1 and after that ac- tive Control: Theory and practice – A survey. Automat-
cording to Case 2 of Table 2; At time t = 1100s and ica, 25:335–347, 1989.
o
t = 3000s the feed temperature ϑ0 changes from 104.9 C [3] A. Gelb. Applied Optimal Estimation. M.I.T. Press,
down to 100.0oC and from 100.0oC up to 115.0oC respec- Cambridge, MA., 1974.
tively. Because of the uncertainties a stationary error in [4] K.-U. Klatt and S. Engell. Rührkesselreaktor mit
the state and disturbance estimation is observed and the Parallel- und Folgereaktion. In VDI-Bericht Nr. 1026,
setpoints cannot be held exactly. But the controlled vari- pages 101–108. VDI-Verlag, Düsseldorf, 1993.
able remains within the required tolerance of ±0.02 mol l
[5] K.-U. Klatt, S. Engell, A. Kremling, and F. Allgöwer.
for both sets of uncertainties and for both the multi-input Testbeispiel: Rührkesselreaktor mit Parallel- und Fol-
and single-input case. gereaktion. In S. Engell, editor, Entwurf Nichtlinearer
Regelungen. Oldenbourg Verlag, 1995.
[6] M. Morari. Robust stability of systems with integral con-
5 Conclusions trol. In Proc. 22nd IEEE Conf. Decision Contr., pages
865–869, Los Angeles, CA, 1983.
In this paper we proposed a benchmark problem for non-
[7] M. Morari, C.E. Garcia, J.H. Lee, and D.M. Prett. Model
linear control system design. The system to be controlled Predictive Control. unpublished, 1994.
is a realistic chemical reactor, that is operated at the
[8] K.R. Muske and J.B. Rawlings. Linear model predictive
point of optimal yield. This operating point is especially
control of unstable processes. J. Proc. Contr., 3(2):85–
desirable for economic reasons. Due to the operation at 96, 1993.
this highly nonlinear operating point, the control prob-
[9] A.A. Patwardhan, J.B. Rawlings, and T.F. Edgar. Non-
lem is very challenging. It can for example be shown that linear model predictive control. Chem. Eng. Commun.,
in the single-input case no linear controller with integral 87:123–141, 1990.
action can robustly stabilize the reactor.
[10] J. Richalet. Industrial applications of model based pre-
In the second part of the paper we gave a particular
dictive control. Automatica, 29(5):1251–1274, 1993.
solution to this benchmark problem based on a nonlinear

3252