# Constrained Fuzzy Predictive Control Using Particle Swarm Optimization

Authors: Oussama Ait Sahed; Kamel Kara; Mohamed Laid Hadjili
Journal: Applied Computational Intelligence and Soft Computing
 (2015)

---

## Abstract

A fuzzy predictive controller using particle swarm optimization (PSO) approach is proposed. The aim is to develop an efficient algorithm that is able to handle the relatively complex optimization problem with minimal computational time. This can be achieved using reduced population size and small number of iterations. In this algorithm, instead of using the uniform distribution as in the conventional PSO algorithm, the initial particles positions are distributed according to the normal distribution law, within the area around the best position. The radius limiting this area is adaptively changed according to the tracking error values. Moreover, the choice of the initial best position is based on prior knowledge about the search space landscape and the fact that in most practical applications the dynamic optimization problem changes are gradual. The efficiency of the proposed control algorithm is evaluated by considering the control of the model of a 4 × 4 Multi-Input Multi-Output industrial boiler. This model is characterized by being nonlinear with high interactions between its inputs and outputs, having a nonminimum phase behaviour, and containing instabilities and time delays. The obtained results are compared to those of the control algorithms based on the conventional PSO and the linear approach.


---

## Body

## 1. Introduction

In the late seventies, a new class of advanced control algorithms, grouped under the denomination Model Predictive Control (MPC), has emerged. The MPC strategy is based on the use of an explicit model to predict the process future behaviour over a finite horizon and then to compute a sequence of future control actions by minimizing a given cost function. The MPC is known to be a very powerful control strategy for many industrial processes [1–3]. Its attraction is due to its ability to handle complex control problems which involve multivariable process interactions, constraints in the system variables, nonminimum phase behaviour, and variable or unknown time delays. Linear MPC (LMPC) techniques, which use a linear model for the process, have been successfully used in many industrial applications, such as refining, chemicals, petrochemicals, air and gas processing, and food processing industries [4–6].Although LMPC algorithms provide satisfactory performance in many applications, against highly nonlinear process, severe degradation in the control performance can occur, unless of course the operating conditions are very close to the steady state around which the model is linearized [7]. Since most practical processes are nonlinear by nature, new efficient nonlinear MPC (NMPC) techniques have to be derived to incorporate nonlinearities and ensure higher control performance. The introduction of nonlinear models will enhance the overall controller performance. However, the relatively simple optimization problem in LMPC algorithms will be transformed to a nonlinear and a nonconvex optimization problem which requires complex and time consuming procedures [8–10]. Thus, the challenge is to derive an efficient optimization algorithm that allows the real-time implementation of the control law.An important step in the design of NMPC algorithms based controllers is the construction of the explicit nonlinear model of the system. The model choice is a critical task; the controller efficiency and performance depend extremely on the accuracy of the used model. In many cases it is even impossible to obtain a suitable model for the underlying process due to its complexity or lack of knowledge of its critical parameters. Different approaches that can provide satisfactory models exist; one which is widely used is the fuzzy logic technique [11]. Fuzzy inference systems (FIS) are universal approximators capable of approximating any continuous function with a certain level of accuracy [12]. Takagi-Sugeno (TS) models, a subdivision of fuzzy models, are particularly suitable for NMPC algorithms [13]. These models are able to express the dynamic nature of systems with characteristics of randomness, large delay time, and strong nonlinearity [14, 15].Another important factor in driving the NMPC based controllers is the resolution of the relatively difficult and time consuming nonlinear and nonconvex optimization problem. Seeing that the analytical approach is impractical, several other suboptimal solutions were proposed [7, 16–18]. Some of them are based on numerical optimization algorithms, which are known to have a slow convergence rate and could easily be trapped in local minima. Other more interesting approaches are the metaheuristic based optimization algorithms, such as the powerful but computationally complex genetic algorithms (GA) and the particle swarm optimization strategy. PSO algorithms require much less computational effort than the GA and their implementation is relatively easy [19]. A lot of papers have considered the PSO to solve the optimization problem of NMPC [20–23]. In these papers, a large number of particles and iterations were used which increases the computational burden of the control algorithm. Recently, several PSO algorithms have been proposed in the literature to solve different optimization problems [24–26]. In most of these algorithms, improving the accuracy and the convergence proprieties of the algorithm is obtained at the cost of increasing its computing complexity. For example, in [26] more than one equation is used to update the particles positions at the same time. Obviously, using several update equations will increase the algorithm computing time and will limit its use in real time applications.In this work, a PSO based constrained fuzzy predictive controller is proposed. The aim is to develop a simple and efficient controller capable of handling the relatively complex optimization problem with minimal computational time. This can be achieved using a reduced population size and a small number of iterations. The proposed control algorithm uses the normal distribution law to distribute the initial particles positions within the area around the best position. The radius limiting this area is adaptively changed according to the tracking error values. Moreover, the choice of the initial best position is based on prior knowledge about the search space landscape and the fact that in most practical applications the dynamic optimization problem changes are gradual. The algorithm efficiency is evaluated by considering the control of the model of a 4 × 4 Multi-Input Multi-Output (MIMO) industrial boiler [27, 28]. This model is characterized by being nonlinear with high interactions between its inputs and outputs, having a nonminimum phase behaviour, containing instabilities and time delays, and including constraints on its variables, disturbances, and uncertainties.The paper is organized as follows. In Section2 the design stages of Takagi-Sugeno fuzzy models are presented. Section 3 presents and describes the fuzzy predictive control problem. The PSO algorithm, used to solve the fuzzy MPC optimization problem, is introduced in Section 4. Section 5 deals with the proposed control algorithm and Section 6 gives the simulation results. Finally, conclusions are drawn in the last section.


## 2. Takagi-Sugeno Fuzzy Modelling

Fuzzy inference systems (FIS) are universal approximators capable of approximating any continuous function with certain level of accuracy [12]. The FIS are based on the so-called IF-THEN rules and are classified into two main categories: Mamdani models and Takagi-Sugeno models. The rules of TS fuzzy models, used in this work, have the following form:(1)
R
j
:  IF
F
is
A
j
THEN
y
j
=
θ
j
T
X
,where F
=
{
y
t
,
y
t
-
1
⋯
y
t
-
n
a
+
1
,
u
t
,
u
t
-
1
⋯
u
(
t
-
m
a
+
1
)
} is the set of n
a
+
m
a premise values, X
T
=
{
y
t
,
y
t
-
1
⋯
y
t
-
n
r
+
1
,
u
t
,
u
t
-
1
⋯
u
t
-
m
r
+
1
,
1
} is the set of m
r inputs and n
r outputs values used in the consequent regressors, A
j
=
{
A
1
,
j
,
A
2
,
j
⋯
A
n
a
,
j
,
B
1
,
j
,
B
2
,
j
⋯
B
m
a
,
j
} is the set of membership functions associated with the antecedents of the j
t
h rule, and θ
j
T
=
{
a
j
,
1
,
a
j
,
2
⋯
a
j
,
n
r
,
b
j
,
1
,
b
j
,
2
⋯
b
j
,
m
r
,
c
j
} is the parameter vector of the j
t
h submodel while y
j is its output (1
≤
j
≤
r; r is the number of fuzzy rules).It is clear, from (1), that the Takagi-Sugeno fuzzy models have fuzzy propositions only in their antecedents while their consequences are linear functions of the antecedents or their variables.Using input/output data extracted from a nonlinear process, a fuzzy model can be constructed to mimic the process behaviour, by defining the parameter matrix as(2)
Θ
=
θ
1
,
θ
2
⋯
θ
rand the normalized membership grade vector as(3)
μ
T
=
μ
1
,
μ
2
⋯
μ
r
.The inferred output of the TS fuzzy model is given by(4)
y
^
t
+
1
=
Θ
μ
T
χ
.Using the input/output representation, the TS model given by (4) can be rewritten as follows:(5)
y
^
t
+
1
=
∑
p
=
1
n
a
a
¯
p
t
y
t
-
p
+
1
+
∑
p
=
1
m
a
b
¯
p
t
u
t
-
p
+
1
+
c
¯
t
,where a
¯
p
t
=
∑
l
=
1
r
μ
l
t
a
l
,
p, b
¯
p
t
=
∑
l
=
1
r
μ
l
t
b
l
,
p, and c
¯
t
=
∑
l
=
1
r
μ
l
t
c
l.The global output of the fuzzy model can be written as(6)
y
^
t
+
1
=
∑
p
=
1
r
μ
p
t
y
p
t
+
1
=
μ
1
t
y
t
,
…
,
μ
1
t
u
t
-
m
r
+
1
,
…
,
μ
r
t
θ
~
=
φ
T
t
θ
~
,where y
p
(
t
+
1
) is the output of the pth submodel and(7)
θ
~
T
=
θ
1
T
θ
2
T
⋯
θ
r
T
.Assuming that a set of N input-output data pairs u
t
,
y
t is available and defining L
=
m
a
x
m
r
,
n
r, a regression matrix which has the following expression can be constructed:(8)
Φ
=
φ
T
L
φ
T
L
+
1
⋮
φ
T
N
-
1
.The vector θ
~ can be calculated by solving the following least square problem:(9)
Y
=
Φ
θ
~
,where Y
=
y
L
+
1
,
…
,
y
(
N
)
T.The resolution of (9), which is the last step in TS fuzzy model construction, deals with the definition of the consequences part of the rules. However, the antecedent, which involves in part the number, the position, the shape, and the distribution of the membership functions, must be first selected. A lot of techniques for the construction of fuzzy models from input-output data are described in the literature; in Espinosa et al. [16] different approaches are well explained.


