12

Applications

This chapter is dedicated to presenting some MPC applications to the control of different real and simulated processes. The rst application presented
corresponds to a self-tuning and a gain scheduling GPC for a distributed collector eld of a solar power plant. In order to illustrate how easily the control
scheme shown in Chapter 5 can be used in any commercial control system,
some applications concerning the control of typical variables such as ows,
temperatures and levels of different processes of a pilot plant are presented.
The description of two applications in the food industry (a sugar renery and
an olive oil mill) are included. Finally the application of an MPC to a highly
nonlinear process (a mobile robot) is also described.

12.1 Solar Power Plant

This section presents an application of an adaptive long-range predictive
controller to a solar power plant and shows how this type of controller, which
normally requires a substantial amount of computation (in the adaptive
case), can easily be implemented with few computation requirements. The
results obtained when applying the controller to the plant are also shown.
The controlled process is the distributed collector eld (Acurex) of the
Solar Platform of Tabernas (Spain). The distributed collector eld consists
mainly of a pipeline through which oil is owing and onto which the solar
radiation is concentrated by means of parabolic mirrors, which follow the
sun by rotating on one axis, in order to heat the oil. It consists of 480 modules
arranged in 20 lines which form 10 parallel loops. A simplied diagram of
the solar collector eld is shown in Figure 12.1. The eld is also provided
with a sun-tracking mechanism which causes the mirrors to revolve around
an axis parallel to that of the pipeline.
On passing through the eld the oil is heated and then introduced into a
storage tank to be used for the generation of electrical energy. The hot oil can
also be used directly for feeding the heat exchanger of a desalination plant.
E. F. Camacho et al., Model Predictive Control
Springer-Verlag London Limited 2007

338

12 Applications
Storage tank
Acurex
field

.........

Pump

To
steam
generator
or
Desalination
Plant
.

Buffer

Fig. 12.1. Schematic diagram of collectors eld

The cold inlet oil to the eld is extracted from the bottom of the storage tank.
Each of the loops is formed by four twelve-module collectors, suitably
connected in series. The loop is 172 metres long, the active part of the loop
(exposed to concentrated radiation) measures 142 metres and the passive
part 30 metres.
The system is provided with a three-way valve which allows the oil to
be recycled in the eld until its outlet temperature is adequate to enter the
storage tank. A more detailed description of the eld can be found in [101].
A fundamental feature of a solar power plant is that the primary energy
source, whilst variable, cannot be manipulated. The intensity of the solar radiation from the sun, in addition to its seasonal and daily cyclic variations,
is also dependent on atmospheric conditions such as cloud cover, humidity,
and air transparency. It is important to be able to maintain a constant outlet
temperature for the uid as the solar conditions change, and the only means
available for achieving this is via adjustment of the uid ow.
The objective of the control system is to maintain the outlet oil temperature at a desired level in spite of disturbances such as changes in the solar
irradiance level (caused by clouds), mirror reectivity or inlet oil temperature. This is accomplished by varying the ow of the uid through the eld.
The eld exhibits a variable delay time that depends on the control variable
(ow). The transfer function of the process varies with factors such as irradiance level, mirror reectance and oil inlet temperature.
The distributed collector eld is a nonlinear system which can be approximated by a linear system when considering small disturbances. The
maintenance of a constant outlet temperature throughout the day as the so-

12.1 Solar Power Plant

339

lar conditions change requires a wide variation in the operational ow level.

This leads to substantial variations in the general dynamic performance and
in particular, from the control viewpoint, gives rise to a system time delay
which varies signicantly. The controller parameters need to be adjusted to
suit the operating conditions, and self-tuning control offers one approach
which can accommodate such a requirement.
Because of the changing dynamics and strong perturbations, this plant
has been used to test different types of controllers [40],[41].
For self-tuning control purposes a simple, linear model is required which
relates changes in uid ow to changes in outlet temperature.
Observations of step responses obtained from the plant indicate that in
the continuous time domain behaviour it can be approximated by a rstorder transfer function with a time delay. Since the time delay d is relatively
small compared to the fundamental time constant , a suitable discrete model
can be constructed by choosing the sample period T equal to the lowest value
of the time delay d . This corresponds to the operating condition where the
ow level is highest. The discrete transfer function model then takes the form
g(z 1 ) = z k

bz 1
1 az 1

and at the high ow level condition, k = 1.

12.1.1 Selftuning GPC Control Strategy
A particular feature of the system is the need to include a series feedforward
term in the control loop [39]. The plant model upon which the self-tuning
controller is based relates changes in outlet temperature to changes in uid
ow only. The outlet temperature of the plant, however, is also inuenced by
changes in system variables such as solar radiation and uid inlet temperature. During estimation, if either of these variables changes it introduces a
change in the system output unrelated to uid ow which is the control input
signal, and in terms of the model, it would result in unnecessary adjustments
of the estimated system parameters.
Since both solar radiation and inlet temperature can be measured, this
problem can be eased by introducing a feedforward term in series to the system, calculated from steady-state relationships, which makes an adjustment
in the uid ow input, aimed at eliminating the change in outlet temperature caused by the variations in solar radiation and inlet temperature. If the
feedforward term perfectly countered the changes in solar radiation and inlet temperature, then the observed outlet temperature changes would only
be caused by changes in the control input signal. Although exact elimination
obviously cannot be achieved, a feedforward term based on steady-state considerations overcomes the major problems inherent in the single-input model
and permits successful estimation of model parameters. The basic idea is to

340

12 Applications
Solar radiation
Inlet oil temperature

r(t)

l r1

+
+ +

1 - a^
b^

+
+

Feed
Plant

y(t)

forward

-1
Z
Ident
b^

^a
Adapt

ly2

Predictor

ly1

Fig. 12.2. Self-tuning control scheme

compute the necessary oil ow to maintain the desired outlet oil temperature
given the inlet oil temperature and the solar radiation. The feedforward signal also provides control benets when disturbances in solar radiation and
uid inlet temperature occur, but the basic reason for its inclusion is that of
preserving the validity of the assumed system model in the self-tuning algorithm; for more details see [45].
In the control scheme, the feedforward term is considered as a part of
the plant, using the setpoint temperature to the feedforward controller as the
control signal.
In this section, the precalculated method described in Chapter 5 is used.
As the dead time d is equal to 1, the control law is given by
u(t) = ly1 y(t + 1 | t) + ly2 y(t) + lr1 r(t)

(12.1)

where the value y(t + 1 | t) is obtained by use of the predictor:

y(t + j + 1 | t) = (1 + a)
y (t + j | t) a
y (t + j 1 | t) + b  u(t + j 1) (12.2)
The proposed control scheme is shown in Figure 12.2. The plant estimated parameters are used to compute the controller coefcients (ly1 , ly2 , lr1 )
via the adaptation mechanism. Notice that in this scheme, the feedforward
term is considered as a part of the plant (the control signal is the setpoint
temperature for the feedforward controller instead of the oil ow). This signal is saturated and ltered before its use in the estimation algorithm. The
controller also has a saturation to limit the increase of the error signal.
As suggested in Chapter 5, a set of GPC parameters were obtained for
(i) = 1, (i) = 5 and N = 15. The pole of the system has been changed with
a 0.0005 step from 0.85 to 0.95, which are the values that guarantee the system
stability if the parameter set estimation is not accurate enough. Notice that

12.1 Solar Power Plant

-2

2.5

-2.5

0.415
0.41
lr1

ly2

ly1

-3

1.5
0.85

-3.5
0.85

0.9
estimated pole

0.95

341

.405
0.4

0.9
estimated pole

0.95

0.395
0.85

0.9
estimated pole

0.95

Fig. 12.3. Controller parameters for = 5

due to the fact that the closed-loop static gain must equal the value unity, the
sum of the three parameters equals zero. The curves shown in Figure 12.3
correspond to the controller parameters ly1 ,ly2 , lr1 for the values of the pole.
The adjustment of analytical functions to the calculated values provide:
/ (1.11 - a
)
ly1 = 0.4338 - 0.6041 a
/ (1.082 - a
)
ly2 = -0.4063 + 0.4386 a
lr1 = ly1 ly2

(12.3)

These expressions give a very good approximation to the true controller

parameters and t the set of computed data with a maximum error of less
than 0.6 per % of the nominal values for the range of interest of the open
loop-pole; for more details see [40].
The parameters of the system model in the control scheme are determined online via recursive least-squares estimation. The estimation algorithm incorporates a variable forgetting factor and only works when the input signal contains dynamic information. These considerations can be taken
into account by checking the following conditions:
| u | A
k=0


| u(k) | B

k=N

If one of these conditions is true, the identier is activated. Otherwise, the

