Informatica 35 (2011 ) 363-374 363
Multivariable Generalized Predictive Control Using an Improved Particle Swarm Optimization Algorithm
Moussa Sedraoui and Samir Abdelmalek
Laboratoire (PI : MIS) Problèmes Inverses: Modélisation, Information et Systèmes, Université 08Mai 1945, Guelma, Algérie E-mail: msedraoui@gmail.com
Sofiane Gherbi
Laboratoire d'Automatique de Skikda (LAS), Université 20 Août 195 5, Skikda, Algérie E-mail: sgherbi@gmail.com
Keywords: improved particle swarm optimization, multivariable generalized predictive control, feasible region Received: November 16, 2010
In this paper, an improvement of the particle swarm optimization (PSO) algorithm is proposed. The aim of this algorithm is to iteratively resolve the cost problem of the Multivariable Generalized Predictive Control (MGPC) method under multiple constraints previously reduced. An ill-conditioned chemical process modelled by an uncertain Multi-Input & Multi-Output (MIMO) plant is controlled in order to verify the validity and the effectiveness of the proposed algorithm. The performances obtained are compared with those given by the MGPC method using the standard PSO algorithm. The simulation results shows that the proposed algorithm outperforms standard PSO algorithm in terms of performance and robustness.
Povzetek: Predstavljena je metoda optimiranja z roji za nadzor s splosnim napovedovanjem in vec spremenljivkami.
1 Introduction
Multivariable generalized predictive control (Morari & Lee, 1999) is a very powerful method. It has been the subject of many researches during the last few years and it was applied successfully in industry, particularly in chemical processes. It is based on MIMO predictive model [1], [2] where the expected behaviour of the system can be predicted in the extended time horizon. The MGPC law is obtained by minimizing linear or nonlinear criterion (Magni, 1999, Duwaish, 2000). This criterion is composed by the sum of the square prediction errors between the predicted and desired outputs, the weighted sum of the square change-controls (control-increments) and others [3], The constraints inclusion (as mathematical inequalities type) distinguishes most clearly MGPC from other process control paradigms as suggested in (Richalet, 1993, Qin 1997, Rawlings, 1999). These constraints are imposed in order to ensure a better stability and performance robustness (Al Hamouz and Duwaish, 2000, Imsland, 2005). The MGPC method formulates the constraint optimization problem at every step time for solving the optimal control move vector [4], At the next sampling time, a new process measurement is received, the process is updated, and a new constraint optimization problem is solved for the next control move vector. An efficient randomized constraint optimization algorithm is suggested to the MGPC method named by
PSO algorithm (Rizvi & al, 2010, Yousuf & al, 2009, Al Duwaish, 2010). This algorithm explores the search space using a population of particles, each one with a particle or an agent, starting from a random velocity vector and a random position vector. Each particle in the swarm represents a candidate solution (treated as a point) in an «-dimensional space for the constraint optimization problem, which adjusts its own "flying" according to the other particles [5], The PSO algorithm can resolve successively various constraint optimization problems, such as linear or non-linear, convex or non-convex problems. Unfortunately, it cannot provide satisfactory results when the MGPC method is applied to poorly modelled processes [6] operating in ill-defined environments. This is, as often, the case when the plant has different gains for the operational range designed by user's trial-and- error. In addition, the PSO algorithm's convergence cannot satisfy multiple time domain specifications if the process (to be controlled) is constrained by a high number of hard constraints (Leandro dos Santos Coelho & al, 2009). Several heuristic algorithms have been developed in recent years to improve the performance and set up the parameters of the PSO algorithm [7],This paper investigates the analysis of the above mentioned problems. Two main contributions are proposed in this paper in order to
366 Informática 35 (2011 ) 363-374
M. Sedraoui et al.
improve the performances of the MGPC method. The first one consists to reduce (if possible) the imposed inequality constraints which are reformulated as boundary constraints. The second one is to resolve the bounds constraints optimization problem by the improved PSO algorithm.
2 Unconstrained MGPC Method
All the considered matrices are in discrete time domain.
A CARIMA (Controller Auto Regression Integrated Moving Average) model for an m inputs and m outputs multivariable process can be expressed by [8]:
A(q-l)ty(t):=B(q-l)M(t-\) + C(q-l)C{t) (1) Where
y2(t) ... ym(t)T
u2(t) ... Um(t)f
C(q~l) are mxm monic polynomial matrices. Set C(q ') equal to the unity diagonal matrix. C(t) is an uncorrected random process
and A(q ') = 1 - f/ '. this form enables to introduce an integrator in the control law. Without lost of generality one can supposed as diagonal polynomial matrix.
}) (t) e 91, ui (t) e 91 denotes respectively, the process output and the control input of the channel number '/'. q ' denotes the backward shift operator. The
role of A(q ') is to ensure an integral action of controller in order to cancel the effect of the step varying output in the channel'/'.
As in all receding horizon predictive control strategies, the control law provides that, for each channel '/', the control-increment A;/( (t) which minimizes the following unconstraint cost problem of the MGPC method [8]:
y(t)eKmA-.= [yi(t) M(0e9Txl :=[m!(0 A(q~l), B(q~l) and
( ;
¿=1
Zhc-
j=i
j It)- w, (t + j)]2 + A, £ [Am, (t + j- l)]2
J=l
(2)
Where
yt (t + j) e 91 is an optimum /-step-ahead prediction of the system output vector on data up to time t. therefore, the expected value of the output vector at time t if the past input vector, the output vector, and the future control sequence are known. Noting that i',0 + ./) is depending to the control-increment Aui from resolving two Diophantine equations (more details are available in the reference [9]).
\r( (t) e 91 is the future set-point or the reference sequence for the output y, (t).
N'2, A; (wiili respect: A ; A \) denotes respectively, the maximum output prediction horizon (assumed equal to \'2 e ) and the maximum control prediction (assumed equal to Y„ e ) for each channel
'/'. Ij e 9l+ denotes the positive parameter weighting the control input for each channel'/'.
3 Classification of Constraints and Problem Formulation
In constrained control, a set of inequality constraints may be set as addition of the control objective and the variation limits of certain variables to the given ranges:
v; <v1(t + j)<vt, with i:=\,2,..m and
j:=NsU...,Ns2. Where
Vj (t + j) e 91 is a variable under restriction,
e 91 and vi e 91 are the lower and the upper boundaries,
A'vl and A'v2 are the lower and the upper constraint horizons respectively.
The two main objectives of constrained predictive control are set-point tracking and prevention / reduction of constraint transgressions. These constraints can be imposed (with respect to the time index) on the control-increment vector, or/and on the control vector as follows: - Constrained on the control-increment: Auj < AUj (t + j) <	(3)
Where / = 1,2,..m and j = 0,...,NU -1. - Constrained on the control:
«,. <u,(t + j) <u,	(4)
Where / = 1,2,..m and j = 0,...,NU -1. By using:
j
Uj (t + j) := Uj (t -1) + ^ AUj (t + k)
k=0
The control constraints (4) becomes as follow: j _
m, - u, (t -1) < ^ A Uj (t + k)<Ui- Uj (t -1)
(5)
(6)
k=0
The constraints on the control vector and the rate of control changes, with respect to the batch index, can be easily combined together:
A^A U<Blnq	(7)
Where
( A U(t) ^ AU (/ + 1)
A U,
(m-Nu)yl ■
denotes the design
KAU(t + Nu -l)y
parameter vector which will be determined later by the PSO algorithm, it contains the future control- increment vector (pjj{t + j))ni of each channel as:
i Au.it + j)^
A U(t+ ])-.=
A u2(t + j)
A«m(t + j)
y j=ox—,Nu-i
MULTIVARIABLE GENERALIZED PREDICTIVE.
Informatica 35 (2011 ) 363-374 365
(Anq )
' (Pinq )(,
-i„q }(4.m.NuHm.Nu, , y^inq 4mJV )xl arc defined by: ( diag(Imym) ^
A •=
inq •
-diag(Il
trWm>
— tri 1(1 n Where
diag(Imxm)e 91' matrix, and tril(I
)
{m.Nu)y{m.Nu)
denotes the unity diagonal
Je«
(.m.Nu)x(m.Nu)
triangular matrix of the unity diagonal matrix (Imxm).
D ._
inq ~
[Ait
i(m-Nu)xl
\p-Ui \m-Nu)-A 1)1
as:
J(AU,t) :=AUt-Q2-AU + Q[-AU + Q0
Where 02 := GT ■ G + A, Q? := 2(r - Wf G and
0o :=(r-w)T .(r-w)
A:=;i-I(m-Nu)y<m-Nu)is diagonal matrix weighting control-increment vector, and
the
^(m-Ai2)xiis the
min J (AU, t) := A UT -Q2 - A U + Q? ■AU + Q0
AU
Is.t:
Alnq-AU<Blnq
(9)
propose a systematic method that determines the minimum set vector of limiting constraints. The lower and the upper bounds of the feasible region are given as below:
For each channel / := L-••.//; . the control-increment AUj (t) is simultaneously constrained by: 1- For the control prediction horizon j = 0 :
denotes the lower
(10)
J(m-Nu )xl
The cost index (2) can be expressed in matrix form
(8) —
(14)
projected set-point vector.
G(m-N2)y(m-Nu), T(m-N2)y\ are the polynomial matrices
which are determined by the recursively resolution of the two Diophantine equations [9],
The cost index (8) and the inequality constraints (7) formulate the following constraint optimization problem as:
Au- < Aui (t) < Aui
- u(t -1) < Aut(t) <ut-u(t-1)
It is easy to see that the new lower and upper bounds are determined by:
Vj(t) < Aut(t) < Vt(t)	(11)
Where
v; (t) \= max{AM; ut -u(t-1)}	(12)
Vj(t) := min{AM, ut-u(f-1)}	(13)
2- For the control prediction horizon / = 1 :
Au. < A ui (t +1) < A ut
\ut -u(t-1) < AUj (t +1) + AUj(t)<ui-u(t-1)
The new lower and upper bounds are determined for j = lby:
v_i(t + \)<Aul(t + \)<v1(t + \)	(15)
Where
v.(? + l):=max{A«. -m(/-1)}-v;(/)}	(16)
v^t + l) := min{AM, {ut-u(t-\)}-Vi(t)}	(17)
This procedure is repeated until the control prediction horizon j = NU -1. Therefore the control-increment Aui (t + Nu-1) is constrained by the new bounds:
Nu-2
vi(t + Nu -1) :=max{Auf -u(t-1)}- ^v^t + k)}
Now, an optimal control vector is given by the PSO algorithm. This algorithm should minimize the objective function (8) under4xmxNu inequality constraints. The computational requirements of the PSO algorithm depend heavily on the number and the type of the constraints to be satisfied. An efficient off-line constraint PSO algorithm, suggested by Ichirio & al, 2009, can resolve this problem [10], Unfortunately, this algorithm is difficult to extend to the MGPC method for two reasons: the first one is due to a large dimension of the inequality constraints which needs excessive computation time. The second one is due to a real-time output feedback implementation of the MGPC method which requires a minimum consuming time. To resolve these above problems, the inequality constraints should be reduced, for each step time, and reformulated as bounds constraints type. Only those constraints (which limit the feasible region) must be taken into account. The efficiency of the PSO algorithm can be increased if the superfluous constraints (which do not limit the feasible region) should be eliminated [11]. In this paper we
k=0
N.,-2
(18)
y, (t + Nu-1) := min{Aut {ut -u(t-1)}- ^v, (t + k)}
k=0
(19)
Then, for each step time value t0,tl,---,tm3K , the feasible region D(i,j,t):=(v^(t + j) vt(t + j))can be
determined by the following proposed algorithm:
3.1 Reduced constraints algorithm
For each point time .. - - -. /m;ix , the feasible region is determined by the folio wings steps:
[Step 1]: Set the first counter i <— 1 which denotes the number of channels, and go to the next step. [Step 2]: Set the second counter j <— 0 which denotes the control horizon prediction, and go to the next step. [Step 3]: Set the parameters hmsK <— 0, hmm <— 0 and go to the next step.
[Step 4]: Build the followings ranges: bound_max, := {Am, {ut — m^(f — 1)} — Amax} .
366 Informática 35 (2011 ) 363-374
M. Sedraoui et al.
bound _min; := {Am, {m; - ui (t -1)} - hmm }.
[Step 5]: Calculate the new upper and the new lower
bounds which limit the control- increment Aui (t + j) by:
vi (t + j) := vam{bound _ max,}
vi(t + j):=max{bound_mmi}
From these above bounds, the feasible region is determined as follow:
D(i, j,t) := (vjXt + j) Vt(f + j))
[Step 6]: Update the parameters: ^ ; hmxi as follows
^max	x+vi(t + j),
^min ^min + vi (f + j) > and go to the next step [Step 7]: Update the second counter j <— j +1 a«<i go back to the step 4 if j -<NU -1. Otherwise, go to the next step.
[Step 8J: Update the first counter / <— / +1 and stop algorithm if i := m . Otherwise go back to the step 2. From this above algorithm, the constraint optimization problem (9) under inequality constraints is reformulated as the bounds optimization problem:
\mmJ(AU,t) := AUT -Q2 -AU + Qf •AU + Q0
j m _ (20)
[s.t: AU<AU<AU
Where AU , AU denotes respectively, the new lower and the new upper bounds vector which limit the feasible region D(m.Nu)x2 := (AU,AU), with:
(ML\m.Ky i:=k(0 v,(' + D - v^t + N,- l))f=12...m
:= f~(/) ^(/ + 1) V^t + N,-l)J1=l2...m
From (20), it is easy to see that the inequality constraints number is reduced to mxNu constraints at each step time. This dramatic reduction has a capital importance for the success of the PSO algorithm. Now, we are able to find the optimal control of the MGPC law. The new constraint optimization problem (20) should be resolved for each step time t := ill.t].---.imc. . its solution vector AU* denotes the optimal design parameter vector. Only the first m rows
of AU* is used to obtain the optimal desired control-increment vector of each channel (' / '). The optimal control vector is obtained by adding the previous control vector to the optimal control-increment vector as follow:
Uj(t) := Uj(t — l) + Am*(/)	(21)
4 Improved PSO Algorithm [6]
Particle swarm optimization algorithm, introduced first by Kennedy and Eberhart in (1995), is one of the modern heuristic algorithms which belong to the category of Swarm Intelligence method (Kennedy, 2001). The PSO algorithm uses a swarm consisting of N e X particles for
each control-increment vector (Aui (t + /)) , n, _ v , to
get an optimal solution Aui (I + j) which minimizes the optimization problem (20). The position of ( ith ) particle and its velocity are respectively denoted as [12]:
A«i,N(t + j))r
Au, (t + j) := Aw, j (t + j) Au,2(t + j)

( t + j)\=\j/ul(f + j) vU2(t + j)

M+J)]
Then, the position of the (ith) particle, Aui (t + j), is based on the following update law:
for £ =l,2,...£max, which indicates the iteration number [12],
^ :=Co^(ff^-Auf)+c2rl(h^ - Am/)
(22)
Au-+1 := Am/ + y/f+1	(23)
Where cl and c2 are respectively, the cognitive (individual) and the social (group) learning rates and are both positive constants. The value of cognitive parameter Cj signifies a particle's attraction to a local best position based on its past experience. The value of social parameter c2 determines the swarm's attraction towards a global best position.
c0 e 9l+ is the inertia weight factor whose value decreases linearly with the iteration number (Shi & Eberhart, 1999) as [13]:
C0 -- ^max
in _n \
max min
(
V	max
and 0„
■I
(24)
Where 0„,,,, and ()„.„„ are the initial and the final
IIldA	Illlll
values of the inertia weight, respectively. The values of ()mc. = 0.9 and 9mm = 0.4 are commonly used [13],
The random numbers r({ and r2J are uniformly distributed in [0,l].
H-estJ andhh^Jm denotes respectively, the best previously obtained position of the ( ith ) particle (the position giving the lower value of the objective criterion) and the best position in the entire swarm at the current iteration £ [10]:
H
best,I
h
bestJ
:= argmin{J(AM/ ),0 < r < £} := argmin{J(AM/),V/'}
Auf
(25)
(26)
From equation (23), some current position of ( i,h ) particle (in each dimension) can exceed the corresponding lower bound or upper bound of the feasible region. Consequently, the given optimal control vector of the MGPC method cannot satisfy some specifications and also some constraints are non-satisfactoriness' in some range time. To avoid, we should improve the convergence of PSO algorithm by adjusting only the corrupted position of ( ith ) particle with the region around the current established solution, if it is too smaller than the corresponding lower bound, its value vt should be replaced. If it is too higher than the
MULTIVARIABLE GENERALIZED PREDICTIVE.
Informatica 35 (2011 ) 363-374 365
corresponding upper bound, then its value is replaced by r, . The proposed modification can be formulated as follows:
Let consider:
Auf (t + j) : The corrupted position of ( /,/, ) particle given at current iteration I := q.
y( (/ + j), Vj (t + j) : The lower bound and upper
bound which are determined by the reduced constraints algorithm. So that, the above corrupted position can be adjusted by using the following inequalities:
fv,- (f + j) : if A uf (t + j) -cv,. (/ + j)
5 Simulation Results and Discussion
In this section, a multivariable generalized predictive control method using a modified particle swarm optimization algorithm is applied to a distillation column which is MIMO plant with two input and output vectors (benchmark problem, see [14]). The two inputs are the reflux and the vapour boil up rate and the outputs are the distillate and the bottom product. The results are compared with those given by the MGPC method using the standard PSO algorithm. The mathematical model is given by [14]:
0.878 -0.864
Auf (t+ ]):=<
\Vl(t + j):if A uf(t + j)yVl(t + j)
(27)
Consequently, from the equation (23), the current velocity should be limited by the following bounds:
v; - All?-1 <iy? <vt- Auf~1	(28)
Now, the modified current positions with their modified velocity are used to improve the next best position and their velocity vector for the next iteration as follow:
= c0vf + c\r\j	- Au? )+ c2rl te* - Au? ]
G(s) :=
Ku
1
755 + 1
1.082 -1.096
0
:[0.8 1.2],
'¿(¿=1,2)
[0.0 l.o].
0
K2e-(31)
Auf+l := Auf + y/f+1
(29)
(30)
The improved PSO algorithm consists of the following steps:
4.1 Proposed algorithm
For each step time t :=	m;i,. the optimal control-
increment is determined by the following steps:
[Step 1]: Determine the lower bound and the upper
boundVjit + j) ,Vj(t + j) which are corresponding the design parameter [Aut (t +	^.
[Step 2]: Initialize random swarm positions and velocities:
initialize a population (array) of particles with random positions and velocities (array) from the search
domain D := (AU,AU).
Set the counter £ <— 1 and go to the next step. [Step 3]: Evaluate the objective criterion (20) and obtain Hfest'e, hbs^fm according to (25) and (26).
[Step 4]: Update of a particle's velocity and its position according to (22) and (23).
[Step 5]: Check each parameter of the particle's position by the following corresponding lower bound and
upper boundvt (t + j), vi (t + j). Replace only those
exceeding these above bounds.
[Step 6]: Update the counter £ <— £ +1 and go back to the step 3 if £ < £max . Otherwise, stop algorithm and take the best position vector as an optimal solution which minimize the constrained optimization problem (20).
¿(/=1,2)
Where
, Kj denotes respectively, the uncertainty temperatures and uncertainty gains of the process.
The time domain specifications are formulated, for the time range t e [0,400] minutes, as below:
a- For the first set-point reference vector: w = (I of , the first and the second output channels yl and y2 must satisfy [14]:
(51):	y1 (t) > 0.9 in more than 30 minutes.
(52):	>', (t)<\.\\ the maximum over-shoot corresponding the first output channel cannot exceed 11% for all range time t e [0,400].
(Sj): 0.99 < i ,(x)< 1.01 : the static error value cannot exceed 1% ((oo) - w1 (oo)| < 1%).
(54):	y2(t) < 0.5 : the maximum over-shoot of the second output channel cannot exceed to 50% for all range time t e [0,400].
(55):For	t—» oo : -0.01 < y2(t) < 0.01: the static error value cannot exceed 1%. From another word:
|_y2 (oo) - w2 (oo)| < 1% .
(56):Closed	loop stability.
(57):Control	signals should be limited by [-200 + 200].
(58):Control-increment	signals should be limited by [-12 +12],
For the set-point reference vector: w = (1 ())' . the sampling time Te = 1 minute is used to determine a CARIMA predictive model of the chemical process for two followings parameters cases [14]:
Kl2 =r12 = 1 and Kl =1.2 ,K2 =0.8,r12 =1.
b-The same previous time domain specifications should be satisfied for the second set-point reference
vector \r = (() l)7 corresponding to the low gains
direction K] = k2 = 0.8 and the same time delay
constants zl = z2 = 1.
366 Informática 35 (2011 ) 363-374
M. Sedraoui et al.
The MGPC method is tuned by choosing: (N'2,N'U,A;)i=u =(8,6,0.01) at time range I :=[().400] minutes.
For each step time: t:=t0,t1,---,400, the feasible region is determined from the following constraints:
-200 < u,(t + j)]=0,..,5 < +200 . ¿=1,2
12 \n (I ■ j).	12.
¿=1,2
From the reduced constraints algorithm (see section 3.1), these above inequality constraints are reduced in order to determine the search space D at each step time. The constrained optimization problem is resolved by standard
and improved PSO algorithms according to the following parameters:
Swarm size: N := 24 .
Maximum iteration: m;i,. :=100.
Cognitive and social learning rates: c, = c2 := 1 .
For the set-point reference vector: w = (1 o)T and
the parameter system's: K12 = z12 = 1, the figures 1.1 to
1.3 shows the results given by the MGPC method using the standard PSO algorithm (dashed curves), and the MGPC method using the improved PSO algorithm (line curves). The tablel summaries the results obtained by the two algorithms.
											
1.1 1 0.9 = 0.6 B- 0.4 0.2											
											
											
											
											
											
											
0 0.5 >. 0.2 B- 0 -0.2	30 50 100 150 200 250 300 350 400 time (minutes)										
		; \									
											
											
	D 50 100 150 200 250 300 350 400 time (minutes)										
Figure 1.1: Set-point tracking results with standard and improved PSO algorithms for w = (1 0)T and Kl2 = z] 2 = 1 .
Figure 1.2: Control effort results with standard and improved PSO algorithms for w = (I 0)T and Kl2=zl2=\.
MULTIVARIABLE GENERALIZED PREDICTIVE.
Informatica 35 (2011 ) 363-374 365
							
							
n v--^						»JWy-r--	
U ¡-\«A J V							
0	50	100	150	200	250	300	350	400
							
J							
▼ Vr-yj		.....					
Ki n							
0	50	100	150	200	250	300	350	400
Figure 1.3: Control-increment results with standard and improved PSO algorithms for u = (l o)T and Kl2 = z] 2 = 1.
Specifications	(SI)	(52)		(S3)	<S4)		(S5)	(S6)	(S7)	(S8)	Decision
<St):	y(30)			yi(-too)			y:(4<m	stable /	Unsatisfactory/ satisfactory	Unsatisfactory/ satisfactory	& reasons
Algorithms:		max(yi)	time		max(y2)	time		unstable	constraints	constraints	
Standard	1.010	1.074	23	1.010	0.343	8	0.006864	stable	unsatisfactory	unsatisfactory	Rejected
PSO									/ /<«,< 349.3 -9.9<m2<351.4	-44<Ju,<93 -44< \u ,<91	algorithm (S7).(S8)
Improved	0.990	1.096	40	0.9975	0.4846	13	0.002512	stable	Satisfactory	Satisfactory	Accepted
PSO									-20.1 <u j<200	-12<Au,<12	algorithm
									-20.3<u2<199	-12<Au2<12	
											
Table 1: Summary of the results (unsatisfactory performances are in bold) for the nominal model and the set-point reference \r = (l O)7 .
According to the figure 1.1, we can see that, the tracking dynamic of set-point reference vector found by MGPC method based on a standard PSO algorithm is better than the other algorithm but unfortunately, the time domain specifications: (S7) and (S8) are not satisfied.
In the figure 1.2, the obtained control signals of the MGPC method based on standard PSO algorithm exceed the constraint ranges at t := (5.6.••• 15} minutes such as:
Mimax (9) = 349 3 and ^2max (9) = 3514 • In addition, the control-increment signals presented in the figure 1.3 also violate the constraint ranges at times:
t := {(2 • • • 5), (7 • • • 11), (13 • • • 22)}	minutes.
Consequently, the performance robustness of this method is very poor in comparison with the MGPC method using the improved PSO algorithm which is capable to satisfy all time domain specifications. These results confirm the usefulness and the robustness of the proposed algorithm.
Figures 2.1, 2.2, 2.3 and table 2 give the results of the MGPC method with the following parametric changes in the process: (/,", = 1.2. K2 =0.8,r12 =l) for
the set-point reference vector w = (l o)r .
According to the figures 2.1 to 2.3, the better results are obtained by the improved PSO algorithm which satisfies all time specifications (Si to S8). These results can be explained by the best stability robustness against the process parametric disturbances. Furthermore, the control and the control-increment signals from the standard PSO algorithm show a dramatic oscillation at transient time region and exceed the constraint ranges. In fact, this algorithm cannot fulfill the three followings time domain specifications: (S2) withmax(Vj) = 11.124%, (S7) with
"lmax = 251-W2max = 374 and(S8)with -26.4 < AMj = 70.76 , -39.7 < Au2 = 99.05 , which can be explained by a high sensitivity to the parametric process variations. Thus, from these figures and table 2, we confirm the superiority of the proposed algorithm.
Figures 3.1, 3.2, 3.3 and table 3 give the results of the MGPC method using both algorithms when low gains directions of the process and set-point reference vector change simultaneously as follows:
(A = 0.8, K2 = 0.8, ru = l), w = (O if .
366 Informática 35 (2011 ) 363-374
M. Sedraoui et al.
Figure 2.1: Set-point tracking results with standard and improved PSO algorithms for w = (l Of and
(£1=1.2,£2=0.8,TU=1).
									
control effort u1 o o o o o o o o o	A								
									
	Y \			....					
						- -1			
									
									
	D 50 100 150 200 250 300 350 400 time (minutes)								
control effort u2 o o o o o o o o o o o									
									
	X \								
	/ \ V	- | .		...	rt m		___		
						v -			
									
									
	50 100 150 200 250 300 350 400 time (minutes)								
Figure 2.2: Control effort results with standard and improved PSO algorithms for w = (l Of and
(£1=1.2,£2=0.8,TU=1).
These above figures clearly show the performance superiority of the proposed PSO algorithm over standard PSO. For this case, the time domain specifications: S2, S7 and S8 are satisfied with the proposed PSO algorithm, while the same specifications are not satisfactoriness with standard PSO. In addition, the obtained outputs by the standard PSO algorithm converge to the set-point references but unfortunately, two other specifications cannot be satisfied at time t = 208 minutes which are:
(S3): |>-2(208) - w2(208)| = 5%.
(S5): ^(208)-Wl(208)| =3%. 6 Conclusion
In this study, we proposed an improvement of the PSO algorithm, it has been introduced and applied to solve the
constrained MGPC problem. In order to find a feasible region, the constraints on the controls and their increments have been previously reduced at each step time, the obtained convergences by improved PSO algorithm are well improved in comparison with the standard PSO algorithm. The efficient of the proposed algorithm is clearly shown and the performances robustness and the stability robustness are guaranteed with little still sensitivity to a set-point references changes and parametric model uncertainties. The results of the proposed algorithm justifies its efficiency and presents quite promising results and can be a subject of an interesting investigations.
MULTIVARIABLE GENERALIZED PREDICTIVE.
Informatica 35 (2011 ) 363-374 365
Figure 2.3: Control-increment results with standard and improved PSO algorithms for ir (l Of and
(A'j = 1.2,A"2 =0.8, Tj 2 =l).
Specifications	m	(S2)		(S3)	(S4)		(S5)	(S6)	(S7)	(S8)	Decision
(St):	yißO)			y ,(400)			y A4 m	stable	Satisfactory/	Satisfactory/	for
								/	unsatisfactory	unsatisfactory	reasons (Si)
		max(yi)	time		max(y2)	time		unstable	constraints	constraints	
Algorithms:											
Standard	1.061	1.1124	24	0.9953	0.2562	8	0.006142	stable	unsatisfactory	unsatisfactory	Rejected
PSO									-19.2<Ui<251	-26.4<du,<70. 76	algorithm
									-29.4<u2<374	-39. 7<dU2<99.05	(S2),(S7),(SJ
Improved	0.920	1.054	40	0.9953	0.493	12	0.006142	stable	Satisfactory	Satisfactory	Accepted
PSO									-3.5<Uj<133.5	-8.38<Au,<12	algorithm
									-6.4<u2<199	-12<Au2<12	
											
Table 2: Summary of the results (unsatisfactory pe
(A 1.2. A - = 0.8, 7] 2 = l)and the set-point reference w = (l
References
[1]	D.W Clarke, C. Mohtadi & P.S Tuffs, Generalized Predictive Control -Part I : The basic algorithm -Part II: extensions and interpretation. Automatica, 23(2), 1987, 137-160.
[2]	D.W Clarke & C. Mohtadi, Properties of generalized predictive control, Automatica, 25(6), 1989, 859-875.
[3]	P. Codron & P. Boucher, Improvement for multivariable generalized predictive control with multiple reference models. Proc. 12th LASTED International Conf. Modeling. Identification. And Control Innsbruck, Austria, Feb. 1993.
[4]	T.S. Chang & D.E. Seborg, A linear approach programming for multivariable feedback control with inequality constraints. International journal of control 37,1983, 583-597.
[5]	M. Zribi, M. Al-Rashed & M. Alrifai. Adaptive decentralized load frequency control of multi-area power systems. Electrical Power and Energy Systems, 27:575583, 2005.
[6]	J.P. Coelho, P.B.M. Oliveira & J.B. Cunlia. Greenhouse air temperature control using the
formances are in bold) for the uncertainty model Of.
particle swarm optimization algorithm. 15th Triennial world congress. Barcelona, Spain. IF AC, 2002.
[7]	L.D.S. Coelho & al. Model-free adaptive control optimization using a chaotic particle swarm approach. Chaos, solutions and fractals. Elsevier. 41, 2009, pp2001-2009.
[8]	P. Codron & P. Boucher. Multivariable generalized predictive control with multiple reference models: A new choice for the reference model, ECC'93. European Control Conf, June 28- July 1993.
[9]	P. Zelinka, B.R. Ilkiv & AG. Kuznetsov. Experimental verification of stabilizing predictive control. Control engineering practice7. 1999, pp 601-610.
[10]	I. Maruta, T. Hyoung & T. Sugie: Fixed-Structure H f Controller: A Meta-Heuristic Approach Using Simple Constrained Particle Swarm Optimization. pp:553-559. Automatica 45(2009).
[11]	E.F. Camacho. Constrained Generalized Predictive Control, IEEE Transaction on Automatic Control AC-38 (2), 1993, 327-332.
366 Informática 35 (2011 ) 363-374	M. Sedraoui et al.
time (minutes)
Figure 3.1: Set-point tracking results with standard and improved PSO algorithms for w = (o if and =0.8 ,K2 =0.8 ,zh2 =l).
time (minutes)
A							
1 \ ï / \							.....
							
							
							
							
0	50	100	150	200	250	300	350	400
time (minutes)
Figure 3.2: Control effort results with standard and improved PSO algorithms for w = (o if and =0.8,K2 =0.8,ru =l) .
[12]	H.N. Al-Duwaish, S.Z. Rizvi & al. PSO based Hammerstein modeling and predictive control of non linear multivariable boiler. Control Engineering practice, pp.1-12, 2010.
[13]	S. R. Singiresu: Engineering Optimization: Theory and Practice, Fourth Edition, by John Wiley & Sons, Inc. Hoboken, New Jersey. pp:708-714, & pp:761-768, 2009.
[14] D.J.N Limebeer. The specification and purpose of a controller design case study, IEEE Conf. Decision Control, Brighton, U.K, 1991, 1579-1580.
MULTIVARIABLE GENERALIZED PREDICTIVE.
Informatica 35 (2011 ) 363-374 365
i								
								
	V							
»I 1								
0	50	100	150	200	250	300	350	400
S								
i	A			i.A. . r				
	~AT -v.-*							
F								
0	50	100	150	200	250	300	350	400
time (minutes)
Figure 3.3: Control-increment results with standard and improved PSO algorithms for w = (0 if and
(A'j = 0.8, A"2 = 0.8, Tj 2 =l) .
Specifications	(SI)	(52)		(S3)	(S4)		(S5)	(S6)	(ST)	m	Decision
(SJ:	$2(30)			.Vi(-f00)			y 2(400)	stable	Satisfactory/	Satisfactory/	for
								/	unsatisfactory	unsatisfactory	reasons
		max(y2)	time		max(yi)	time		unstable	constraints	constraints	
Algorithms:											
Standard	1.046	1.167	21	0.9985	0.202	8	-0.003	stable	unsatisfactory	unsatisfactory	Rejected
PSO									-379.7<u,<56	-97.7<Au1<62.4	algorithm
									-383<U2<55.3	-99.2<AU2<62.5	(SJ.fSJ.fSs)
Improved	1.032	1.10	40	1.005	0.345	12	0.0009	stable	Satisfactory	Satisfactory	Accepted
PSO									-199<u1<20.3	-12<Auj<l 2'	algorithm
									-200<Ui<19.9	-12<Au:<12	
											
Table 3: Summary of the results (unsatisfactory performances are in bold) for the uncertainty model
(A = 0.8, A"2 = 0.8, t12 = l)and the change set-point reference w = (0 if.