## 3. Fuzzy Predictive Control Principle

We can always distinguish two common features present in any given MPC based algorithm. The explicit model used to predict the future process behaviour and the optimization problem from which a control sequence is derived.The main steps of the MPC strategy are given below:(1)
The constructed model is used to predict the future behaviour of the process over a given prediction horizon.
(2)
The reference trajectory must be defined over the prediction horizon.
(3)
The control sequence is obtained by minimizing a given cost function.
(4)
Only the first element of the control sequence is applied on the system.The previous steps are repeated, at every sampling time, according to the receding horizon idea. These steps are illustrated by the block diagram of Figure1.Figure 1
Basic structure of MPC technique.The optimization problem to be solved at every sampling time usually has the following form:(10)
min
⁡
Δ
u
t
J
u
t
,
y
^
t
,
w
t
=
∑
j
=
N
1
N
2
y
^
t
+
j
∣
t
-
w
t
+
j
T
Q
y
^
t
+
j
∣
t
-
w
t
+
j
+
∑
j
=
1
N
u
Δ
u
T
t
+
j
-
1
R
Δ
u
t
+
j
-
1
subject to:
Δ
u
t
+
j
-
1
=
0
,
for
j
>
N
u
y
min
⁡
<
y
^
t
+
j
∣
t
<
y
max
⁡
for
N
1
≤
j
≤
N
2
u
min
⁡
<
u
t
+
j
<
u
max
⁡
for
1
≤
j
≤
N
u
Δ
u
min
⁡
<
Δ
u
t
+
j
-
1
<
Δ
u
max
⁡
for
j
≤
N
u
,where Δ
u
(
t
)
=
u
(
t
)
-
u
(
t
-
1
) is the control increment, J is the cost function, N
1 and N
2 are, respectively, the minimum and maximum prediction horizons, N
u is the control horizon, Q is a positive definite matrix and R is a positive semidefinite matrix which are called weight matrices, w
(
t
) is the reference trajectory, u
(
t
) represents the system inputs, and y
^
(
t
+
j
∣
t
) is a j-step ahead prediction of the system outputs on data up to time t. The predicted outputs can recursively be generated using the process fuzzy model or obtained using a multistep prediction model of the process.Given that the predicted outputs are nonlinear with respect to the control inputs, the optimization problem given by (10) is then a constrained nonlinear and nonconvex optimization problem, the solution of which is difficult and generally expensive in computing time. Different approaches were investigated to solve this problem, such as the numerical optimization techniques [29–31], the metaheuristic based optimization algorithms [7, 15, 22, 32], the linearization of the process fuzzy model [18], and the use of particular model structures to obtain a convex form for the cost function [33].


## 4. PSO Algorithm

The particle swarm optimization is a modern and versatile population based optimization technique, in many respects, similar to evolutionary algorithms [19]. It was first introduced by Kennedy and Eberhart in 1995, as a solution for the optimization problem of a single objective continuous problem. It is based on the observations made upon the collective behaviour of social living entities, such as ants, bees, and flocks of birds. Basically, each particle will be moving around the search space looking for an adequate solution to the optimization problem. By evaluating the fitness of each one of them, their suitability as a possible solution will be assessed. In the next iteration, each particle will change its position to provide more accurate solution, by considering the history of the particle or its best position, its neighborhood history, and mainly the best position found by the whole swarm.The particles movements are bound by the following equations:(11)
v
m
+
1
=
v
m
+
a
m
+
1
,
x
m
+
1
=
x
m
+
v
m
+
1
,where v, a, x are, respectively, the particle, velocity, acceleration, and position and m is the iteration index.In the Clerc-Kennedy PSO [34], considered as the conventional PSO algorithm, the particle acceleration has the following expression:(12)
a
=
χ
c
ϵ
1
P
g
-
x
+
c
ϵ
2
P
-
x
-
1
-
χ
v
,where χ
=
0.729843788, c
=
2.05, ϵ
1 and ϵ
2 are random numbers from the uniform distribution U
[
0,1
], P is the particle personal best position, and P
g is the global best position of all particles. Parameter χ guarantees a decreasing velocity for each particle as the number of iterations increases. This property will improve the convergence quality as the particles approach the solution; their movement will be more and more limited, which help to fine-tune the found solution.The PSO population size is extremely dependent on the optimization problem at hand; a compromise must be reached between the size of the population and the global performance. Empirically, it has been shown that a population of about 20 to 30 particles is generally sufficient in handling almost all the classical optimization problems [35]. Some other optimization problems could be efficiently solved by a smaller population or rather require a larger one. Several papers have dealt with this problem. The author in [35] describes an adaptive PSO algorithm where the swarm population size is obtained through different strategies that regulate the population by removing badly performing particles and creating new ones. In [36] the VarPSO has been proposed, where the PSO population size is varied during the run, based on the concepts of particle age and neighborhood.One of the major drawbacks of the original PSO algorithm is its inability to consistently improve a solution as the number of iterations increases [37]. Several improvements were proposed to address this weakness, such as the tent-map PSO algorithm introduced by Song et al. [38] and the self-organizing hierarchical particle swarm optimization [39]. In [34], a constriction factor χ
χ
<
1 is used to progressively contract the particle movement, as the iterations increase. Different distribution laws like exponential and Gaussian have already been used for the fine-tuning of PSO parameters [40, 41]. But for initializing the particles most of the proposed algorithms use the uniform distribution. The authors in [42] have investigated the possibility of using different distribution law to initialize the particles. They have shown that using other distribution laws gives better performances for most of the cases. In the remaining ones, the performances were the same. This is due to the nature of the uniform distribution, which sometimes gives a wrong impression of the relative performance of algorithms as it is shown in [43].One of the reasons for the poor performance of the basic version of PSO algorithm may be attributed to the dispersion of initial population points in the search space. If the population does not cover the search space efficiently, it may not be able to locate the appropriate solution points, thereby missing the global optimum [44]. This difficulty may be minimized to a great extent by selecting a well-organized distribution law.The PSO algorithm used to solve the optimization problem of the fuzzy predictive control uses the Gaussian distribution to initialize the swarm and is based on the ideas that take advantage of both prior knowledge about the search space landscape and the fact that in most practical applications the dynamic optimization problems changes are gradual [19]. In most cases, it will be logical to assume that once a solution was found, it is better to track it rather than to look for it every time the optimization problem changes. In this algorithm a group of candidate solutions is generated and their fitness is evaluated. The fittest one is adopted by the PSO population as their initial global best position before even initializing the PSO algorithm; then the particles population is used to enhance this solution by favouring the search of an adaptive immediate neighbourhood around the initial solution. This allows increasing the efficiency of the algorithm in regard to the execution time and enhances the convergence quality of the solution. Contrary to the conventional PSO algorithms that involve the use of a large population size to effectively probe the whole search space, the downside, of course, is the huge computing requirement necessary to manage this population.


## 5. Fuzzy Predictive Controller Design Based on the PSO Algorithm

### 5.1. Notation and Principle

Without excluding the use of other nonlinear models, the TS fuzzy model described in [14] is used, as the explicit model to predict the future values of the system. In the case of other nonlinear models, the same steps can be used in solving the optimization problem.Let us define the notation adopted in the formulation of the proposed control algorithm.We consider a MIMO system withm inputs and n outputs.The position of theh
t
h particle, at the sampling time t, is represented by(13)
X
h
t
=
u
h
t
,
u
h
t
+
1
,
…
,
u
h
t
+
N
u
-
1
,
h
=
1
,
…
,
N
pop
,where N
pop is the number of particles in the population (swarm), and (14)
u
h
t
+
j
=
u
1
h
t
+
j
,
…
,
u
m
h
t
+
j
T
,
j
=
0
,
…
,
N
u
-
1
.The position of each particle is adaptively changed according to (11) and (12).The cost functionJ is used as the fitness function (f
fitness
·
=
J
·).The control sequenceU
t which is the optimal solution of the optimization problem is then given by the global best position (X
g
best); that is,(15)
U
t
=
u
t
,
u
t
+
1
,
…
,
u
t
+
N
u
-
1
=
X
g
best
t
,where u
t
+
j
=
u
1
(
t
+
j
)
,
u
2
(
t
+
j
)
,
…
,
u
m
(
t
+
j
)
T, j
=
0
,
…
,
N
u
-
1, is the system control input vector.The block diagram of the developed control algorithm is illustrated in Figure2.Figure 2
Block diagram of the proposed NMPC algorithm.The different steps of the algorithm are drawn based on the following points:(1)
The solution obtained at the sampling timet is taken as an initial solution for the next sampling time t
+
1:(16)
X
1
0
=
U
t
.
This is a good choice when the system is operating in the steady state.
(2)
Another candidate initial solution for the sampling timet
+
1 is given by(17)
X
2
0
=
u
t
+
1
,
u
t
+
2
,
…
,
u
t
+
N
u
,
whereu
t
+
1
,
…
,
u
t
+
N
u
-
1 is the control sequence obtained at the previous sampling time and u
t
+
N
u
=
u
t
+
N
u
-
1.
(3)
Both initial solutions chosen in points 1 and 2 are evaluated using the fitness function. The fittest one is chosen as the initial global best positionX
best
0
t
=
u
0
t
,
u
0
t
+
1
,
…
,
u
0
t
+
N
u
-
1 of the entire PSO population.
(4)
Using the Gaussian distribution, the particles (X
h, h
=
1
,
…
,
N
pop) are distributed around X
g
best
0 within the radius r
d
=
r
d
1
⋯
r
d
m
T according to the following expression:(18)
u
i
h
t
+
j
=
u
i
0
t
+
j
+
r
d
i
y
j
k
m
N
u
dist
⁡
y
j
=
0
,
…
,
N
u
-
1
;
i
=
1
,
…
,
m
,
wherey
j is random value from the normal distribution N
[
0,1
], k
m is a random value from the normal distribution N
[
0
,
1
/
3
], and dist
⁡
y
=
∑
j
=
1
N
u
y
j
2.
The values of the radiusr
d, at the sampling time t
+
1, are given by(19)
r
d
=
u
t
+
1
-
u
t
.To explore the whole search space, the particles positions are updated according to (11) and (12).Using a Gaussian distribution will ensure that the density of the particles gets higher around the chosen initial global best position; thus the optimal solution can be reached in few iterations, using a small number of particles. This is very interesting, in the case of small changes of the system outputs, if the system is operating in the steady state, or if we want to fine-tune the solution. The use of a uniform distribution favors the diversity of the particles and increases the probability to find better solution relatively far from the previous one. Since the number of PSO particles and the allowed iterations number should not be large, their initial positions carry a lot of influence in the quality of the solution.To keep the algorithm from converging toward a local optimum, a minimal valuer
d
min
⁡
=
r
d
1
m
i
n
⋯
r
d
m
min
⁡
T given by the following expression is imposed to r
d:(20)
r
d
min
⁡
=
α
1
w
t
-
y
t
+
α
2
,where α
1 is an n
×
m matrix of scaling parameters and α
2 (m
×
1) is used to impose a minimum value to r
d
min
⁡ regardless of the tracking error (the value of r
d
min
⁡ can have a fixed value by setting α
1 to zero).If the current solution is far from the global optimum, the corresponding tracking error will be large and consequently the values ofr
d
min
⁡ will also be large. When r
d
i
min
⁡
>
r
d
i   (i
=
1
,
…
,
m), the particles are distributed around X
g
best
0 within r
d
i
min
⁡ instead of r
d
i. This will allow the algorithm to escape the current local optimum by looking for a better solution far from the actual one.(5)
The termination criterion is the maximal number of iterations.