last set of estimated parameters is used. Typical values of A, B and N chosen
from simulation studies are: A = 9, 7 B 9 and N = 5. The covariance matrix P (k) is also monitored by checking that its diagonal elements
are kept within limits; otherwise P (k) is reset to a diagonal matrix having
the corresponding elements saturated to the violated limit.
With respect to the adaptation mechanism, it only works when the estimated parameters are contained within the ranges (0.85 a
0.95 and
))
0.9 kest 1.2, where kest is the estimated static gain of the plant b/(1 a
in order to avoid instability in cases of nonconvergence of the estimator. A
backup controller is used in situations in which these conditions are not accomplished (for example, when daily operation starts).

342

12 Applications

In each sampling period the self-tuning regulator consists of the following steps:
1. Estimate the parameters of a linear model by measuring the inlet and
outlet values of the process.
2. Adjust the parameters of the controller using Expressions (12.3).
3. Compute y(t + d | t) using the predictor in (12.2).
4. Calculate the control signal using (12.1).
5. Supervise the correct working of the control.
Plant Results
Figure 12.4 shows the outlet oil temperature and reference when the proposed self-tuning generalised predictive controller is applied to the plant.
The value of the control weighting was made equal to 5 and, as can be
seen, a fast response to changes in the setpoint is obtained (the rising time observed is approximately 7 minutes). When smaller overshoots are required,
the control weighting factor has to be increased.
The evolution of the irradiation for the test can be seen in the same gure
and it corresponds to a day with small scattered clouds. The oil ow changed
from 4.5 l/s to 7 l/s, and the controller could maintain performance in spite
of the changes in the dynamics of the process caused by ow changes.
12.1.2 Gain Scheduling Generalized Predictive Control
There are many situations in which it is known how the dynamics of a
process change with the operating conditions. It is then possible to change
the controller parameters taking into account the current operating point of
the system. Gain scheduling is a control scheme with open-loop adaptation,
which can be seen as a feedback control system in which the feedback gains
are adjusted by a feedforward compensation. Gain scheduling control is a
nonlinear feedback of a special type: it possesses a linear controller whose
parameters are modied depending on the operating conditions in a prespecied manner.
The working principle of this kind of controller is simple; it is based on
the possibility of nding auxiliary variables which guarantee a good correlation with process changing dynamics. In this way, it is possible to reduce
the effects of variations in the plant dynamics by adequately modifying the
controller parameters as functions of auxiliary variables.
An essential problem is the determination of the auxiliary variables. In
the case studied here, the behaviour and changes in the system dynamics
mainly depend on the oil ow if very strong disturbances are not acting on
the system (due to the existence of the feedforward controller in series with
the plant). The oil ow has been the variable used to select the controller

12.1 Solar Power Plant

343

set point/outlet temperatures ( C)

190.0

180.0

170.0

160.0
11.0

11.5

12.0

12.5
13.0
local time (hours)

13.5

14.0

14.5

7.0

oil flow (l/s)

6.5

6.0

5.5

5.0

4.5
11.0

11.5

12.0

12.5
13.0
local time (hours)

13.5

14.0

14.5

730.0

Solar radiation (W/m2 )

710.0

690.0

670.0

650.0

630.0
11.0

11.5

12.0

12.5
13.0
local time (hours)

13.5

14.0

14.5

Fig. 12.4. Adaptive GPC: plant outlet oil temperature, ow and solar radiation

parameters (a low-pass lter is used to avoid the inclusion of additional dynamics due to sudden variations in the controller parameters).
Once the auxiliary variables have been determined, the controller parameters have to be calculated at a number of operating points, using an
adequate controller design algorithm, which in this case is the GPC methodology. When coping with gain scheduling control schema, stability and performance of the controlled system are usually evaluated by simulation studies [142]. A crucial point here is the transition between different operating

344

12 Applications

points. In those cases in which a nonsatisfactory behaviour is obtained, the

number of inputs to the table of controller parameters must be augmented.
As has been mentioned, it is important to point out that no feedback exists
from the behaviour of the controlled system to the controller parameters. So
this control scheme is not considered as an adaptive one, but rather as a special case of a nonlinear controller.
The main disadvantages of gain scheduling controllers are:
It is an open-loop compenzation: there is no way to compensate for a
wrong election of the controller parameters within the table.
Another important inconvenience is that the design stage of the strategy often consumes too much time and effort. The controller parameters
must be calculated for enough operating points, and the behaviour of the
controlled system has to be checked under very different operating conditions.
Its main advantage is the ease of changing controller parameters in spite
of changes in process dynamics. As classical examples of applications of this
kind of controller, the following control elds can be mentioned: design of
ship steering autopilots, pH control, combustion control, engine control, design of ight autopilots, etc. [12].
Plant Models and Fixed Parameter Controllers
The frequency response of the plant has been obtained by performing a PRBS
PRBS test in different operating conditions, using both the plant and a nonlinear distributed parameter model1 [25]. In this way, different linear models
were obtained from input-output data in different working conditions. These
models relate changes in the oil ow to those of the outlet oil temperature,
and can take into account the antiresonance characteristics of the plant if they
are adequately adjusted. The control structure proposed is shown in Figure
12.5. As can be seen, the output of the generalized predictive controller is the
input (trf f ) of the series compensation controller, which also uses the solar
radiation, inlet oil temperature and reectivity to compute the value of the
oil ow, which is sent to the pump controller.
The controller parameters were obtained from a linear model of the plant.
From input-output data of the plant, the degrees of the polynomials A and B
and the delay (of a CARIMA plant model) that minimizes Akaikes Information Theoretic Criterion (AIC) were found to be na = 2, nb = 8 and d = 0. By a
least squares estimation algorithm, the following polynomials were obtained
using input-output data of one test with oil ow of around 6 l/s:
1

The software simulation package for the solar distributed collector eld with real
data from the plant can be obtained by contacting the authors or by accessing
http://www.esi.us.es/ eduardo.

12.1 Solar Power Plant

Inlet Temperature (Tin)

set
point
sp

GPC
Controller

trff

l[0]...l[12]

Series
compensator

345

Radiation

oil
flow

Plant

tout

U[1]

-1

U[2]

..
.
.
.
.

-1

U[9]

-1

Y[1]

-1

Y[2]

Fig. 12.5. Control scheme using high-order models

A(z 1 ) = 1 1.5681z 1 + 0.5934z 2

B(z 1 ) = 0.0612 + 0.0018z 1 0.0171z 2 + 0.0046z 3 + 0.0005z 4
+0.0101z 5 0.0064z 6 0.015z 7 0.0156z 8
The most adequate value for the control horizon (N = 15) was calculated
taking into account the values of the fundamental time constant and the sampling period used for control purposes. In this case, N1 = 1 and N2 = 15. The
value of was determined by simulation studies using the nonlinear model
and was found to be = 6 (fast) and = 7 (without overshoot). For smaller
values of , faster and more oscillatory responses were obtained. Following
the design procedure of the GPC methodology, the controller parameters corresponding to = 7 were obtained (Table 12.1).
Table 12.1. Fixed GPC controller coefcients
l[0]
l[1]
l[2]
l[3]
l[4]
l[5]
l[6]
2.4483 6.8216 4.7091 0.0644 0.0526 0.0084 0.0629
l[7]
l[8]
l[9]
l[10]
l[11]
l[12]
0.0161 0.0311 0.0631 0.0231 1.0553 0.3358

346

12 Applications

260.0

820.0

780.0
direct solar radiation (W/m2)

set point and outlet oil temperatures (C)

800.0

250.0

240.0

760.0
740.0
720.0
700.0
680.0
660.0
640.0
620.0

230.0
11.6

12.1

12.6

13.1
13.6
local time (hours)

14.1

14.6

set point and outlet oil temperatures o C

600.0
11.6

12.1

12.6

13.1
13.6
local time (hours)

14.1

14.6

direct solar radiation Wm2 

Fig. 12.6. Test with the xed GPC high-order controller

The control law can be written by

trf f = l[2]tout + l[1]y[1] + l[0]y[2] + l[6]u[6] + l[3]u[9] + l[4]u[8] + l[5]u[7]
+l[7]u[5] + l[8]u[4] + l[9]u[3] + l[10]u[2] + l[11]u[1] + l[12]sp (12.4)
where
trf f : reference temperature for the feedforward controller;
tout : outlet temperature of the eld;
sp: setpoint temperature;
l[i]: controller parameters;
y[i]: outlet temperature of the eld at sampling time (t i)and
u[i]: reference temperature for the feedforward controller at sampling time
(t i).
With these values, the behaviour of this xed parameter controller was analyzed in operation with the distributed solar collector eld. The outlet oil
temperature of the eld evolution and corresponding setpoint can be seen
in Figure 12.6. The evolution of the solar radiation during this test can also
be seen in Figure 12.6. Although direct solar radiation goes from 810 W/m2
to 610 W/m2 , the eld was working in midow conditions because the setpoint was also changed from 258o C to 230o C. When operating conditions in
the eld change, the dynamics of the plant also change and the controller
should be redesigned to cope with the control objectives.
The dynamics of the eld are mainly dictated by the oil ow, which depends on the general operating conditions of the eld: solar radiation, reectivity, oil inlet temperature, ambient temperature and outlet oil temperature setpoint. These changes in plant dynamics are illustrated in Figure 12.7,
where the frequency response of the nonlinear distributed parameter dynamic model of the eld can be seen. The curves shown in Figure 12.7 were
obtained by a spectral analysis of the input-output signals of the model at
different operating points (PRBS signals were used for the input).

12.1 Solar Power Plant

347

AMPLITUDE PLOT

amplitude

10

-1

10

flow=2.8 l/s
flow=5.2 l/s
flow=7.9 l/s
flow=9.3 l/s

-2

10

-3

10

-2

10
frequency (rad/s)

-1

10

Fig. 12.7. Frequency response of the eld under different operating conditions

As can be seen, the frequency response changes signicantly for different

operating conditions. The steady-state gain changes for different operating
points, as does the location of the antiresonance modes.
Taking into account the frequency response of the plant and the different
linear models obtained from it, it is clear that a self-tuning controller based on
this type of model is very difcult to implement. The fundamental reason is
the fact that estimation of the model parameters requires a lot of computation
when the number of estimated parameters increases and the convergence of
the estimation process is seldom accurate and fast enough.
Another way of coping with changing dynamics is to use a gain scheduling controller, making the controller parameters dependent on some variables which indicate the operating conditions.
With the input-output data used to obtain the frequency responses shown
in Figure 12.7 and using the method and type of model previously described
for the case of a high-order xed parameter controller, process (a[i] and b[i])
and controller (l[i]) parameters were obtained for several oil ow conditions
(q1 = 2.8 l/s, q2 = 5.2 l/s, q3 = 7.9 l/s and q4 = 9.3 l/s), using different values of the weighting factor . Tables 12.2 and 12.3 contain model and control
parameters, respectively, for a weighting factor = 6. A value of = 7 has
also been used to obtain responses without overshoot.
The controller parameters applied in real operation are obtained by using a linear interpolation with the data given in Table 12.3. It is important
to point out that to avoid the injection of disturbances during the controller
gain adjustment, it is necessary to use a smoothing mechanism of the transition surfaces of the controller gains. In this case, a linear interpolation in
combination with a rst-order lter has been used, given a modied ow
Q(t) = .95 Q(t 1) + .05 q(t) (where q(t) is the value of the oil ow at
instant t and Q(t) is the ltered value used for controller parameter adjustment). The linear interpolation has also been successfully applied in [97]. Another kind of gain scheduling approach can be obtained by switching from
one set of controller parameters to another depending on the ow condi-

348

12 Applications
Table 12.2. Coefcients of polynomials A(z 1 ) and B(z 1 ) for different ows
a[1]
a[2]
b[0]
b[1]
b[2]
b[3]
b[4]
b[5]
b[6]
b[7]
b[8]

q1
1.7820
0.81090
0.00140
0.03990
0.0182
0.0083
0.00060
.00001
0.00130
0.00160
0.00450

q2
1.438
0.5526
0.0313
0.0660
.0272
0.0071
0.0118
0.0138
0.0098
0.0027
.0054

q3
1.414
0.5074
0.0687
0.0767
.0392
0.0127
0.0060
.0133
.0156
.0073
0.0037

q4
1.524
0.7270
0.0820
0.0719
.0474
0.0349
0.0098
.0031
0.0111
0.0171
0.0200

Table 12.3. GPC controller coefcients in several operating points ( = 6)

l[0]
l[1]
l[2]
l[3]
l[4]
l[5]
l[6]
l[7]
l[8]
l[9]
l[10]
l[11]
l[12]

7.0481
16.2223
9.5455
0.03910
0.00980
0.00560
0.0070
0.0016
0.0793
0.1575
0.36470
0.82620
0.37130

1.4224
3.84390
2.7794
0.0139
0.00830
0.02610
0.03390
0.02480
0.00880
0.0822
0.16410
0.83010
0.35800

1.1840
3.48440
2.6527
0.00860
0.0184
0.0352
0.0239
0.02630
0.03980
0.0869
0.19600
0.89360
0.35230

1.3603
3.02280
2.0142
0.03740
0.02730
0.01080
0.0197
0.00460
0.05070
0.1098
0.12480
0.87390
0.35170

tions, without interpolating between controller parameters. The set of controller parameters c can be obtained by choosing one of the sets ci in Table
12.3, related to ow conditions qi (i = 1, 2, 3, 4):
if

qi + qi+1
qi1 + qi
<q
then c = ci , i = 2, 3
2
2
if q q1 then c = c1
if q q4 then c = c4

The control structure is similar to the one obtained for the xed controller
previously studied. The optimal realization of the gain scheduling controller
consists of calculating the controller parameters under a number of operating conditions and suppose that the values of the controller coefcients are

12.1 Solar Power Plant

349

222.0

set point/outlet temperatures (C)

220.0

218.0

216.0

214.0

212.0

210.0

208.0
15.0

15.2

15.5
15.8
local time (hours)

16.0

15.5
15.8
local time (hours)

16.0

Setpoint and oil outlet temperatures o C

820.0
800.0

beam radiation (W/m2)

780.0
760.0
740.0
720.0
700.0
680.0
660.0
640.0
15.0

15.2

Direct solar radiation (W/m2 )

Fig. 12.8. Test with the gain scheduling GPC controller, = 7

constant between different operating conditions, generating a control surface

based on an optimization criterion which takes into account the tracking error and control effort. It is evident that if the procedure is applied at many
working points, an optimum controller will be achieved for those operating
conditions if there is a high correlation between the process dynamics and
the auxiliary variable. The drawback to this solution is that the design process becomes tedious. This is one of the main reasons for including a linear
interpolation between the controller parameters.
Plant Results
In the case of real tests, similar results were obtained and depending on the
operating point, disturbances due to passing clouds, inlet oil temperature
variations, etc., different performance was achieved.
Figure 12.8 shows the results of one of these tests with the gain scheduling GPC with = 7. The operating conditions correspond to a clear afternoon with the solar radiation changing from 800 W/m2 to 660 W/m2 and oil
ow changing from 3.75 l/s to 2 l/s. As can be seen, the effect of the antireso-

350

12 Applications

set point/outlet temperatures (C)

290.0

285.0

280.0

275.0

270.0
11.7

12.2

12.7

13.2

13.7
14.2
local time (hours)

14.7

15.2

15.7

Setpoint and outlet oil temperature ( C)

870.0

860.0

beam radiation (W/m2)

850.0

840.0

830.0

820.0

810.0

800.0
11.7

12.2

12.7

13.2
13.7
14.2
local time (hours)

14.7

15.2

15.7

Direct solar radiation (W/m )

210.0

inlet oil temperature (C)

200.0

190.0

180.0

170.0
11.7

12.2

12.7

13.2

13.7
14.2
local time (hours)

14.7

15.2

15.7

Inlet oil temperature ( C)

Fig. 12.9. Test with the gain scheduling GPC controller, = 7

nance modes does not appear in the response, due to the use of an extended
high order model which accounts for these system characteristics.
Figure 12.9 shows the result of a test with a weighting factor = 7 corresponding to a day of intermittent scattered clouds which produce large
changes in the solar radiation level and the inlet oil temperature changing

12.1 Solar Power Plant

351

240.0
730.0

690.0
direct solar radiation (W/m2)

set point and outlet oil temperatures (C)

710.0
235.0

230.0

225.0

220.0

670.0
650.0
630.0
610.0
590.0
570.0
550.0

215.0
11.8

12.3

12.8

13.3
13.8
local time (hours)

14.3

14.8

set point and outlet oil temperatures  C

o

530.0
11.8

12.3

12.8

13.3
13.8
local time (hours)

14.3

14.8

direct solar radiation Wm 

2

Fig. 12.10. Test with the gain scheduling GPC controller, = 7

from 170o C to 207o C. As can be seen, the outlet oil temperature follows the
setpoint in spite of changing operating conditions and the high level of noise
in the radiation level produced by clouds.
The results of a test corresponding to a day with sudden changes in the
solar radiation caused by clouds can be seen in Figure 12.10. As can be seen,
the controller (also designed with = 7) is able to handle different operating
conditions and the sudden perturbations caused by the clouds. After the tests
presented using a weighting factor = 7, two new test campaigns were
carried out to test the behaviour of the controller with a weighting factor
= 6. In the rst campaign, the evaluation of the controller performance was
considered. In the second, the behaviour of the controller operating under
extreme working conditions was studied.
Figure 12.11 shows the results obtained in the operation on a day with
normal levels of solar radiation, but on which a wide range of operating conditions is covered (oil ow changing between 2 l/s and 8.8 l/s) by performing several setpoint changes. At the start of operation there is an overshoot
of 6o C, due to the irregular conditions of the oil owing through the pipes
because the operation starts with a high temperature level at the bottom of
the storage tank. After the initial transient, it can be observed that the controlled system quickly responds to setpoint changes under the whole range
of operating conditions with a negligible overshoot. The rise time is about 6
minutes with a setpoint change of 15 degrees, as can be seen in Figure 12.11,
with smooth changes in the control signal, constituting one of the best controllers implemented at the plant. It is important to note that the controller
behaves well even with great setpoint changes.

352

12 Applications
920.0

245.0
235.0

900.0

temperature (C)

215.0
205.0
195.0
185.0

inlet oil temperature

outlet oil temperature
set point temperature

175.0

direct solar radiation (W/m2)

225.0

880.0

860.0

840.0

165.0
155.0
145.0
10.1 10.6 11.1 11.6 12.1 12.6 13.1 13.6 14.1 14.6 15.1 15.6 16.1
local time (hours)

820.0
10.1 10.6 11.1 11.6 12.1 12.6 13.1 13.6 14.1 14.6 15.1 15.6 16.1
local time (hours)

Fig. 12.11. Test with the gain scheduling GPC controller, = 6

12.2 Pilot Plant

In order to show how GPC can be implemented on an industrial SCADA, applications to the control of the most typical magnitudes found in the process
industry (ow, level, temperatures) are introduced in this section.
The tests are carried out on a pilot plant existing in the Departamento de
Ingeniera de Sistemas y Automatica of the University of Seville. The pilot plant
is provided with industrial instrumentation and is used as a testbed for new
control strategies which can be implemented on an industrial SCADA connected to it. This plant is basically a system using water as the working uid
in which various thermodynamic processes with interchange of mass and energy can take place. It essentially consists of a tank with internal heating with
a series of input-output pipes and recirculation circuit with a heat exchanger.
The design of the plant allows for various control strategies to be tested
in a large number of loops. Depending on the conguration chosen, it is possible to control the types of magnitudes most frequently found in the process
industry such as temperature, ow, pressure and level. For this, four actuators are available: three automatic valves and one electric heater that heats
the interior of the tank. Later some of the possible loops are chosen (considered as being independent) for implanting the GPC controllers.
12.2.1 Plant Description
A diagram of the plant which shows its main elements as well as the localization of the various instruments is given in Figure 12.12.
The main elements are:
feed circuit. The plant has two input pipes, a cold water one (at air temperature) and a hot water one (at about 70o C) with nominal ow and
temperature conditions of 10 l/min and 2 bar for the cold water and 5
l/min and 1 bar for the hot. The temperatures and ows of the inputs
are measured with thermocouples and orice plates, respectively, with
motorized valves for regulating the input ows.

12.2 Pilot Plant

Hot water

FT1

FT2

Cold water

TT1

V4

TT2

V5

353

Heat exchanger
TT4

PT1
LT1

V8

TT5

Tank
Heater

FT3

Pump
TT3

Waste

Fig. 12.12. Diagram of the pilot plant

tank. It has a height of 1 m and an interior diameter of 20 cm, it is thermically insulated, and it has an approximate volume of 31 l. It can work
pressurized (up to a limit of 4 bar) or at atmospheric pressure, depending
on the position of the vent valve. In its interior there is a 15 kW electric
resistance for heating, and an overow, an output pipe and another pipe
for recirculating the water through the exchanger.
recirculation circuit. The hot water in the tank can be cooled by entering
cold water through the cooling circuit. This circuit is composed of a centrifugal pump that circulates the hot water from the bottom of the tank
through a tube bundle heat exchanger returning at a lower temperature
at its top.
12.2.2 Plant Control
To control the installation there is an O RSI Integral Cube distributed control
system, composed of a controller and a supervisor connected by a local data
highway. The former is in charge of carrying out the digital control and analogous routines whilst the latter acts as a programming and communication
platform with the operator. On this distributed control system the GPC algorithms seen before will be implemented. This control system constitutes
a typical example of an industrial controller, having the most normal characteristics of medium-size systems to be found in the market today. As in
most control computers the calculation facilities are limited and there is little

354

12 Applications

time available for carrying out the control algorithm because of the attention
called for by other operations. It is thus an excellent platform for implanting
precalculated GPC in industrial elds.
From all the possible loops that could be controlled the results obtained
in certain situations will be shown. These are: control of the cold water ow
F T2 with valve V5 , control of the output temperature of the heat exchanger
T T4 with valve V8 , control of the tank level LT1 with the cold water ow by
valve V5 and control of the tank temperature T T5 with the resistance.
12.2.3 Flow Control
The control of the cold water ow has been chosen as an example of regulating a simple loop. Because all the water supplied to the plant comes from
only one pressure group, the variations affecting the hot water ow or the
cold water ow of the heat exchanger will affect the cold water ow as disturbances. Regulating the cold water ow is not only important as a control
loop but it may be necessary as an auxiliary variable to control the temperature or level in the tank.
The dynamics of this loop are mainly governed by the regulation valve.
This is a motorized valve with a positioner with a total open time of 110 seconds, thus causing slow dynamics in the ow variation. The ow behaviour
will approximate that of a rst-order system with delay.
First the parameters identifying the process are obtained using the reaction curve, and then the coefcients of the GPC are found using the method
described in Chapter 5. In order to do this, working with the valve 70% open
(ow of 3.98 l/min), a step of up to 80% is produced, obtaining after the
transition a stationary value of 6.33 l/min. From the data obtained it can be
calculated that
K = 0.25

= 10.5 seconds

d = 10 seconds

when a sampling time T = 2 seconds is used, the parameters for the corresponding discrete model are:
a = 0.8265

b = 0.043

d=5

The control signal can easily be computed using the expression

u(t) = u(t 1) + (ly1 y(t + 5) + ly2 y(t + 4) + lr1 r(t))/K
where u(t) is the position of the valve V5 and y(t) is the value of the ow
F T2 . Using the approximation formulas (5.10) with = 0.8, the controller
gains result as:
ly1 = 2.931
ly2 = 1.864
lr1 = 1.067

12.2 Pilot Plant

355

Cold water flow (l/min)

100

200

300
Time (seconds)

400

500

600

100

200

300
Time (Seconds)

400

500

600

90.0

Valve V5

80.0

70.0

60.0

Fig. 12.13. Flow control

The behaviour of the closed loop when setpoint changes are produced
can be seen in Figure 12.13.
The disturbances appearing in the ow are sometimes produced by
changes in other ows of the installation and at other times by electrical
noise produced by equipment (mainly robots) located nearby. The setpoint
was changed from 4 to 6.3 litres per minute. The measured ow follows the
step changes quite rapidly (taking into account the slow valve dynamics)
with an overshoot of 13%. The manual valve at the cold water inlet was partially closed in order to introduce an external disturbance. As can be seen,
the controller rejected the external perturbation quite rapidly.
12.2.4 Temperature Control at the Exchanger Output
The heat exchanger can be considered to be an independent process within
the plant. The exchanger reduces the temperature of the recirculation water,
driven by the pump, using a constant ow of cold water for this. The way of
controlling the output temperature is by varying the ow of the recirculation
water with the motorized valve V8 ; thus the desired temperature is obtained
by variations in the ow. In brief, the heat exchanger is nothing more than a
tube bundle with hot water inside that exchanges heat with the exterior cold
water. It can thus be considered as being formed of a large number of rstorder elements that together act as a rst-order system with pure dead time

356

12 Applications
40

Temperature and Reference

38

Cooling water
valve closed

36

34

Resistor duty

32

cycle decreased

30

100

200

300

400
500
600
Time (Seconds)

700

800

900

1000

100

200

300

400
500
600
Time (Seconds)

700

800

900

1000

100

90

Valve V8

80

70

60

50

40

30

Fig. 12.14. Behaviour of the heat exchanger

(see Chapter 5). Thus the T T4 -V8 system will be approached by a transfer
function of this type.
Following the procedure used until now the system parameters which
will be used for the control law are calculated. Some of the results obtained
are shown in Figure 12.14. It should be born in mind that as the exchanger
is not independent from the rest of the plant, its output temperature affects,
through the tank, that of the input, producing changes in the operating conditions. In spite of this its behaviour is reasonably good.
The setpoint was changed from 38o C to 34o C. As can be seen in gure
12.14, the heat exchanger outlet temperature evolved to the new setpoint
quite smoothly without exhibiting oscillations. Two different types of external disturbances were introduced. First the manual valve of the refrigerating
cold water was closed for a few seconds. As was expected, the outlet temperature of the heat exchanger increased very rapidly because of this strong
external perturbation and then it was taken back to the desired value by
the GPC. The second perturbation is caused by decreasing the duty cycle of
the resistor in the tank, thus decreasing the inlet hot water temperature and
changing the heat exchanger operating point. As can be seen, the GPC rejects
almost completely this perturbation, caused by a change in its dynamics.

12.2 Pilot Plant

357

70
65

Temperature of the tank

60
55
50
45
40
35
30

200

400

600
800
Time (Seconds)

1000

1200

1400

200

400

600
800
Time (Seconds)

1000

1200

1400

100

Resistance (%)

80

60

40

20

Fig. 12.15. Evolution of the tank temperature

12.2.5 Temperature Control in the Tank

The next example chosen is also that of a very typical case in the process
industry: the temperature of the liquid in a tank. The manipulated variable
in this case is the duty cycle of the heating resistor.
The process has integral effect and was identied around the nominal
operating conditions (50o C). The following model was obtained:
G(s) =

0.41
e50s
s(1 + 50s)

The GPC was applied with a sampling time T = 10 seconds, = 1.2 and
N = 15. As in the previous case, the controller parameters were computed
by the formulas given in Chapter 5 for integrating processes. The results obtained are shown in Figure 12.15. A perturbation (simulating a major failure
of the actuator) was introduced. As can be seen, after the initial drop in the
temperature of the tank, caused by the lack of actuation, the control system
is able to take the tank temperature to the desired value with a very smooth
transient. A change in the setpoint from 50o C to 60o C was then introduced.
The temperature of the tank evolves between both setpoints without big oscillations.

358

12 Applications
90

85

Level (%)

80

75

70

65

60

200

400

600
800
Time (Seconds)

1000

1200

1400

200

400

600
800
Time (Seconds)

1000

1200

1400

70

60

Valve V5

50

40

30

20

10

Fig. 12.16. Evolution of the level

12.2.6 Level Control

Level is another of the most common variables in the process industry. In the
pilot plant, the level of the tank can be controlled by the input ows (cold
or hot water). In this example, the valve V5 will be used. The system was
identied around the nominal operating point (70%) by the reaction curve
method. The model transfer function is:
1.12 45s
e
G(s) =
1 + 87s
For a sampling time T = 10 seconds, the dead time results to be noninteger, thus the controller parameters must be calculated as shown in Section 5.2. The results obtained when working with a weighting factor of 1
and a prediction horizon of 15 are shown in Figure 12.16. The setpoint was
changed from 70 to 75%. As can be seen, the level of the tank moves smoothly
between both setpoints. An external perturbation was introduced to test the
disturbance rejection capabilities of the GPC. The perturbation consisted of
opening the hot water inlet and thus increasing the level of the tank. As can
be seen, the perturbation is rejected after a well-damped transient.
12.2.7 Remarks
The main objective of the control examples presented in this section was to
show how easily GPC can be implemented on a commercial distributed con-

12.3 Model Predictive Control in a Sugar Renery

359

trol system using the implementation technique presented in Chapter 5. The

GPC s were implemented without difculty using the programming language
(ITER) of the Integral Cube distributed control system.
Although comparing the results obtained by GPC with those obtained
using other control techniques was not one of the objectives, GPC has been
shown to produce better results than the traditional PID on the examples
treated. In all the processes, the results obtained by PID, tuned by the ZieglerNichols open-loop tuning rules, were very oscillatory, better results were obtained [203] after a long commissioning period where optimal PID parameters
were found. The commissioning of the GPC controllers was done in virtually
no time; they worked from the word go. The results obtained by GPC were
superior in all cases to the ones obtained by the PID controllers as reported
in [203].

12.3 Model Predictive Control in a Sugar Renery

This section shows an application of Precomputed GPC to a process in a sugar
renery. The implementation was carried out by the authors in collaboration

(Valladolid, Spain)
with the rm PROCISA. The renery is located in Penael
and belongs to Ebro Agricolas. The controller runs in a ORSI Integral Cube
control system, where the GPC has been included as a library routine which
can be incorporated in a control system as easily as the built-in PID routine.
The factory produces sugar from sugar beet by means of a series of processes
such as precipitation, crystallization, etc. The process to be controlled is the
temperature control of the descummed juice in the diffusion.
In order to extract the sugar from the beet it is necessary to dilute the
saccharose contained in the tuber tissue in water to form a juice from which
sugar for consumption is obtained.
The juice is obtained in a process known as diffusion. Once the beet has
been cut into pieces (called chunks) to increase the interchangeable surface,
it enters the macerator (which revolves at a velocity of 1 rpm) where it is
mixed with part of the juice coming from the diffusion process (see Figure
12.17). Part of the juice inside the macerator is recirculated to be heated by
means of steam and in this way it maintains the appropriate temperature for
maceration. The juice from the maceration process passes into the diffusor (a
slowly revolving pipe 25 m long with a diameter of 6 m) where it is mixed
with water and all the available sugar content is extracted, leaving the pulp
as a subproduct. The juice coming out of the diffusor is recirculated to the
macerator, from which the juice already prepared is extracted for the next
process.
For the diffusor to work correctly it is necessary to supply thermal energy to the juice during maceration. To obtain this objective, part of the juice
from the macerator (150 m3 /h) is made to recirculate through a battery of exchangers; within these the steam proceeding from the general services of the

360

12 Applications
Heat exchangers
Control
valve
Steam
TT Temperature

Beet chunks

Water

Product

Macerator

Diffusor
Pulp

Recirculation

Fig. 12.17. Diffusion process in a sugar renery

factory provides the heat needed to obtain optimum maceration. Therefore

the controller must adjust the steam valve (u) to achieve a determined return
temperature to the macerator (y).
The system response is seriously disturbed by changes in the steam pressure, which are frequent because the steam used in the exchangers has to be
shared with other processes which can function in a noncontinuous manner.
The process is basically a thermal exchange between the steam and the
juice in the pipes of the exchanger, with overdamped behaviour and delay
associated to the transportation time of the juice through pipes about 200
metres long. These considerations, together with the observation of the development of the system in certain situations, justify the use of a rst-order
model with delay.
A model was identied by its step response. Starting from the conditions
of 82.42 C and the valve at 57 %, the valve was closed to 37 % in order
to observe the evolution; the new stationary state is obtained at 78.61 C.
The values of gain, time constant and delay can easily be obtained from the
response (as seen in previous examples):
K=

0
C
82.42 78.61
= 0.1905
57 37
%

= 5 min

d = 1 min 45 s

However, it is seen that the system reacts differently when heated to

when cooled, the delay being much greater in the rst case. A similar test
changing the valve to 57 % again provides values of
K = 0.15

= 5 min 20 s

d = 4 min 50 s

Although an adaptive strategy could be used (with the consequent computational cost), a xed parameter controller was employed, showing, at the

12.3 Model Predictive Control in a Sugar Renery

361

110.0

90.0

70.0

50.0
Valve
Temperature
30.0

10.0
0.0

1.0

2.0

3.0
Time (hours)

4.0

5.0

6.0

1.0

2.0

3.0
Time (hours)

4.0

5.0

6.0

1.1

1.0

Steam pressure

0.9

0.8

0.7

0.6

0.5
0.0

Fig. 12.18. System response in the presence of external disturbances

same time, the robustness of the method when using the T-polynomial in the
presence of modelling errors. The error in the delay, which is the most dangerous, appears in this case. The following values of the model were chosen
for this
K = 0.18

= 300 s

d = 190 s

and sampling time of T = 60 s.

It should be noticed that there are great variations in the delay (that produced on heating is about three times greater than that on cooling), due to
which it is necessary to introduce a lter T (z 1 ) as suggested in [210] to increase the robustness.
The following gures show various moments in operating the temperature control. The behaviour of the controller rejecting the disturbances (sudden variations in the steam pressure and macerator load) can be seen in Figure 12.18. On the other hand, Figure 12.19 shows the response to a setpoint
change in the juice temperature.
It should be emphasized that this controller worked satisfactorily and
without interruption until the end of the years campaign, being handled
without difculty by the plant operators who continuously changed the values of the model throughout the operation time. Following many operational
days the operators themselves concluded that a satisfactory model was given

362

12 Applications
80.0
79.0

Temperature

78.0
77.0
76.0
75.0
74.0
73.0
72.0
0.0

0.5

1.0
Time (hours)

1.5

2.0

0.5

1.0
Time (hours)

1.5

2.0

90.0

Control valve

70.0

50.0

30.0
0.0

Fig. 12.19. Setpoint change

by:

K = 0.25

= 250 s

d = 220 s

= 0.1

The sampling time was set to 50 seconds and the following robust lter was
used
T (z 1 ) = A(z 1 )(1 0.8az 1 )
with a equal to the discrete process pole.

12.4 Olive Oil Mill

This example describes the application of a predictive controller that deals
with measurable disturbances in the extraction process in an olive oil mill.
The application focuses on the thermal part of the process, where the raw
material is prepared for the mechanical separation. The system under consideration can be viewed as composed of several changing level stirred tanks.
The example shows the development of the controller based upon a model
obtained from rst principles combined with experimental results and validated with real data. Strong disturbances and large time delays appear in the
process, so predictive control strategies have been tested under simulation
and have been implemented on the real plant. A study about the consideration of different models for the estimation of measurable disturbances along

12.4 Olive Oil Mill

363

Clean olives
Heating water
Addition water
LT

LT

Extraction

LT

Paste pump
Mill

Thermomixer

Olive
Oil

By-product

Fig. 12.20. Process

the prediction horizon shows that good performance can be obtained by use
of an appropriate model. A new idea that can improve periodic disturbance
rejection in Model Based Predictive controllers is also presented.
The process is composed of several operations: recollection and reception
of raw material (olives), washing, preparation, extraction, and storage of the
produced oil. Figure 12.20 shows the most important phases of the process:
preparation and extraction.
The preparation phase is crucial for the whole process; it consists of two
subprocesses. The rst is olive crushing by a special mill, whose objective
is to destroy the olive cells where oil is stored. The second aims at homogenising the paste by stirring it while its temperature is kept constant at a
specic value (around 35 o C). This is performed by a machine called a thermomixer, which homogenises the three phases of the paste (oil, water and a
by-product) while it exchanges energy with surrounding pipes of hot water.
This is done to facilitate oil extraction in the following process: mechanical
separation in the decanter. This case study is focused on thermomixer control
since homogenisation is crucial in the entire process. Bad operational conditions in the thermomixer can dramatically reduce the quality and quantity of
the nal product.
There are three main obstacles that appear when trying to maintain the
optimal operating conditions in the thermomixer. The rst is the existence of
large delays (around one and a half hours) because of the thermal nature of
the process. The predictive controller treats the delays in a convenient way.
The second obstacle is caused by the on-off mechanism feeding the paste,
so the inlet paste ow does not reach a constant value. These changes introduce periodic variations in level and therefore temperature changes since
the quantity of product inside the machine varies. As the level can be eas-

364

12 Applications

ily measured, it can be considered as a measurable disturbance and hence

can be taken into account by the predictive algorithm as a feedforward action. The third difculty usually takes place at the beginning or the end of the
campaign. The process is frequently interrupted because of the heterogeneity
and low quantity (and often, low quality) of raw material. When the process
is stopped and restarted, the temperature inside the thermomixer increases
rapidly.
12.4.1 Plant Description
The system considered corresponds to a thermomixer, whose main objective
is to homogenise the three phases of the paste (oil, water and the by-product)
and keep it at a desired temperature to facilitate oil extraction. Heating the
paste is achieved by means of hot water circulating through a jacket. The
machine is divided into different (usually two, three or four) tanks or bodies, each with revolving blades to facilitate homogenising. The bodies are
composed of semicylinders about 3 metres long with a diameter of 1 metre.
Paste is dropped over one side of the rst body and pushed by the revolving blades which make the paste fall down to the second body through the
overow and so on. The existence of several bodies allows a gradual temperature increment along the thermomixer, since abrupt changes in paste temperature would affect the quality of the end product. Each body has its own
water jacket. The water circuit is connected in parallel, so the same quantity
of water ows for each jacket. We can only set the total water ow (which is
shared equally amongst all the bodies).
The paste is heated to facilitate mixing since the paste turns more liquid
as the temperature rises. However, there exists a maximum temperature at
which olive oil loses quality (avour, fragrance, etc.) due to the oxidation
process and the loss of volatile components. The heating water comes from
a boiler that supplies hot water to several processes in the factory, so it is
affected by load changes and is therefore another disturbance.
Therefore, the outlet temperature presents oscillations at the frequency
of level variations with changes produced by heating water variation. The
controller must be able to reduce the effect of these disturbances as much
as possible. Level and water temperature can easily be measured and future
level evolution can be predicted as shown later. Figure 12.21 shows level evolution in the last body and the random variations of inlet water temperature.
It is clear that a more efcient performance could be obtained by a good design in the feeding system that keeps the level constant; this is the kind of
level control system that is installed in most olive oil mills for cost reasons.
Therefore, to address the current problem, this disturbance has been considered as something external to the proposed control solution.

12.4 Olive Oil Mill

365

Measured level in the third body (%)

100
95
90
85
80
75
70

0.5

1.5

2.5

3.5

4.5

4.5

Measured inlet water temperature (C)

80
75
70
65
60
55
50

0.5

1.5

2.5
Time(hours)

3.5

Fig. 12.21. Measurable disturbances

12.4.2 Process Modelling and Validation

The attainment of the process model is described in [32], where thermodynamical equations together with input-output data allows construction of a
nonlinear simulator that can be used to test the controller. The model was
validated using real data obtained from an olive oil mill located in Malaga
(Spain). These data were used to estimate many of the parameters that appear in the model which are not perfectly known, since they depend on several circumstances: type and moisture of olives, soil in the heating circuit,
etc.
A linear model must be developed to design a simple linear predictive
controller that will run on a Programmable Logic Controller (PLC) with low
computational capabilities. The linear model is obtained by performing tests
in the nonlinear model, where all manipulated variables can be changed independently to see their inuence on temperature behaviour. The models
needed for control must predict correctly the nal paste temperature (the
paste that leaves the last body of the thermomixer) as a function of hot water
ow, level, and temperature of the heating water.
For the CARIMA model, numbers that give temperature as a function of
the hot water ow are

366

12 Applications
45
Real
Linear model

Temperature (C)

40

35

30

25

20
10

11

12

13

14

15
Time (hours)

16

17

18

19

20

Fig. 12.22. Validation linear prediction model

B(z 1 ) = 4.315 105 z 5

A(z 1 ) = 1 0.92z 1
temperature with respect to level:
B(z 1 ) = 0.0028z 3
A(z 1 ) = 1 0.8825z 1
and temperature with respect to water temperature
B(z 1 ) = 0.0093z 3
A(z 1 ) = 1 0.926z 1
The model used in simulation has a sample time of 100 seconds. Figure
12.22 shows the predictions of temperature using a linear model with a prediction horizon of seven, which is the one that will be used by the predictive
controller, compared with real data. This corresponds to a new set of data not
used to build the model.
12.4.3 Controller Synthesis
The control objective is to maintain the operating conditions in the thermomixer, that is equivalent to keep the temperature of the last body as constant as possible in spite of disturbances (level and variations in hot water
temperature). The manipulated variable is the hot water ow.
The process is characterized by big dead times in temperature dynamics.
Moreover, the effect of disturbances, mainly level, on the controlled variable

12.4 Olive Oil Mill

367

shows fast dynamics (mainly at high production rates). This fact makes disturbance rejection more difcult and implies the necessity of a control strategy that eliminates the effect of level disturbance, at least in the nominal case.
MPC can be an interesting candidate to control this system and disturbance rejection capabilities can be improved by the estimation of the future
disturbance values in the controlled variable prediction model. That is, the
current value of the disturbance is known at the sampling time, but its future
evolution along the horizon can also be estimated since it affects the process
predicted output (temperature). This is a slight improvement with respect
to standard MPC algorithms, which consider that disturbances are kept constant (and equal to their current value) in the future. The information that
provides this future evolution is very important in this case, allowing the
controller to anticipate its inuence on the process output.
In this application, as the main disturbance acting on the output (level)
exhibits a predictable behaviour, the control law is calculated considering an
Auto-Regressive (AR) model of disturbance. Periodic disturbance rejection
is also treated in [29]. Therefore the AR model is not used to estimate the
actual values of the disturbances, but utilized as an instrument to improve
the prediction of future values of the disturbances. So the prediction of the
controlled variable must also be enhanced. The identied AR polynomial is:
A(z 1 ) = 1 1.589z 1 + 1.696z 2 1.178z 3 + 0.697z 4
When a disturbance model is included in the prediction, the control law
that minimizes the general cost function
J=

p


[
y (t + i) w(t + i)] +

i=1

m1


[u(t + i)]

(12.5)

i=0

is given by (see [31])



u = (GT G + I)1 GT w (fu + fd )

(12.6)

where

u is the vector of future control action increments;

fu is the calculated free response without disturbances;
w is the reference trajectory and
fd is the estimated value of the free response due to measurable disturbances and can be calculated using the following equation
fd (t + i) =

i1

k=0

+ i k) +
gk d(t

N


gk d(t + i k)

(12.7)

k=i

where gk are the samples of the truncated step response and d(t) is the
increment of the disturbance signal at time t and the AR prediction is given
by

368

12 Applications

=
d(t)

N


ai d(t i)

(12.8)

i=1

The terms that depend on future values of the disturbances are separated
from the directly measured ones. The rst part in the sum is assumed to be
zero in standard MPC formulations. This means that the polynomial AR is
equal to 1 z 1 , which is equivalent to the assumption that future values of
disturbance equal its current value.
Therefore, the control algorithm (MPC with Auto-Regressive model for
the measurable disturbances or MPC - AR) is reduced to the following:
1. Compute the free response as in a basic MPC algorithm.
2. Add the term due to measurable AR disturbance as shown in (12.7 - 12.8).
3. Calculate the control law with Equation (12.6).
The controller presented here has been devised to be implemented in
low-cost control equipment, as is the case of a PLC. This can be done since
only part of the algorithm is calculated in real time. The most demanding
part of MPC (the optimization that corresponds to Equation (12.6)) is done
beforehand, since the model is known and the only part to be calculated at
every sampling time is the free response of the system. Should the model
change with time, the control law parameters will be updated by a higherlevel routine running on the SCADA computer.
12.4.4 Experimental Results
Several tests have been done on the real plant. The experiments have been
carried out in different situations since it is very difcult to obtain the same
operational conditions in different days because of the continuously changing raw material. Several tunings of the controller (using different weighting
factors and different disturbance models) have been tested. The controller
runs on a PLC, where the control law can be switched to an existing PID controller for comparison. The PID was tuned manually by trained personnel:
see [32].
Table 12.4 shows the root mean square error of the temperature, comparing a PID with feedforward to the proposed MPC. Results with changes in
temperature setpoint are also presented.
Table 12.4. Fits obtained with PID and proposed MPC in real plant
Error (o C) PID MPC-AR
Regulation 2.1
0.5
Tracking 2.9
0.9

Figure 12.23 shows the controlled temperature in the last body when the
setpoint changes in real tests. MPC - AR is able to reach the set point faster

12.4 Olive Oil Mill

    

40

3 body temperature (C)

369

38

36

34

32

30

10

11

12

13

10

11

12

13

14

15

16

17

18

19

15

16

17

18

19

80

Valve (%)

70
60
50
40
30
20

   
 

40

3 body temperature (C)

14
Time (hours)

=0.001

38

36

34

32

30
12
90

13

14

15

16

17

18

19

13

14

15

16
Time (hours)

17

18

19

80

Valve (%)

70
60
50
40
30
20
12

Fig. 12.23. Changes in the setpoint. PID control (two upper graphs) and MPC-AR
control (two lower graphs).

than PID, which ever continuously oscillates around the setpoint. The control
action of PID is obviously slower, and it cannot be accelerated much more
without losing stability.
Figure 12.24 presents two different tests of the designed MPC in the real
plant, performed during the intermediate dates of the campaign, in which

370

12 Applications

Temperature and set point (C)

MPC, aggressive control

40
38
36
34
32
30
19

20

21

22

23

24

Temperature and set point (C)

MPC, smooth control

40
38
36
34
32
30
0

3
Time (hours)

Fig. 12.24. MPC with different values for

the process is often stabilized. The rst corresponds to a controller with a

small value for the control weighting factor () and no measurable disturbances considered. The second shows the same controller with a bigger value
of , giving a slower control. The parameter denes the aggressiveness of
the controller. If is set to a very small value, the closed-loop system behaves faster, but it loses robustness. In this application, the ts achieved in
both cases are similar, because a faster control induces more high-frequency
noise. Although the t in the rst case is slightly better, it is preferable to
have a more robust and smoother control.
The biggest benets of the inclusion of AR model for measurable disturbances can be obtained at the end of the olive season, when the operational
conditions are not stationary and effect level evolution. In the rst case (Figure 12.25, top graph) the behaviour of the temperature is not good, with considerable oscillations around the setpoint, although it is the same controller
that has shown good performance during the intermediate dates of the campaign. The performance is clearly improved with the AR model (Figure 12.25,
bottom graph), under the same conditions (similar evolution of level and water temperature), showing that the proposed algorithm is a viable solution to
the problem.
Plant results have shown that an MPC considering the prediction of future level variations can be a recommendable solution for the problem that
exists in the real plant.

12.5 Mobile Robot

371

Temperature and set point (C)

MPC controller, RMS=1.35

40
38
36
34
32
30
0

0.5

1.5

2.5

3.5

3.5

Temperature and set point (C)

MPCAR controller, RMS=0.60

40
38
36
34
32
30
0

0.5

1.5

2
Time (hours)

2.5

Fig. 12.25. MPC at the end of olive season

12.5 Mobile Robot

This section shows the application of NMPC to mobile robot navigation. This
implementation tries to solve one of the main problems in the development
of autonomous mobile robots: the problem of path tracking in an environment with unexpected obstacles. MPC is a suitable technique to apply to this
problem because the objective is to drive future outputs (robot position and
heading) close to the desired values in some way, bearing in mind the control
activity required to do so.
If the control signal is constrained and the system model is nonlinear
the MPC results in a very complex and time-consuming problem. A neural
network is used to solve the problem, as described in Section 9.4.
12.5.1 Problem Denition
The problem of driving a mobile robot to follow a previously calculated desired path has been addressed in [77] where an NN is used to solve the NMPC
problem online with the following objective function

372

12 Applications

J1 (N1 , N2 , Nu ) =

N2


[Y (t + i|t) Yd (t + i)]2

i=N1

Nu


1 ([r (t + i 1)]2 + [l (t + i 1)]2 ) +

i=1

Nu


2 [r (t + i 1) l (t + i 1)]2

i=1

where Y (t + i|t) = {
x(t + i|t), y(t + i|t)} is an i step prediction of the robot
position made at instant t; r and l are the right and left angular velocities
of the two driving wheels, which are the control variables; and 1 , 2 and
are constant weighting factors. The rst term in J penalizes the position
error; the second term penalizes the acceleration and the third penalizes the
robot angular velocity. These last two terms ensure smooth robot guidance.
When unexpected static obstacles are taken into account, a new term
must be added to the cost function in order to penalize the proximity between the robot and the obstacles, which are detected with an ultrasound
proximity system placed on board the mobile robot. Therefore, the function
to be minimized now is:

J(N1 , N2 , Nu ) = J1 (N1 , N2 , Nu ) +

N
FO


N2


j=1 i=N1

[distf (Y (t + i|t), F Oj )]2

The new term is a potential function term, where distf () is a measurement in t + i of the distance between the robot and a xed obstacle F Oj ,
which is considered to have a polygonal geometry in the plane. A more precise description of this function is presented here. A block diagram of the
system is shown in Figure 12.26. Notice that the consideration of unexpected
obstacles makes the objective function not quadratic and thus the computational burden is increased.
In the following N1 and N2 will be considered to be N1 = d + 1 and
N2 = N , and Nu will be given a value of N2 d. So the controller has only
one free parameter N . The predictive problem, formulated under these circumstances, has to be solved with numerical optimization methods, which
are not acceptable for real-time control. The controller is implemented using
a neural network scheme, which allows real time.
12.5.2 Prediction Model
For an MBPC formulation, a model of the mobile platform is needed to predict
the future positions and headings of the robot. As a testbed for the experiments, a TRC LABMATE mobile robot has been used.

12.5 Mobile Robot

Predictive controller

Past control
actions

Future predicted

Non-linear

^
outputs Y(t+i)

model

Environment
sensor
information

Future obstacle
positions prediction

Future
references
Yd (t+i)

373

+
-

Future

U(t)

Minimization

errors

Future control
actions
Objetive
function

Mobile
Robot

Y(t)

constraints

Fig. 12.26. The predictive controller scheme

A model of the LABMATE mobile robot, which takes into account the dead
time produced by communication with the host processor, was obtained by
using kinematic equations and identication tests. The following kinematic
model (which corresponds to a differential-drive vehicle) is used for computing the predictions
(k + 1) = (k) + AT
V
(sin((k) + AT ) sin((k)))
A
V
y(k + 1) = y(k) (cos((k) + AT ) cos((k)))
A
r (k 1) l (k 1)
A=R
2W
r (k 1) + l (k 1)
V =R
2
where x, y and are the position and heading of the robot in a xed
reference frame (see Figure 12.27), T is the sample interval and W is the halfdistance between wheels, which value has been estimated to be 168 mm (Figure 12.27). V is the linear velocity of the mobile robot, A is the steering speed,
and r (k 1), l (k 1) and R are the right and left wheel angular velocities
(which are considered to be constant for each sample interval) and the wheel
radius, respectively. In the case of a linear trajectory (A = 0), the equations
of motion are given by:
x(k + 1) = x(k) +

374

12 Applications
Yf

Yr
Xr

PG
2
W
x

Xf

Fig. 12.27. Reference frame

(k + 1) = (k)
x(k + 1) = x(k) + V T cos (k)
y(k + 1) = y(k) + V T sen (k)
Using the maximum acceleration value, the velocities of both wheels
have been considered to be constant for each sample period.
12.5.3 Parametrization of the Desired Path
The reference path is given to the MPC controller as a set of straight lines and
circular arcs. The MPC approach needs the desired positions and headings of
the mobile platform at the next N time instants. So, given the current position
and heading of the robot, it is necessary to parameterize the desired path for
the next N periods of time to calculate the N future path points desired. As is
shown in Figure 12.28, the desired point for the current instant (Xd (k), Yd (k))
is obtained rst. It is located at the intersection between the desired path and
its perpendicular, traced from the actual robot position (Xr (k), Yr (k)). The
next N points are spaced equally on the path, with a separation between
them of S, which is a design parameter.
12.5.4 Potential Function for Considering Fixed Obstacles
As stated earlier, the xed obstacles are considered to have polygonal geometry in the plane, and the surfaces of the obstacles are considered to be
perpendicular to the moving plane of the mobile robot.
Two different potential functions have been used for the polygonal geometry case: one for the convex polygon and the other for the concave polygon:

12.5 Mobile Robot

375

(k
y(k)

d (k)
90

y (k)

. . . .

y d (k+1)

x(k)

. . . . . . .x (k+N)

x (k) xd(k+1)
d

Fig. 12.28. Parametrization of the desired path

potential function for a convex polygon [94]. A convex region will be described by a set of inequalities
g(x) 0, g Lm , x Rn
where L is the set of linear functions and n is the space dimension (in this
case n = 2). The function
f (x) =

N
SF


gi (x) + |gi (x)|

i=1

is zero inside the convex region and increases linearly out of it, as the
distance to the frontier is augmented. N SF is the number of segments of
the obstacle frontier. The following potential function is used
1

pcvx (x) = [ + f (x)]

=
+

N
SF
i=1

(gi (x) + |gi (x)|)

where is a small constant that limits the value of pcvx (x) inside the convex region; pcvx (x) reaches its maximum value 1 inside the region occupied by the obstacle, and decreases with the distance between the robot
and the obstacle. A graphic example of this function is shown in Figure
12.29, where two rectangular static obstacles are present in the proximity
of the robot.

376

12 Applications

0.5
0.4
0.3
0.2
0.1
12

-2
10

0
8

6
x

4
y

0
-2

10
12

Fig. 12.29. Convex regions potential function

potential function for a concave polygon. A concave region will be described by a set of inequalities
v(x) 0, v Lm , x Rn
The potential function used in this case is
pccv (x) =

1
+ gccv (x)

where is a small constant and gccv (x) is the minimum of the distances
between the robot position and every straight line that denes the obstacle frontier. This function has the same characteristics as pccx (x).
12.5.5 The Neural Network Approach
As was mentioned before, the minimization of the cost function J has to be
carried out by a numerical optimization method which requires too much
computation to be used in real time. A neural network solution is proposed,
which guarantees real time for the robot control.
The modules of the control scheme used in this work (see Figure 12.30)
are:
articial neural network controller. The NN architecture chosen here is a
Multilayer Perceptron, with one hidden layer (see Figure 12.31).
The input layer consists of twelve neurons. The rst two inputs correspond to the previous linear and angular velocities of the robot. The next

12.5 Mobile Robot

Global

377

Past control
actions

path

planning

Sensored enviromental
Future references at
information

absolute coordinates

Transformation of
future references into
relative coordinates

Input

Output

Vector

Vector

Future

referenc. Generation

U(t)

Mobile

Y(t)

robot

Generation

Position and
heading estimation

Fig. 12.30. Predictive neural network scheme for mobile robot navigation

V(k-1)

(k)

....

d (k)

....

(k-1)

V(k)

1/(k)

(k)

....

....

....

S0 (k)

S5 (k)
S6 (k)

Fig. 12.31. Neural network scheme

three inputs are associated to the parameterization of the desired trajectory over the prediction horizon. To reduce the number of inputs, the
parameters given to the network are the distance d from the robot guide
point to the path, the angle between the robot heading and the path
orientation and an average of the inverse of the curvature of the future
desired points (1/) (see Figure 12.28). The last seven inputs correspond
to the distances measured directly by the ring of sonar sensors. This fact
avoids the high-level process that usually has to be carried out with sen-

378

12 Applications
Global path
planning

SIMULATED
SIMULATED

Robot position
and heading

PREDICTIVE
OUTPUTS
CONTROLLER

Obstacle position

SONAR
Sonar model

MEASURES

A.N.N.
+

A.N.N. weights

Modification of

OUTPUTS

BACKPROPAGATION

ERROR

MODULE

Fig. 12.32. ANN supervised training scheme

sor data to provide useful environment information, which is a difcult

and high computation time-consuming phase. The output layer consists
of two nodes which correspond to the control commands (the linear and
angular robot velocities).
input vector generation module. A symmetry analysis is made here to
reduce the number of training patterns needed to provide good performance of the neural network controller [78]. Also a normalization is made
which leads to better performance at the NN training stage.
output vector generation module. This performs the inverse transformation of that made at the symmetry input module when required.
reference path coordinates transformation module. The desired path coordinates are transformed from a global reference system to a local reference system, attached to the mobile robot. This avoids the use of additional NN inputs for the robot position and heading, which are implicitly
given to the NN in the reference path.
past control actions. These are needed for NN to consider the delay time
of the robot system.
sonar range measurements. Their measurements are directly used as inputs for NN. ANN learns from the input patterns set, where different situations of static obstacles are present. Thus, there is no need for a high-level
sonar measurement preprocess, which guarantees real time.

12.5.6 Training Phase

The controller has been trained using a classic supervised training scheme as
the backpropagation algorithm (see Figure 12.32). The training patterns set

12.5 Mobile Robot

379

have been obtained from an ofine simulation system. For the minimization
phase a Powell iterative algorithm has been used, where constraints on the
control variables are considered. Also, the sonar system measurements have
been simulated using a sonar model where the objects in the environment
are described as a set of geometric primitives such as planes, cylinders, edges
and corners.
12.5.7 Results
The proposed control structure has been tested with the L ABMATE mobile
robot.
The sampling interval T was given a value of 2 s. The value of N chosen
for the MPC was made equal to 7, thus N1 , N2 and Nu were given the values
2, 7 and 6, respectively, and the weighting factors were given the following
values: 1 = 35, 2 = 5 and = 0.5. The maximum and minimum linear and
angular velocities were given the following values, respectively: 0 m/s, 0.4
m/s, -20 o /s and 20 o /s. For s, a value of 0.15 m was chosen, which leads to
an average linear robot velocity of 0.25 m/s.
Figure 12.33 shows some of the experimental results obtained in the laboratory when applying the proposed algorithm to the LABMATE mobile robot.
Although it is drawn as straight lines in the gures, the real environment
includes laboratory objects such as chairs, tables, etc. This shows the robustness of the neural controller as it has been trained with only a set of geometric
primitives such as planes, corners, etc.
Figure 12.33(a) shows the desired trajectory and the real trajectory through
the laboratory where unexpected static obstacles have been positioned. It is
important to notice how the mobile robot returns to the reference path after
an unexpected obstacle has been avoided.
Figure 12.33(b) shows an experiment where an unexpected obstacle is
situated in the path that the robot must follow to avoid a previous unexpected obstacle. Again, the controller performance is quite good. Finally, Figure 12.33(c) shows another test for a path where small curvature radii are
specied. The tracking error observed in Figure 12.33(c) is due to saturation
in the angular velocity and to the penalization chosen for the control actions.
More details about the implementation of the controller can be found in
[78], [79].

380

12 Applications
4.0

2.0

Y (m)

0.0

-2.0

-4.0

-6.0
-2.0

0.0

2.0

a)

4.0

6.0

X (m)

3.0

Y (m)

2.0

1.0

0.0

-1.0
0.0

2.0

b)

4.0

6.0

X (m)

1.0

Y (m)

-1.0

-3.0

-5.0
-3.0

-1.0

c)

1.0

3.0

5.0

X (m)

Fig. 12.33. Experimental results of the neural predictive controller for mobile robot
navigation