### 5.2. Constraints Handling

One of the strong points of predictive control is its ability to handle the control constraints directly and efficiently in the design stage by solving a constrained optimization problem. Different types of constraints could be distinguished depending on their nature. Mainly, we have constraints on the inputs, presented by actuators with limited range of action and a limited slew rate, and constraints on the outputs due to environmental, safety, and economic goals.Several approaches were proposed for handling constraints using evolutionary algorithms [45, 46]. They can be grouped into four categories [47]: methods based on preserving feasibility of solutions, methods based on penalty functions, methods that differentiate the feasible and infeasible solutions, and hybrid methods. The proposed control algorithm is based, in the constraints handling, on the hybrid approach. Indeed, the output constraints are handled using a penalty function while the input constraints are handled using the method based on preserving feasibility of solutions.


#### 5.2.1. Outputs Constraints

One practical and popular approach is to soften the constraints [48]; this means that a violation of the constraints is allowed, but under severe and special conditions. The easy way to soften the outputs constraints is to add new variables, called slack variables, to the cost function that will heavily penalize any deviation or violation of the constraints in order to force the optimizer to violate these constraints only when it is absolutely necessary. The output-dependent weight function Γ
y
y was chosen to soften the output constraints. It has the following expression [5]:(21)
Γ
y
y
=
Γ
y
1
y
1
0
⋯
0
0
Γ
y
2
y
2
⋯
0
0
⋮
⋱
0
0
⋯
0
Γ
y
n
y
n
,
Γ
y
i
y
i
=
Γ
y
i
0
1
+
ε
i
y
i
-
y
i
min
⁡
2
if
y
i
≤
y
i
min
⁡
Γ
y
i
0
if
y
i
min
⁡
≤
y
i
≤
y
i
max
⁡
Γ
y
i
0
1
+
ε
i
y
i
-
y
i
max
⁡
2
if
y
i
≥
y
i
max
⁡
,where i
=
1
,
…
,
n and y
t
=
y
1
(
t
)
,
y
2
(
t
)
,
…
,
y
n
(
t
)
T
.The function depicted in Figure3 is the output-dependent weight function (Γ
y
i) used in the softening of the constraints.Figure 3
Weight functionΓ
y
i
y.ε
i is used to define the degree of softening: ε
i
=
∞ indicates hard constraint while ε
i
=
0 indicates no constraint. y
m
i
n
=
y
1
m
i
n
,
…
,
y
n
m
i
n
T and y
m
a
x
=
y
1
m
a
x
,
…
,
y
n
m
a
x
T limit the range where the outputs are allowed to exist.The functionΓ
y
i will be also used on unconstrained outputs (y
i) by setting y
i
m
i
n and y
i
m
a
x, respectively, to -
∞ and +
∞.By introducing the weight functionΓ
y
y, the optimization problem will be reformulated to have the following expression:(22)
min
⁡
Δ
u
t
J
u
t
,
y
^
t
,
w
t
=
∑
j
=
N
1
N
2
y
^
t
+
j
∣
t
-
w
t
+
j
T
Γ
y
y
^
t
+
j
∣
t
-
w
t
+
j
+
∑
j
=
1
N
u
Δ
u
T
t
+
j
-
1
R
Δ
u
t
+
j
-
1
subject to:
Δ
u
t
+
j
-
1
=
0
,
for
j
>
N
u
u
min
⁡
<
u
t
+
j
<
u
max
⁡
for
1
≤
j
≤
N
u
Δ
u
min
⁡
<
Δ
u
t
+
j
-
1
<
Δ
u
max
⁡
for
j
≤
N
u
.


#### 5.2.2. Inputs Constraints

The PSO algorithm, as any other population based algorithms, can systematically handle the input constraints by searching for solutions in the space for which these constraints are satisfied. The lower and upper bounds of this space are obtained as follows.The constraints on the inputs and their increments are typically given by(23)
u
min
⁡
<
u
t
<
u
max
⁡
(24)
Δ
u
min
⁡
<
Δ
u
t
<
Δ
u
max
⁡
,where u
m
i
n
=
u
1
m
i
n
,
…
,
u
m
m
i
n
T and u
m
a
x
=
u
1
m
a
x
,
…
,
u
m
m
a
x
T limit the range of the inputs, while Δ
u
m
i
n
=
Δ
u
1
m
i
n
,
…
,
Δ
u
m
m
i
n
T and Δ
u
m
a
x
=
Δ
u
1
m
a
x
,
…
,
Δ
u
m
m
a
x
T limit the range of the inputs increments.For any given particleh
(
h
=
1
,
…
,
N
pop
), we can write, from (24),(25)
Δ
u
i
min
⁡
≤
Δ
u
i
h
t
+
j
∣
t
≤
Δ
u
i
max
⁡
⟹
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
≤
u
i
h
t
+
j
∣
t
=
u
i
h
t
+
j
-
1
∣
t
+
Δ
u
i
h
t
+
j
∣
t
≤
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
.For j
=
0, u
i
h
t
-
1
∣
t is the already calculated control action u
i
(
t
-
1
), and Δ
u
i
h
t
+
j
∣
t
=
u
i
h
t
+
j
∣
t
-
u
i
h
t
+
j
-
1
∣
t.On the other hand, from (23), the constraints on the inputs magnitude could also be written as follows:(26)
u
i
min
⁡
≤
u
i
h
t
+
j
∣
t
≤
u
i
max
⁡
,
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
.Hence, from (25) and (26), the lower and upper bounds of the search space will be given by(27)
L
i
min
⁡
h
t
+
j
≤
u
i
h
t
+
j
∣
t
≤
L
i
max
⁡
h
t
+
j
,
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
,where(28)
L
i
min
⁡
h
t
+
j
=
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
if
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
>
u
i
min
⁡
u
i
min
⁡
Otherwise,
L
i
max
⁡
h
t
+
j
=
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
if
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
<
u
i
max
⁡
u
i
max
⁡
Otherwise.


### 5.3. The Control Algorithm

The basic steps of the proposed control algorithm are summarized as follows.Step 1.
Initialize the control algorithm by choosing suitable values for the design and the PSO algorithm parameters.Step 2.
According to points 1 and 2, determine the initial solutionsX
1
0 and X
2
0.Step 3.
Specify the desired reference trajectory over the prediction horizon and compute the future values of the system outputs:y
^
t
+
j, N
1
≤
j
≤
N
2.Step 4.
Use the cost function(
J
) given by (22) to evaluate the fitness of the initial solutions X
1
0 and X
2
0. The fittest one will be adopted by the PSO population as their initial global best position X
g
best
0.
SetX
g
best
=
X
g
best
0.Step 5.
Evaluate the value ofr
d using (19) and the value of r
d
min
⁡ using (20):(29)
If
r
d
i
<
r
d
i
min
⁡
then
r
d
i
=
r
d
i
min
⁡
,
i
=
1
,
…
,
m
.Step 6.
For each particleh,randomly initialize its velocityv
h,
using (18), distribute the particle X
h around X
g
best within the radius r
d,
setP
h
=
X
h, and evaluate f
fitness
P
h (P
h is the personal best position of particle X
h)(30)
If
f
fitness
P
h
<
f
fitness
X
g
best
then
X
g
best
=
P
h
.Step 7 (repeat).
For each particleh,using (11) and (12), update X
h,
evaluatef
fitness
(
X
h
),(31)
If
f
fitness
X
h
<
f
fitness
P
h
P
h
=
X
h
if
f
fitness
P
h
<
f
fitness
X
g
best
then
X
g
best
=
P
h
,until termination criterion is reached.Step 8.
The optimal solutionU
t
=
u
t
,
u
t
+
1
⋯
u
t
+
N
u
-
1 is the global best position X
g
best.For the next sampling time go to Step2.


## 5.1. Notation and Principle

Without excluding the use of other nonlinear models, the TS fuzzy model described in [14] is used, as the explicit model to predict the future values of the system. In the case of other nonlinear models, the same steps can be used in solving the optimization problem.Let us define the notation adopted in the formulation of the proposed control algorithm.We consider a MIMO system withm inputs and n outputs.The position of theh
t
h particle, at the sampling time t, is represented by(13)
X
h
t
=
u
h
t
,
u
h
t
+
1
,
…
,
u
h
t
+
N
u
-
1
,
h
=
1
,
…
,
N
pop
,where N
pop is the number of particles in the population (swarm), and (14)
u
h
t
+
j
=
u
1
h
t
+
j
,
…
,
u
m
h
t
+
j
T
,
j
=
0
,
…
,
N
u
-
1
.The position of each particle is adaptively changed according to (11) and (12).The cost functionJ is used as the fitness function (f
fitness
·
=
J
·).The control sequenceU
t which is the optimal solution of the optimization problem is then given by the global best position (X
g
best); that is,(15)
U
t
=
u
t
,
u
t
+
1
,
…
,
u
t
+
N
u
-
1
=
X
g
best
t
,where u
t
+
j
=
u
1
(
t
+
j
)
,
u
2
(
t
+
j
)
,
…
,
u
m
(
t
+
j
)
T, j
=
0
,
…
,
N
u
-
1, is the system control input vector.The block diagram of the developed control algorithm is illustrated in Figure2.Figure 2
Block diagram of the proposed NMPC algorithm.The different steps of the algorithm are drawn based on the following points:(1)
The solution obtained at the sampling timet is taken as an initial solution for the next sampling time t
+
1:(16)
X
1
0
=
U
t
.
This is a good choice when the system is operating in the steady state.
(2)
Another candidate initial solution for the sampling timet
+
1 is given by(17)
X
2
0
=
u
t
+
1
,
u
t
+
2
,
…
,
u
t
+
N
u
,
whereu
t
+
1
,
…
,
u
t
+
N
u
-
1 is the control sequence obtained at the previous sampling time and u
t
+
N
u
=
u
t
+
N
u
-
1.
(3)
Both initial solutions chosen in points 1 and 2 are evaluated using the fitness function. The fittest one is chosen as the initial global best positionX
best
0
t
=
u
0
t
,
u
0
t
+
1
,
…
,
u
0
t
+
N
u
-
1 of the entire PSO population.
(4)
Using the Gaussian distribution, the particles (X
h, h
=
1
,
…
,
N
pop) are distributed around X
g
best
0 within the radius r
d
=
r
d
1
⋯
r
d
m
T according to the following expression:(18)
u
i
h
t
+
j
=
u
i
0
t
+
j
+
r
d
i
y
j
k
m
N
u
dist
⁡
y
j
=
0
,
…
,
N
u
-
1
;
i
=
1
,
…
,
m
,
wherey
j is random value from the normal distribution N
[
0,1
], k
m is a random value from the normal distribution N
[
0
,
1
/
3
], and dist
⁡
y
=
∑
j
=
1
N
u
y
j
2.
The values of the radiusr
d, at the sampling time t
+
1, are given by(19)
r
d
=
u
t
+
1
-
u
t
.To explore the whole search space, the particles positions are updated according to (11) and (12).Using a Gaussian distribution will ensure that the density of the particles gets higher around the chosen initial global best position; thus the optimal solution can be reached in few iterations, using a small number of particles. This is very interesting, in the case of small changes of the system outputs, if the system is operating in the steady state, or if we want to fine-tune the solution. The use of a uniform distribution favors the diversity of the particles and increases the probability to find better solution relatively far from the previous one. Since the number of PSO particles and the allowed iterations number should not be large, their initial positions carry a lot of influence in the quality of the solution.To keep the algorithm from converging toward a local optimum, a minimal valuer
d
min
⁡
=
r
d
1
m
i
n
⋯
r
d
m
min
⁡
T given by the following expression is imposed to r
d:(20)
r
d
min
⁡
=
α
1
w
t
-
y
t
+
α
2
,where α
1 is an n
×
m matrix of scaling parameters and α
2 (m
×
1) is used to impose a minimum value to r
d
min
⁡ regardless of the tracking error (the value of r
d
min
⁡ can have a fixed value by setting α
1 to zero).If the current solution is far from the global optimum, the corresponding tracking error will be large and consequently the values ofr
d
min
⁡ will also be large. When r
d
i
min
⁡
>
r
d
i   (i
=
1
,
…
,
m), the particles are distributed around X
g
best
0 within r
d
i
min
⁡ instead of r
d
i. This will allow the algorithm to escape the current local optimum by looking for a better solution far from the actual one.(5)
The termination criterion is the maximal number of iterations.


## 5.2. Constraints Handling

One of the strong points of predictive control is its ability to handle the control constraints directly and efficiently in the design stage by solving a constrained optimization problem. Different types of constraints could be distinguished depending on their nature. Mainly, we have constraints on the inputs, presented by actuators with limited range of action and a limited slew rate, and constraints on the outputs due to environmental, safety, and economic goals.Several approaches were proposed for handling constraints using evolutionary algorithms [45, 46]. They can be grouped into four categories [47]: methods based on preserving feasibility of solutions, methods based on penalty functions, methods that differentiate the feasible and infeasible solutions, and hybrid methods. The proposed control algorithm is based, in the constraints handling, on the hybrid approach. Indeed, the output constraints are handled using a penalty function while the input constraints are handled using the method based on preserving feasibility of solutions.


### 5.2.1. Outputs Constraints

One practical and popular approach is to soften the constraints [48]; this means that a violation of the constraints is allowed, but under severe and special conditions. The easy way to soften the outputs constraints is to add new variables, called slack variables, to the cost function that will heavily penalize any deviation or violation of the constraints in order to force the optimizer to violate these constraints only when it is absolutely necessary. The output-dependent weight function Γ
y
y was chosen to soften the output constraints. It has the following expression [5]:(21)
Γ
y
y
=
Γ
y
1
y
1
0
⋯
0
0
Γ
y
2
y
2
⋯
0
0
⋮
⋱
0
0
⋯
0
Γ
y
n
y
n
,
Γ
y
i
y
i
=
Γ
y
i
0
1
+
ε
i
y
i
-
y
i
min
⁡
2
if
y
i
≤
y
i
min
⁡
Γ
y
i
0
if
y
i
min
⁡
≤
y
i
≤
y
i
max
⁡
Γ
y
i
0
1
+
ε
i
y
i
-
y
i
max
⁡
2
if
y
i
≥
y
i
max
⁡
,where i
=
1
,
…
,
n and y
t
=
y
1
(
t
)
,
y
2
(
t
)
,
…
,
y
n
(
t
)
T
.The function depicted in Figure3 is the output-dependent weight function (Γ
y
i) used in the softening of the constraints.Figure 3
Weight functionΓ
y
i
y.ε
i is used to define the degree of softening: ε
i
=
∞ indicates hard constraint while ε
i
=
0 indicates no constraint. y
m
i
n
=
y
1
m
i
n
,
…
,
y
n
m
i
n
T and y
m
a
x
=
y
1
m
a
x
,
…
,
y
n
m
a
x
T limit the range where the outputs are allowed to exist.The functionΓ
y
i will be also used on unconstrained outputs (y
i) by setting y
i
m
i
n and y
i
m
a
x, respectively, to -
∞ and +
∞.By introducing the weight functionΓ
y
y, the optimization problem will be reformulated to have the following expression:(22)
min
⁡
Δ
u
t
J
u
t
,
y
^
t
,
w
t
=
∑
j
=
N
1
N
2
y
^
t
+
j
∣
t
-
w
t
+
j
T
Γ
y
y
^
t
+
j
∣
t
-
w
t
+
j
+
∑
j
=
1
N
u
Δ
u
T
t
+
j
-
1
R
Δ
u
t
+
j
-
1
subject to:
Δ
u
t
+
j
-
1
=
0
,
for
j
>
N
u
u
min
⁡
<
u
t
+
j
<
u
max
⁡
for
1
≤
j
≤
N
u
Δ
u
min
⁡
<
Δ
u
t
+
j
-
1
<
Δ
u
max
⁡
for
j
≤
N
u
.


### 5.2.2. Inputs Constraints

The PSO algorithm, as any other population based algorithms, can systematically handle the input constraints by searching for solutions in the space for which these constraints are satisfied. The lower and upper bounds of this space are obtained as follows.The constraints on the inputs and their increments are typically given by(23)
u
min
⁡
<
u
t
<
u
max
⁡
(24)
Δ
u
min
⁡
<
Δ
u
t
<
Δ
u
max
⁡
,where u
m
i
n
=
u
1
m
i
n
,
…
,
u
m
m
i
n
T and u
m
a
x
=
u
1
m
a
x
,
…
,
u
m
m
a
x
T limit the range of the inputs, while Δ
u
m
i
n
=
Δ
u
1
m
i
n
,
…
,
Δ
u
m
m
i
n
T and Δ
u
m
a
x
=
Δ
u
1
m
a
x
,
…
,
Δ
u
m
m
a
x
T limit the range of the inputs increments.For any given particleh
(
h
=
1
,
…
,
N
pop
), we can write, from (24),(25)
Δ
u
i
min
⁡
≤
Δ
u
i
h
t
+
j
∣
t
≤
Δ
u
i
max
⁡
⟹
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
≤
u
i
h
t
+
j
∣
t
=
u
i
h
t
+
j
-
1
∣
t
+
Δ
u
i
h
t
+
j
∣
t
≤
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
.For j
=
0, u
i
h
t
-
1
∣
t is the already calculated control action u
i
(
t
-
1
), and Δ
u
i
h
t
+
j
∣
t
=
u
i
h
t
+
j
∣
t
-
u
i
h
t
+
j
-
1
∣
t.On the other hand, from (23), the constraints on the inputs magnitude could also be written as follows:(26)
u
i
min
⁡
≤
u
i
h
t
+
j
∣
t
≤
u
i
max
⁡
,
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
.Hence, from (25) and (26), the lower and upper bounds of the search space will be given by(27)
L
i
min
⁡
h
t
+
j
≤
u
i
h
t
+
j
∣
t
≤
L
i
max
⁡
h
t
+
j
,
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
,where(28)
L
i
min
⁡
h
t
+
j
=
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
if
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
>
u
i
min
⁡
u
i
min
⁡
Otherwise,
L
i
max
⁡
h
t
+
j
=
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
if
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
<
u
i
max
⁡
u
i
max
⁡
Otherwise.


## 5.2.1. Outputs Constraints

One practical and popular approach is to soften the constraints [48]; this means that a violation of the constraints is allowed, but under severe and special conditions. The easy way to soften the outputs constraints is to add new variables, called slack variables, to the cost function that will heavily penalize any deviation or violation of the constraints in order to force the optimizer to violate these constraints only when it is absolutely necessary. The output-dependent weight function Γ
y
y was chosen to soften the output constraints. It has the following expression [5]:(21)
Γ
y
y
=
Γ
y
1
y
1
0
⋯
0
0
Γ
y
2
y
2
⋯
0
0
⋮
⋱
0
0
⋯
0
Γ
y
n
y
n
,
Γ
y
i
y
i
=
Γ
y
i
0
1
+
ε
i
y
i
-
y
i
min
⁡
2
if
y
i
≤
y
i
min
⁡
Γ
y
i
0
if
y
i
min
⁡
≤
y
i
≤
y
i
max
⁡
Γ
y
i
0
1
+
ε
i
y
i
-
y
i
max
⁡
2
if
y
i
≥
y
i
max
⁡
,where i
=
1
,
…
,
n and y
t
=
y
1
(
t
)
,
y
2
(
t
)
,
…
,
y
n
(
t
)
T
.The function depicted in Figure3 is the output-dependent weight function (Γ
y
i) used in the softening of the constraints.Figure 3
Weight functionΓ
y
i
y.ε
i is used to define the degree of softening: ε
i
=
∞ indicates hard constraint while ε
i
=
0 indicates no constraint. y
m
i
n
=
y
1
m
i
n
,
…
,
y
n
m
i
n
T and y
m
a
x
=
y
1
m
a
x
,
…
,
y
n
m
a
x
T limit the range where the outputs are allowed to exist.The functionΓ
y
i will be also used on unconstrained outputs (y
i) by setting y
i
m
i
n and y
i
m
a
x, respectively, to -
∞ and +
∞.By introducing the weight functionΓ
y
y, the optimization problem will be reformulated to have the following expression:(22)
min
⁡
Δ
u
t
J
u
t
,
y
^
t
,
w
t
=
∑
j
=
N
1
N
2
y
^
t
+
j
∣
t
-
w
t
+
j
T
Γ
y
y
^
t
+
j
∣
t
-
w
t
+
j
+
∑
j
=
1
N
u
Δ
u
T
t
+
j
-
1
R
Δ
u
t
+
j
-
1
subject to:
Δ
u
t
+
j
-
1
=
0
,
for
j
>
N
u
u
min
⁡
<
u
t
+
j
<
u
max
⁡
for
1
≤
j
≤
N
u
Δ
u
min
⁡
<
Δ
u
t
+
j
-
1
<
Δ
u
max
⁡
for
j
≤
N
u
.


## 5.2.2. Inputs Constraints

The PSO algorithm, as any other population based algorithms, can systematically handle the input constraints by searching for solutions in the space for which these constraints are satisfied. The lower and upper bounds of this space are obtained as follows.The constraints on the inputs and their increments are typically given by(23)
u
min
⁡
<
u
t
<
u
max
⁡
(24)
Δ
u
min
⁡
<
Δ
u
t
<
Δ
u
max
⁡
,where u
m
i
n
=
u
1
m
i
n
,
…
,
u
m
m
i
n
T and u
m
a
x
=
u
1
m
a
x
,
…
,
u
m
m
a
x
T limit the range of the inputs, while Δ
u
m
i
n
=
Δ
u
1
m
i
n
,
…
,
Δ
u
m
m
i
n
T and Δ
u
m
a
x
=
Δ
u
1
m
a
x
,
…
,
Δ
u
m
m
a
x
T limit the range of the inputs increments.For any given particleh
(
h
=
1
,
…
,
N
pop
), we can write, from (24),(25)
Δ
u
i
min
⁡
≤
Δ
u
i
h
t
+
j
∣
t
≤
Δ
u
i
max
⁡
⟹
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
≤
u
i
h
t
+
j
∣
t
=
u
i
h
t
+
j
-
1
∣
t
+
Δ
u
i
h
t
+
j
∣
t
≤
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
.For j
=
0, u
i
h
t
-
1
∣
t is the already calculated control action u
i
(
t
-
1
), and Δ
u
i
h
t
+
j
∣
t
=
u
i
h
t
+
j
∣
t
-
u
i
h
t
+
j
-
1
∣
t.On the other hand, from (23), the constraints on the inputs magnitude could also be written as follows:(26)
u
i
min
⁡
≤
u
i
h
t
+
j
∣
t
≤
u
i
max
⁡
,
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
.Hence, from (25) and (26), the lower and upper bounds of the search space will be given by(27)
L
i
min
⁡
h
t
+
j
≤
u
i
h
t
+
j
∣
t
≤
L
i
max
⁡
h
t
+
j
,
i
=
1
,
…
,
m
;
j
=
0
,
…
,
N
u
-
1
,where(28)
L
i
min
⁡
h
t
+
j
=
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
if
Δ
u
i
min
⁡
+
u
i
h
t
+
j
-
1
∣
t
>
u
i
min
⁡
u
i
min
⁡
Otherwise,
L
i
max
⁡
h
t
+
j
=
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
if
Δ
u
i
max
⁡
+
u
i
h
t
+
j
-
1
∣
t
<
u
i
max
⁡
u
i
max
⁡
Otherwise.


## 5.3. The Control Algorithm

The basic steps of the proposed control algorithm are summarized as follows.Step 1.
Initialize the control algorithm by choosing suitable values for the design and the PSO algorithm parameters.Step 2.
According to points 1 and 2, determine the initial solutionsX
1
0 and X
2
0.Step 3.
Specify the desired reference trajectory over the prediction horizon and compute the future values of the system outputs:y
^
t
+
j, N
1
≤
j
≤
N
2.Step 4.
Use the cost function(
J
) given by (22) to evaluate the fitness of the initial solutions X
1
0 and X
2
0. The fittest one will be adopted by the PSO population as their initial global best position X
g
best
0.
SetX
g
best
=
X
g
best
0.Step 5.
Evaluate the value ofr
d using (19) and the value of r
d
min
⁡ using (20):(29)
If
r
d
i
<
r
d
i
min
⁡
then
r
d
i
=
r
d
i
min
⁡
,
i
=
1
,
…
,
m
.Step 6.
For each particleh,randomly initialize its velocityv
h,
using (18), distribute the particle X
h around X
g
best within the radius r
d,
setP
h
=
X
h, and evaluate f
fitness
P
h (P
h is the personal best position of particle X
h)(30)
If
f
fitness
P
h
<
f
fitness
X
g
best
then
X
g
best
=
P
h
.Step 7 (repeat).
For each particleh,using (11) and (12), update X
h,
evaluatef
fitness
(
X
h
),(31)
If
f
fitness
X
h
<
f
fitness
P
h
P
h
=
X
h
if
f
fitness
P
h
<
f
fitness
X
g
best
then
X
g
best
=
P
h
,until termination criterion is reached.Step 8.
The optimal solutionU
t
=
u
t
,
u
t
+
1
⋯
u
t
+
N
u
-
1 is the global best position X
g
best.For the next sampling time go to Step2.


## 6. Case Study

### 6.1. Process Description

To assess the performance of the proposed control algorithm, the control of the industrial boiler model developed by Pellegrinetti and Bentsman [28] is considered. This model represents a dual fuel (oil and gas) boiler in the Abbott Power Plant in Champaign, Illinois, which could be used for both heating and electric generation. It is a MIMO process, which has four manipulated inputs (fuel flow rate, air flow rate, feed water flow rate, and a valve to regulate the steam demand) and four measured outputs (pressure, oxygen level in the flue gas, drum water level, and steam flow). The control objectives are to maintain the pressure, the level of the water in the drum, and the oxygen level in the flue gas at their desired values, regardless of the steam flow rate provided by the boiler and any other disturbances, such as the fluctuations in the heating value of the fuel or the variations in the ambient temperature. This mathematical model includes perturbations and measurement noises and is described by the following equations:(32)
x
˙
1
t
=
c
11
x
4
t
x
1
9
/
8
t
+
c
12
u
1
t
-
τ
1
-
c
13
u
3
t
-
τ
3
+
c
14
,
x
˙
2
t
=
-
c
21
x
2
t
+
c
22
u
2
t
-
τ
2
-
c
23
u
1
t
-
τ
1
-
c
24
u
1
t
-
τ
1
x
2
t
c
25
u
2
t
-
τ
2
+
c
26
u
1
t
-
τ
1
,
x
˙
3
t
=
c
31
x
1
t
-
c
32
x
4
t
x
1
t
+
c
33
u
3
t
-
τ
3
,
x
˙
4
t
=
-
c
41
x
4
t
+
c
42
u
1
t
-
τ
1
+
c
43
+
c
41
u
4
t
,
y
1
t
=
c
51
x
1
t
-
τ
4
+
n
1
t
,
y
2
t
=
c
61
x
2
t
-
τ
5
+
n
2
t
,
y
3
t
=
c
70
x
1
t
-
τ
6
+
c
71
x
3
t
-
τ
6
+
c
72
x
4
t
-
τ
6
x
1
t
-
τ
6
+
c
73
u
3
t
-
τ
3
-
τ
6
+
c
74
u
1
t
-
τ
1
-
τ
6
+
c
75
x
1
t
-
τ
6
+
c
76
1
-
c
77
x
3
t
-
τ
6
x
3
t
-
τ
6
x
1
t
-
τ
6
+
c
78
+
c
79
+
n
3
t
,
y
4
t
=
c
81
x
4
t
-
τ
7
+
c
82
x
1
t
-
τ
7
+
n
4
t
,where x
˙
1 is the drum pressure state K
g
f
/
c
m
2, y
1 is the measured drum pressure (PSI), y
2 and x
˙
2 are the measured excess oxygen level and its state, respectively (%), x
˙
3 is the system fluid density K
g
/
m
3, y
3 is the drum water level (in), y
4 is the steam flow rate K
g
/
s, u
1, u
2, and u
3 are, respectively, the fuel, air, and feed water flow rates which take values in the interval 0
1, and x
˙
4 is the exogenous variable related to the steam demanded. The variables n
i are the outputs of the first-order colored noise models driven by zero mean and unit variance white noise and are considered to be unmeasured output disturbances. The remaining coefficients values are given in Table 1. Consider(33)
n
1
=
0.75
s
+
0.1
s
+
0.001
w
1
,
n
2
=
0.019
s
+
0.001
s
+
0.024
w
2
,
n
3
=
0.105
s
+
0.038
s
+
0.01
w
3
,
n
4
=
0.01
s
+
0.0001
s
+
0.001
w
4
,where n
i, i
=
1
,
…
,
4, are colored noise and w
i, i
=
1
,
…
,
4, is the unit variance white noise.Table 1
Model equations coefficients.
c
11
=
-
0.00478
c
12
=
0.28
c
13
=
0.01348
c
14
=
0.02493
c
21
=
0.1540357
c
22
=
103.5462
c
23
=
107.4835
c
24
=
1.9515
c
25
=
29.04
c
26
=
1.824
c
31
=
0.00533176
c
32
=
0.025195
c
33
=
0.7317058
c
41
=
0.04
c
42
=
0.0299886
c
43
=
0.018088
c
51
=
14.214
c
61
=
1.00
c
81
=
0.85663
c
82
=
-
0.18128
c
70
=
-
0.1048569
c
71
=
0.15479
c
72
=
0.4954961
c
73
=
-
0.20797
c
74
=
1.2720
c
75
=
-
324212.7805
c
76
=
-
99556.24778
c
77
=
0.001185
c
78
=
-
1704.50476
c
79
=
-
103.7351
τ
1
=
2
τ
2
=
2
τ
3
=
3
τ
4
=
3
τ
5
=
4
τ
6
=
10
τ
7
=
2The process is supposed to be initially operating around the following states:(34)
x
0
=
22.5,2.5,621.17,0.6941
T
,
y
0
=
320,2.5,0
,
9.3053
T
,
u
0
=
0.3227,0.39503,0.37404,0
T
.The linear representation of the boiler around the operating state given by (34) is described by(35)
x
˙
=
A
b
x
+
B
b
u
,
y
=
C
b
x
+
D
b
u
,where(36)
A
b
=
-
0.005508
0
0
-
0.15872
0
-
0.20625
0
0
-
0.012156
0
0
-
0.566887
0
0
0
-
0.04
,
B
b
=
0.28
0
-
0.01348
0
-
9.374893
7.658411
0
0
0
0
0.731706
0
0.029989
0
0
0.04
,
C
b
=
14.214
0
0
0
0
1
0
0
0.322072
0
0.149410
11.148662
0.413307
0
0
19.274175
,
D
b
=
0
0
0
0
0
0
0
0
1.272
0
-
0.20797
0
0
0
0
0
.There is a direct feedthrough between the outputy
3 and the inputs u
1 and u
3, expressed by a nonzero D matrix in the state representation. This property will create some problems when considering the linear MPC controller [49]. In order to have D
b = 0, a delay will be introduced to the inputs u
1 and u
3 with regard to the output y
3.This process is inherently unstable. Essentially, a stabilisation scheme of the water level must be introduced to make any identification approaches possible by incorporating a proportional feedforward action (C2) of 0.0403 and a PI controller (C1) withK
p
=
0.258 and T
i
=
1.1026
e
-
4, as it is illustrated in Figure 4.Figure 4
Stabilization scheme for the boiler.


### 6.2. Process Fuzzy Identification

The first step in the design of a fuzzy NMPC controller is the construction of a nonlinear model of the process. Four fuzzy multi-inputs single output models were constructed to predict the future behaviour of the process. Two triangular membership functions per input, a sampling period of 3 s, and the models given by the following equations were used:(37)
y
^
1
=
f
1
y
^
1
k
-
1
,
u
1
k
-
5
,
x
^
4
k
-
3
,
y
^
3
=
f
2
y
^
3
k
-
1
,
u
1
k
-
12
,
r
Level
k
-
13
,
u
4
k
-
10
,
x
4
k
-
10
,
y
^
4
=
f
3
y
^
4
k
-
1
,
u
1
k
-
4
,
u
4
k
-
2
,
x
4
k
-
2
,
x
^
4
=
f
4
x
^
4
k
-
1
,
u
1
k
-
2
,
u
4
k
-
1
,where f
i, i
=
1
,
…
,
4, represent the fuzzy relationships between the inputs and the outputs.Figure5 shows the validation results of these models.Figure 5
Validations of the fuzzy models.


### 6.3. Controller Implementation and Simulation Results

It has been shown in [28] that there is a nearly linear relationship between the fuel/air ratio (FAR
=
u
2
/
u
1) and the oxygen level. So, the oxygen level y
2 can be kept constant at a desired value by maintaining the corresponding FAR value constant. Given that u
1 affects all the systems outputs, varying u
2 to maintain the FAR value constant is the appropriate way. Thus, the value of u
2 can be directly deduced from that of u
1. To keep the desired oxygen level in flue gas y
2 around 2.5% the FAR value must be equal to 1.2241 and u
1 equals 0.8169 u
2.The control objective is to force the steam flow ratey
4 to track a desired reference trajectory, while fulfilling the following constraints on the manipulated inputs,(38)
Δ
u
1
Max
⁡
=
0.1
,
Δ
u
2
Max
⁡
=
0.1
,
Δ
u
r
Level
Max
⁡
=
1
-
4
<
r
Level
<
4
,
Δ
u
4
Max
⁡
=
0.1
-
0.2
<
u
4
<
0.35
,
0
<
u
3
<
1
,and the following constraints on the outputs,(39)
y
1
=
320
,
y
2
=
2.5
,
y
3
=
0
.The following parameters values were used in implementing the proposed control algorithm:(40)
R
=
100
0
0
0
0
0
0
0
0
0
2
0
0
0
0
130
,
α
1
=
1.65
e
-
4
0
9.54
e
-
5
3
e
-
3
0
0
0
0
1.65
e
-
3
0
1.05
e
-
2
2.76
e
-
2
1.28
e
-
4
0
1.246
e
-
4
3
e
-
3
.And N
u
=
2, number of particles = 5, number of iterations = 5, and α
2
=
0
0
0
0
T. The values of the other parameters are gathered in Table 2.Table 2
Design parameters.
N
1
N
2
ε
i
Γ
0
y
min
⁡
y
max
⁡
y
^
1
1
10
1
0.0085
318
322
y
^
3
10
20
0.7
0.0065
−1.5
1.5
y
^
4
1
10
0
1.006
/
/
x
^
4
1
10
0
0
/
/In addition to the proposed algorithm, two other control strategies were considered:Using the Conventional PSO. In this algorithm, the particles are randomly distributed in the whole search space and the solution of the previous sampling time is taken as a possible solution for the current sampling time. The PSO parameters values χ
=
0.729843788 and c
=
2.05, adopted by Clerc and Kennedy [34], are used in both PSO based controllers.Using the Linear MPC Strategy. In this strategy, only the constraints on the inputs variables and their increments are considered and the optimization problem is converted into a quadratic programming (QP) problem.To implement this strategy, a linear representation for the process given in Figure4 must be used. To remove the direct feedthrough between the output y
3 and inputs u
1 and u
3, the structure given in Figure 4 is slightly modified by introducing a unit delay in these inputs (Figure 6). The obtained system is linearized around the operating state given by (34) and then discretized using a sampling period of 3 s. The following discrete state space representation is obtained:(41)
x
˙
k
+
1
=
A
x
k
+
B
u
k
,
y
k
+
1
=
C
x
k
,where(42)
A
=
0.9861
0
0.001485
-
0.3706
0
0.5386
0
0
-
0.1736
0
0.9187
-
5.665
0
0
0
0.8869
,
B
=
0.8172
0
-
0.009937
-
0.02264
-
20.97
17.13
0
0
-
0.3384
0
0.5441
-
0.3525
0.08478
0
0
0.1131
,
C
=
14.21
0
0
0
0
1
0
0
0.3221
0
0.1494
11.15
0.4133
0
0
19.27
.Figure 6
Stabilization scheme for the boiler without direct feedthrough.The design parameters values are chosen as follows:(43)
R
=
100
0
0
0
0
0
0
0
0
0
2
0
0
0
0
130
,
Q
=
0.4
0
0
0
0
0
0
0
0
0
0.5
0
0
0
0
2.006
.And N
1
=
1, N
2
=
15, and N
u
=
2.Figure7 presents the obtained results when a population size of 5 particles and 5 iterations is used. A zoom of the area lying between the 17 and 35 minutes of the steam flow rate is given by Figure 8. The corresponding control signals are depicted in Figure 9.Figure 7
Response of the boiler using the three controllers (red solid line: proposed controller; blue dotted line: conventional PSO based controller; black dashed dotted line: linear MPC; and blue dashed line: reference trajectory).Figure 8
Zoomed area between the instants 17.5 and 35 min of the steam flow rate.Figure 9
Control signals for the three controllers (red solid line: proposed controller; blue dotted line: conventional PSO based controller; and black dashed dotted line: linear MPC).The obtained results show that the proposed controller outperforms the linear MPC and the conventional PSO based controller. The outputy
2 is almost identical for the three controllers, due to the fact that the FAR is always kept constant regardless of the used controller.To compare the performance of the considered algorithms when different values of the population size and the number of iterations are used, the Mean Cost Value (MCV) criterion given by (44) is adopted: (44)
MCV
=
1
N
∑
t
=
1
N
J
U
t
,
N
is the number of samples
.For each value of the population size and the number of iterations both control algorithms are executed fifteen times and the corresponding value of the MCV criterion is computed. The averages values of the MCV criterion, over the fifteen executions, corresponding to each population size and iterations number are depicted in Figure 10. It is clear, from this figure, that the proposed PSO based controller always outperforms the conventional PSO based controller, regardless of the selected population size and the number of iterations.Figure 10
MCV average values.The execution time of a control algorithm is a very important parameter to evaluate its computing efficiency and real-time implementation feasibility. Table3 gives the execution time of the considered control algorithms for different population size and number of iterations; they were executed on an Intel Core i5 3.10 GHz based machine.Table 3
Execution time of the considered algorithms.
Configuration (population × iteration)
Proposed PSO based controller (ms)
Conventional PSO based controller (ms)
5 × 5
223.26
207.86
5 × 10
398.42
385.98
5 × 15
572.94
560.08
5 × 20
746.91
736.00
5 × 25
922.58
916.40
10 × 10
771.95
768.84Table3 indicates that the execution time of the proposed algorithm is slightly higher than that of the conventional algorithm. Although, when the population size or the iteration number increases, the execution times of both algorithms get closer. It was expected to get slightly higher computing time with the proposed controller as it requires more operations than the conventional PSO based controller. However, this small increase in the computation requirement has produced a decent performance increase.


## 6.1. Process Description

To assess the performance of the proposed control algorithm, the control of the industrial boiler model developed by Pellegrinetti and Bentsman [28] is considered. This model represents a dual fuel (oil and gas) boiler in the Abbott Power Plant in Champaign, Illinois, which could be used for both heating and electric generation. It is a MIMO process, which has four manipulated inputs (fuel flow rate, air flow rate, feed water flow rate, and a valve to regulate the steam demand) and four measured outputs (pressure, oxygen level in the flue gas, drum water level, and steam flow). The control objectives are to maintain the pressure, the level of the water in the drum, and the oxygen level in the flue gas at their desired values, regardless of the steam flow rate provided by the boiler and any other disturbances, such as the fluctuations in the heating value of the fuel or the variations in the ambient temperature. This mathematical model includes perturbations and measurement noises and is described by the following equations:(32)
x
˙
1
t
=
c
11
x
4
t
x
1
9
/
8
t
+
c
12
u
1
t
-
τ
1
-
c
13
u
3
t
-
τ
3
+
c
14
,
x
˙
2
t
=
-
c
21
x
2
t
+
c
22
u
2
t
-
τ
2
-
c
23
u
1
t
-
τ
1
-
c
24
u
1
t
-
τ
1
x
2
t
c
25
u
2
t
-
τ
2
+
c
26
u
1
t
-
τ
1
,
x
˙
3
t
=
c
31
x
1
t
-
c
32
x
4
t
x
1
t
+
c
33
u
3
t
-
τ
3
,
x
˙
4
t
=
-
c
41
x
4
t
+
c
42
u
1
t
-
τ
1
+
c
43
+
c
41
u
4
t
,
y
1
t
=
c
51
x
1
t
-
τ
4
+
n
1
t
,
y
2
t
=
c
61
x
2
t
-
τ
5
+
n
2
t
,
y
3
t
=
c
70
x
1
t
-
τ
6
+
c
71
x
3
t
-
τ
6
+
c
72
x
4
t
-
τ
6
x
1
t
-
τ
6
+
c
73
u
3
t
-
τ
3
-
τ
6
+
c
74
u
1
t
-
τ
1
-
τ
6
+
c
75
x
1
t
-
τ
6
+
c
76
1
-
c
77
x
3
t
-
τ
6
x
3
t
-
τ
6
x
1
t
-
τ
6
+
c
78
+
c
79
+
n
3
t
,
y
4
t
=
c
81
x
4
t
-
τ
7
+
c
82
x
1
t
-
τ
7
+
n
4
t
,where x
˙
1 is the drum pressure state K
g
f
/
c
m
2, y
1 is the measured drum pressure (PSI), y
2 and x
˙
2 are the measured excess oxygen level and its state, respectively (%), x
˙
3 is the system fluid density K
g
/
m
3, y
3 is the drum water level (in), y
4 is the steam flow rate K
g
/
s, u
1, u
2, and u
3 are, respectively, the fuel, air, and feed water flow rates which take values in the interval 0
1, and x
˙
4 is the exogenous variable related to the steam demanded. The variables n
i are the outputs of the first-order colored noise models driven by zero mean and unit variance white noise and are considered to be unmeasured output disturbances. The remaining coefficients values are given in Table 1. Consider(33)
n
1
=
0.75
s
+
0.1
s
+
0.001
w
1
,
n
2
=
0.019
s
+
0.001
s
+
0.024
w
2
,
n
3
=
0.105
s
+
0.038
s
+
0.01
w
3
,
n
4
=
0.01
s
+
0.0001
s
+
0.001
w
4
,where n
i, i
=
1
,
…
,
4, are colored noise and w
i, i
=
1
,
…
,
4, is the unit variance white noise.Table 1
Model equations coefficients.
c
11
=
-
0.00478
c
12
=
0.28
c
13
=
0.01348
c
14
=
0.02493
c
21
=
0.1540357
c
22
=
103.5462
c
23
=
107.4835
c
24
=
1.9515
c
25
=
29.04
c
26
=
1.824
c
31
=
0.00533176
c
32
=
0.025195
c
33
=
0.7317058
c
41
=
0.04
c
42
=
0.0299886
c
43
=
0.018088
c
51
=
14.214
c
61
=
1.00
c
81
=
0.85663
c
82
=
-
0.18128
c
70
=
-
0.1048569
c
71
=
0.15479
c
72
=
0.4954961
c
73
=
-
0.20797
c
74
=
1.2720
c
75
=
-
324212.7805
c
76
=
-
99556.24778
c
77
=
0.001185
c
78
=
-
1704.50476
c
79
=
-
103.7351
τ
1
=
2
τ
2
=
2
τ
3
=
3
τ
4
=
3
τ
5
=
4
τ
6
=
10
τ
7
=
2The process is supposed to be initially operating around the following states:(34)
x
0
=
22.5,2.5,621.17,0.6941
T
,
y
0
=
320,2.5,0
,
9.3053
T
,
u
0
=
0.3227,0.39503,0.37404,0
T
.The linear representation of the boiler around the operating state given by (34) is described by(35)
x
˙
=
A
b
x
+
B
b
u
,
y
=
C
b
x
+
D
b
u
,where(36)
A
b
=
-
0.005508
0
0
-
0.15872
0
-
0.20625
0
0
-
0.012156
0
0
-
0.566887
0
0
0
-
0.04
,
B
b
=
0.28
0
-
0.01348
0
-
9.374893
7.658411
0
0
0
0
0.731706
0
0.029989
0
0
0.04
,
C
b
=
14.214
0
0
0
0
1
0
0
0.322072
0
0.149410
11.148662
0.413307
0
0
19.274175
,
D
b
=
0
0
0
0
0
0
0
0
1.272
0
-
0.20797
0
0
0
0
0
.There is a direct feedthrough between the outputy
3 and the inputs u
1 and u
3, expressed by a nonzero D matrix in the state representation. This property will create some problems when considering the linear MPC controller [49]. In order to have D
b = 0, a delay will be introduced to the inputs u
1 and u
3 with regard to the output y
3.This process is inherently unstable. Essentially, a stabilisation scheme of the water level must be introduced to make any identification approaches possible by incorporating a proportional feedforward action (C2) of 0.0403 and a PI controller (C1) withK
p
=
0.258 and T
i
=
1.1026
e
-
4, as it is illustrated in Figure 4.Figure 4
Stabilization scheme for the boiler.


## 6.2. Process Fuzzy Identification

The first step in the design of a fuzzy NMPC controller is the construction of a nonlinear model of the process. Four fuzzy multi-inputs single output models were constructed to predict the future behaviour of the process. Two triangular membership functions per input, a sampling period of 3 s, and the models given by the following equations were used:(37)
y
^
1
=
f
1
y
^
1
k
-
1
,
u
1
k
-
5
,
x
^
4
k
-
3
,
y
^
3
=
f
2
y
^
3
k
-
1
,
u
1
k
-
12
,
r
Level
k
-
13
,
u
4
k
-
10
,
x
4
k
-
10
,
y
^
4
=
f
3
y
^
4
k
-
1
,
u
1
k
-
4
,
u
4
k
-
2
,
x
4
k
-
2
,
x
^
4
=
f
4
x
^
4
k
-
1
,
u
1
k
-
2
,
u
4
k
-
1
,where f
i, i
=
1
,
…
,
4, represent the fuzzy relationships between the inputs and the outputs.Figure5 shows the validation results of these models.Figure 5
Validations of the fuzzy models.


## 6.3. Controller Implementation and Simulation Results

It has been shown in [28] that there is a nearly linear relationship between the fuel/air ratio (FAR
=
u
2
/
u
1) and the oxygen level. So, the oxygen level y
2 can be kept constant at a desired value by maintaining the corresponding FAR value constant. Given that u
1 affects all the systems outputs, varying u
2 to maintain the FAR value constant is the appropriate way. Thus, the value of u
2 can be directly deduced from that of u
1. To keep the desired oxygen level in flue gas y
2 around 2.5% the FAR value must be equal to 1.2241 and u
1 equals 0.8169 u
2.The control objective is to force the steam flow ratey
4 to track a desired reference trajectory, while fulfilling the following constraints on the manipulated inputs,(38)
Δ
u
1
Max
⁡
=
0.1
,
Δ
u
2
Max
⁡
=
0.1
,
Δ
u
r
Level
Max
⁡
=
1
-
4
<
r
Level
<
4
,
Δ
u
4
Max
⁡
=
0.1
-
0.2
<
u
4
<
0.35
,
0
<
u
3
<
1
,and the following constraints on the outputs,(39)
y
1
=
320
,
y
2
=
2.5
,
y
3
=
0
.The following parameters values were used in implementing the proposed control algorithm:(40)
R
=
100
0
0
0
0
0
0
0
0
0
2
0
0
0
0
130
,
α
1
=
1.65
e
-
4
0
9.54
e
-
5
3
e
-
3
0
0
0
0
1.65
e
-
3
0
1.05
e
-
2
2.76
e
-
2
1.28
e
-
4
0
1.246
e
-
4
3
e
-
3
.And N
u
=
2, number of particles = 5, number of iterations = 5, and α
2
=
0
0
0
0
T. The values of the other parameters are gathered in Table 2.Table 2
Design parameters.
N
1
N
2
ε
i
Γ
0
y
min
⁡
y
max
⁡
y
^
1
1
10
1
0.0085
318
322
y
^
3
10
20
0.7
0.0065
−1.5
1.5
y
^
4
1
10
0
1.006
/
/
x
^
4
1
10
0
0
/
/In addition to the proposed algorithm, two other control strategies were considered:Using the Conventional PSO. In this algorithm, the particles are randomly distributed in the whole search space and the solution of the previous sampling time is taken as a possible solution for the current sampling time. The PSO parameters values χ
=
0.729843788 and c
=
2.05, adopted by Clerc and Kennedy [34], are used in both PSO based controllers.Using the Linear MPC Strategy. In this strategy, only the constraints on the inputs variables and their increments are considered and the optimization problem is converted into a quadratic programming (QP) problem.To implement this strategy, a linear representation for the process given in Figure4 must be used. To remove the direct feedthrough between the output y
3 and inputs u
1 and u
3, the structure given in Figure 4 is slightly modified by introducing a unit delay in these inputs (Figure 6). The obtained system is linearized around the operating state given by (34) and then discretized using a sampling period of 3 s. The following discrete state space representation is obtained:(41)
x
˙
k
+
1
=
A
x
k
+
B
u
k
,
y
k
+
1
=
C
x
k
,where(42)
A
=
0.9861
0
0.001485
-
0.3706
0
0.5386
0
0
-
0.1736
0
0.9187
-
5.665
0
0
0
0.8869
,
B
=
0.8172
0
-
0.009937
-
0.02264
-
20.97
17.13
0
0
-
0.3384
0
0.5441
-
0.3525
0.08478
0
0
0.1131
,
C
=
14.21
0
0
0
0
1
0
0
0.3221
0
0.1494
11.15
0.4133
0
0
19.27
.Figure 6
Stabilization scheme for the boiler without direct feedthrough.The design parameters values are chosen as follows:(43)
R
=
100
0
0
0
0
0
0
0
0
0
2
0
0
0
0
130
,
Q
=
0.4
0
0
0
0
0
0
0
0
0
0.5
0
0
0
0
2.006
.And N
1
=
1, N
2
=
15, and N
u
=
2.Figure7 presents the obtained results when a population size of 5 particles and 5 iterations is used. A zoom of the area lying between the 17 and 35 minutes of the steam flow rate is given by Figure 8. The corresponding control signals are depicted in Figure 9.Figure 7
Response of the boiler using the three controllers (red solid line: proposed controller; blue dotted line: conventional PSO based controller; black dashed dotted line: linear MPC; and blue dashed line: reference trajectory).Figure 8
Zoomed area between the instants 17.5 and 35 min of the steam flow rate.Figure 9
Control signals for the three controllers (red solid line: proposed controller; blue dotted line: conventional PSO based controller; and black dashed dotted line: linear MPC).The obtained results show that the proposed controller outperforms the linear MPC and the conventional PSO based controller. The outputy
2 is almost identical for the three controllers, due to the fact that the FAR is always kept constant regardless of the used controller.To compare the performance of the considered algorithms when different values of the population size and the number of iterations are used, the Mean Cost Value (MCV) criterion given by (44) is adopted: (44)
MCV
=
1
N
∑
t
=
1
N
J
U
t
,
N
is the number of samples
.For each value of the population size and the number of iterations both control algorithms are executed fifteen times and the corresponding value of the MCV criterion is computed. The averages values of the MCV criterion, over the fifteen executions, corresponding to each population size and iterations number are depicted in Figure 10. It is clear, from this figure, that the proposed PSO based controller always outperforms the conventional PSO based controller, regardless of the selected population size and the number of iterations.Figure 10
MCV average values.The execution time of a control algorithm is a very important parameter to evaluate its computing efficiency and real-time implementation feasibility. Table3 gives the execution time of the considered control algorithms for different population size and number of iterations; they were executed on an Intel Core i5 3.10 GHz based machine.Table 3
Execution time of the considered algorithms.
Configuration (population × iteration)
Proposed PSO based controller (ms)
Conventional PSO based controller (ms)
5 × 5
223.26
207.86
5 × 10
398.42
385.98
5 × 15
572.94
560.08
5 × 20
746.91
736.00
5 × 25
922.58
916.40
10 × 10
771.95
768.84Table3 indicates that the execution time of the proposed algorithm is slightly higher than that of the conventional algorithm. Although, when the population size or the iteration number increases, the execution times of both algorithms get closer. It was expected to get slightly higher computing time with the proposed controller as it requires more operations than the conventional PSO based controller. However, this small increase in the computation requirement has produced a decent performance increase.


## 7. Conclusion

A particle swarm optimization algorithm that uses the Gaussian distribution was proposed to solve the constrained fuzzy predictive control problem. A Takagi-Sugeno fuzzy model was used to predict the future values of the system outputs. Fuzzy systems are known for their simplicity and their ability to approximate complex nonlinear processes. However, when using a nonlinear model the optimization problem to be solved is then a nonlinear and a nonconvex one. The solution of this problem using numerical methods is expensive in computational time, especially when considering constraints. Particle swarm optimization algorithms are powerful metaheuristic based algorithms that require much less computational effort than the genetic algorithms, and their implementation is relatively easy. In this work, the PSO algorithm was adapted to find an efficient solution to the fuzzy predictive control optimization problem with a low computational burden. This was done essentially by reducing both the PSO population size and the number of iterations needed to obtain a suitable solution.The efficiency and the control accuracy of the proposed algorithm were investigated and compared to other control strategies by considering the control of a highly nonlinear MIMO process, namely, an industrial boiler. The obtained simulation results have proved the efficiency and the control accuracy of the proposed algorithm, when considering highly nonlinear systems.



---

*Source: 437943-2015-06-28.xml*