Nonlinear Robust Identification Using Multiobjective ... - CiteSeerX
Nonlinear Robust Identification Using Multiobjective Evolutionary Algorithms J.M. Herrero, X. Blasco, M. Mart´ınez, and C. Ramos Predictive Control and Heuristic Optimization Group, Department of Systems Engineering and Control, Polytechnic University of Valencia [email protected], http://ctl-predictivo.upv.es Abstract. In this article, a procedure to estimate a nonlinear models set (Θp ) in a robust identification context, is presented. The estimated models are Pareto optimal when several identification error norms are considered simultaneously. A new multiobjective evolutionary algorithm −M OEA has been designed to converge towards ΘP , a reduced but well distributed representation of ΘP since the algorithm achieves good convergence and distribution of the Pareto front J(Θ). Finally, an experimental application of the −M OEA algorithm to the nonlinear robust identification of a scale furnace is presented. The model has three unknown parameters and ∞ and 1 norms are been taken into account.
1
Introduction
When modelling dynamic systems by first principles, the problem of parametric uncertainty always appears. Model parameters are never exactly known and system identification can be used in order to identify their values from measured input-output data. Uncertainty can be caused mainly by measurement noise and model error (e.g. unmodelled dynamics) and it always appears as an error between model and process outputs (identification error for a specific experiment). Well established identification techniques exist for linear models [6], [10], but not for nonlinear ones. Most of these techniques rely on cost function optimization where an error norm is considered. In particular, if the cost function is convex, local optimizers offer the solution, but in the nonconvex case, global optimization becomes necessary. In general, nonlinear models produce nonconvex optimization problems, and in this case, Evolutionary Algorithms [3], [5] offer a good solution. Furthermore, considering several norms of the identification error simultaneously (e.g. ∞ , 1 , 2 , F air, Huber, T urkey), the quality of the estimated models can be improved [8]. Thus, identification of nonlinear processes is stated as a multiobjective optimization problem. To solve this problem, a Multiobjective Evolutionary Algorithm ( −M OEA), based on -dominance concept [1], [7], has been developed. The algorithm converges to a well distributed sample of the Pareto Front. The paper is organized as follows. In section 2, robust identification problem is posed when different error norms are considered simultaneously. Section 3 des´ J. Mira and J.R. Alvarez (Eds.): IWINAC 2005, LNCS 3562, pp. 231–241, 2005. c Springer-Verlag Berlin Heidelberg 2005
232
J.M. Herrero et al.
cribes the −M OEA algorithms to find the Pareto optimal set. Finally, section 4 presents experimental results when −M OEA is applied to the identification of a thermal process represented by a nonlinear model with three unknown parameters.
2
Robust Identification Approach
In this work the following structure is assumed for the nonlinear model: ˙ x(t) = f (x(t), u(t), θ),
ˆ (t) = g(x(t), u(t), θ) y
(1)
where – – – – –
f (.), g(.) are the nonlinear functions of the model, θ ∈ D ⊂ Rq is the vector of unknown model parameters, x(t) ∈ Rn is the vector of model states, u(t) ∈ Rm is the vector of model inputs, ˆ (t) ∈ Rl is the vector of model outputs. y ˆ Let E(θ) = Y − Y(θ), where
– E(θ) is the identification error, – Y are the process output measurements , [y(0), y(T )...y(N T )], when the inputs U = [u(0), u(T )...u(N T )] are applied to the process, ˆ are the simulated1 outputs [ˆ ˆ (T )...ˆ y(0), y y(N T )], when the same inputs – Y: U are applied to the model. Denote E(θ)pi as the pi -norm of the identification error, with i ∈ A := [1, 2 . . . s]. When several norms of the identification error E(θ) are considered simultaneously, model identification can be posed as a multiobjective optimization problem min J(θ) θ∈D
(2)
where J(θ) = {J1 (θ), J2 (θ), . . . , Js (θ)} = {E(θ)p1 , E(θ)p2 , . . . , E(θ)ps } . Consequently, there is not a unique optimal model and a Pareto optimal set ΘP (solutions where no-one dominates others) must be found. Pareto dominance is defined as follows. A model θ1 , with cost function value J(θ1 ) dominates another model θ2 with cost function value J(θ2 ), denoted by J(θ1 ) ≺ J(θ2 ), iff ∀i ∈ A, Ji (θ1 ) ≤ Ji (θ2 ) ∧ ∃k ∈ A : Jk (θ1 ) < Jk (θ2 ) . 1
Model outputs are calculated integrating equations (1). T is the sample time and N is the number of measurements.
Nonlinear Robust Identification
233
Therefore the Pareto optimal set ΘP , is given by ˜ ≺ J(θ)} . ΘP = {θ ∈ D | θ˜ ∈ D : J(θ)
(3)
ΘP is unique and normally includes infinite models. Hence a set Θ∗P , with a finite number of elements from ΘP , will be obtained2 .
3
−M OEA
The variable MOEA ( −M OEA) is a multiobjective evolutionary algorithm based on -dominance concept. Consider the cost function space splitted up in a fixed number of boxes (for each dimension, n boxi cells of i wide)3 . That grid preserves diversity of J(Θ∗P ) since one box can be occupied by just one solution. This fact avoids the algorithm converging to just one point or area inside the cost function space (Fig. 1). The concept of -dominance is defined as follows. For a model θ, boxi (θ) is defined by4 Ji (θ) − Jimin boxi (θ) = ceil . i Let box(J(θ)) = {box1 (θ), . . . , boxs (θ)}. A model θ1 with cost function value J(θ1 ) -dominates the model θ2 with cost function value J(θ2 ), denoted by J(θ1 ) ≺ J(θ2 ), iff box(θ1 ) ≺ box(θ2 ) ∨ (box(θ1 ) = box(θ2 ) ∧ J(θ1 ) ≺ J(θ2 )) . Hence, a set Θ∗P is -Pareto iff Θ∗P ⊆ ΘP ∧ (box(θ1 ) = box(J(θ2 )), ∀θ1 , θ2 ∈ Θ∗P , θ1 = θ2 .
(4)
−M OEA algorithm obtains the -Pareto front J(Θ∗P ), a well-distributed sample of the Pareto front J(ΘP ). The algorithm, which adjusts the width dynamically, is composed of three populations (see Fig. 2). 1. Main population P (t) explores the searching space D during the algorithm iterations (t). Population size is N indP . 2. Archive A(t) stores the solution (Θ∗P ). Its size N indA can be variable and it will never be higher than N ind max A which depends on the number of the boxes (n box). 3. Auxiliary population G(t). Its size is N indG , which must be an even number. The pseudocode of the −M OEA algorithm is given by 2 3 4
Notice that ΘP∗ is not unique. i = (Jimax − Jimin )/n boxi . The function ceil(x) returns the smallest integer greater or equal than x.
234
J.M. Herrero et al. box
max
J2
J(qP)
Grey area is e-dominated by
J(qi)
e2
J(q*P)
qi
n_box1=10
are the points
n_box2=10 min
J2
min
max
e1
J1
J1
Fig. 1. -dominance concept. -Pareto Front J(Θ∗P ) in a bidimensional problem dominate update
e-dominate storeini
A(t)
create
P(t)
G(t)
store e-dominate
Fig. 2. −M OEA algorithm structure
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.
t := 0 A(t) := ∅ P (t) := ini random(D) eval(P (t)) A(t) := storeini (P (t), A(t)) while t < n iterations do G(t) := create(P (t), A(t)) eval(G(t)) A(t + 1) := store(G(t), A(t)) P (t + 1) := update(G(t), P (t)) t := t + 1 end while
Nonlinear Robust Identification
235
Main steps of the algorithm are detailed: Step 3. P (0) is initialized with N indP individuals (models) randomly selected from the searching space D. Step 4 and 8. Function eval calculates cost function value (2) for each individual in P (t) (step 4) and G(t) (step 8). Step 5. Function storeini checks individuals of P (t) that might be included in the archive A(t) as follows: 1. P(t) individuals non-dominated are detected, θN D . 2. Cost function space limits are calculated from J(θN D ). 3. Individuals of θN D are analyzed, one by one, and those not -dominated by individuals of A(t), will be included in A(t). Step 7. Each iteration, function create creates G(t) as follows: 1. Two individuals are randomly selected, one from P (t), θ1 , and another from A(t), θ2 . 2. θ1 and θ2 are crossed over by extended lineal recombination technique creating two new individuals θ1 and θ2 . 3. θ1 and θ2 are included in G(t). This procedure is repeated N indG /2 times until G(t) will be filled up. Step 9. Function store checks, one by one, which individuals of G(t) must be included in A(t) based on their location in the cost function space (see Fig. 3). Thus ∀θG ∈ G(t)
Z2
max
J2
Z3
Z1
min
J2
min
Z4
J1
Z3
max
J1
Fig. 3. Cost function space
236
J.M. Herrero et al.
1. If θG lies on the area Z1 and is not -dominated by any individual of A(t), it will be included in A(t)5 . Individuals of A(t) which are -dominated by θG will be eliminated. 2. If θG lies on the area Z2 then it is not included in the archive, since it is dominated by all individuals of A(t). 3. If θG lies on the area Z3, new cost function limits and i widths are recalculated. One by one the individuals of A(t) are again analyzed using the same procedure as storeini function. Finally, θG is included in A(t). 4. If θG lies on the area Z4, all individuals of A(t) are deleted since all of them are -dominated by θG . θG is included and cost function space limits are J(θG ). Step 10. Function update updates P (t) with individuals from G(t). Every individual θG from G(t) is randomly compared with an individual θP from P (t). If J(θG ) ≺ J(θ) then θG replaces θ, on the other hand, θ is maintained. Finally, individuals from A(t) compound the solution Θ∗P to the multiobjective minimization problem.
4
Robust Identification Example
Consider a scale furnace with a resistance placed inside. A fan continuously introduces air from outside (air circulation) while energy is supplied by an actuator controlled by voltage. Using a data acquisition system, resistance temperature and air temperature are measured when voltage is applied to the process. Fig. 4 shows the input signal applied and the output signal measured for an experiment of length N = 6000. These signals will be used for the robust identification problem. The dynamics of the resistance temperature can be modelled by 1 k1 u(t)2 − k2 ((x1 (t) − Ta (t)) + Of f set , 1000 x˙2 (t) = (1/k3 )(x1 (t) − x2 (t)) ,
x˙1 (t) =
(5)
yˆ(t) = x2 (t) , where – – – – – 5
x˙ 1 (t), x˙ 2 (t) are the model states, u(t) is the input voltage with rank 0 - 100 (%), yˆ(t) is the resistance temperature (o C) (model output), Ta (t) is the temperature inside furnace (o C), θ = [k1 , k2 , k3 ] are the model parameters to be identified,
When the individual is not -dominated and its box is occupied, that individual lying farthest from the box centre will be eliminated.
Nonlinear Robust Identification
237
Output y(ºC)
120 100 80 60 40 20
0
1000
2000
3000
4000
5000
6000
0
1000
2000
3000
4000
5000
6000
0
1000
2000
3000 Time(sec.)
4000
5000
6000
Input u(%)
80 60 40 20
Ta(ºC)
19 18 17 16
Fig. 4. Output process y(t), input process u(t) and disturbance Ta (t). T = 1sec. (sample time), N = 6000 (experiment length)
– Of f set is the correction to ensure zero steady-state error at a particular operating point 6 . To determine Θ∗P , two norms p1 = ∞ and p2 = 1 have been considered, to limit maximum and average identification error respectively. Therefore two cost functions have to be minimized to solve the problem stated at (2) J1 (θ) = E(θ)∞ , E(θ)1 . J2 (θ) = N The parameters of the −M OEA algorithm are set to:
(6) (7)
– The searching space D is k1 ∈ [0.05 . . . 0.12], k2 ∈ [3.0 . . . 8.0] and k3 ∈ [1.0 . . . 25.0]. – N indP = 8000, N indG = 8, and n iterations = 3000. – Cost funcion space parameter n box = [10, 10] (number of divisions for each dimension). At the end of the optimization process, the solution Θ∗P is A(3000), which contains six optimal models. Fig. 5 shows Θ∗P with its projections on the planes (k1 , k3 ), (k2 , k3 ) and (k2 , k1 ) respectively. 6
Before applying input (see Fig. 4) the process is in steady-state. Assuming that 1 (k1 u(0)2 − k2 (x1 (0) − Ta (0)), is not necessary x˙ 1 (0) = x˙ 2 (0) = 0 and Of f set = − 1000 to identify x1 (0) and x2 (0) since x1 (0) = x2 (0) = y(0).
238
J.M. Herrero et al. 8 25
7 6
15 k2
k3
20 10
5
5 8 6 4
0.06
k2
0.08
0.1
4
0.12
3
0.06
k1
20
20
15
15
0.1
0.12
k3
25
k3
25
0.08 k1
10
10
5
5 0.06
0.08 k1
0.1
0.12
3
4
5
6
7
8
k2
Fig. 5. The Pareto optimal set Θ∗P and its projections on (k1 , k3 ), (k2 , k3 ) and (k2 , k1 ) planes respectively
Fig. 6 shows associated Pareto front J(Θ∗P ). Notice that the algorithm has achieved an appropriate distribution of the Pareto front7 due to the gridding of the cost function space. The −M OEA algorithm succeeds in characterizing the front despite its short length and the differences among models. The Pareto optimal set Θ∗P is a sample of non-dominated solutions well distributed over the Pareto front J(ΘP ). Since the final objective is to obtain a unique model (θ∗ ), a way to select it has to be established8 . In this work, the proposal is based on the minimum distance to an ideal value cost function space. The ideal value is obtained from the extremes of the sampled Pareto front J(Θ∗P ): θ1 = arg min J1 (θ) = [0.078, 4.92, 4.77]; J1 (θ1 ) = 0.42, J2 (θ1 ) = 2.9 , θ∈D
θ∞ = arg min J2 (θ) = [0.082, 5.17, 9.44]; J1 (θ∞ ) = 0.480, J2 (θ∞ ) = 1.72 . θ∈D
Hence, Jideal = (J1min , J2min ) = (0.42, 1.72) . 7 8
If a better characterization of Pareto front is required, n box parameter should be increased. In multiobjective literature this is not a trivial issue and there is several methodologies regarding it.
Nonlinear Robust Identification
239
Archive 2.8
2.6
J2
2.4
2.2
J(θ *)
2
1.8
J(θ ideal) 0.425
0.43
0.435
0.44
0.445
0.45
0.455
0.46
0.465
0.47
0.475
J1
Fig. 6. The Pareto front J(Θ∗P ) 120 100
(a)
80 60 40 20
0
1000
2000
3000 Time(sec.)
4000
5000
6000
120 100
(b)
80 60 40 20 3800
3850
3900
3950
4000
4050 4100 Time(sec.)
4150
4200
4250
4300
Fig. 7. (a) For θ∗ model, real output y(t) (solid) and simulated output yˆ(t) (dash). (b) Notice differences between both outputs
240
J.M. Herrero et al.
The proposed solution θ∗ is the nearest model (in the cost function space) to the ideal model (see Fig. 6): θ∗ = [0.08, 5, 7.7]; J1 (θ∗ ) = 0.43, J2 (θ∗ ) = 2 . The estimated nonlinear model θ∗ can be validated by comparing its simulated output against measured output data (see Fig. 7).
5
Conclusions
A multiobjective evolutionary algorithm −M OEA, based on -dominance concept, has been developed for robust identification of nonlinear processes. Consequently, robust identification is posed as a multiobjective optimization problem and −M OEA estimates the nonlinear model set Θ∗P by assuming, simultaneously, the existence of several bounds in identification error. J(Θ∗P ) results in a well distributed sample of the optimal Pareto front J(ΘP ). The algorithm presents the following features: – Assuming parametric uncertainty, all kind of processes can be identified if their outputs can be calculated by model simulation. Differentiability respect to the unknown parameters is not necessary. – The algorithm dynamically adjusts the Pareto front precision without increasing the archive size N ind max A. – The algorithm adapts the extremes of the Pareto front with independence of the parameter n box. Consequently, good distribution and convergence to the Pareto front is achieved.
Acknowledgements This work has been partially financed by AGL 2002-04108-C02-01 and DPI 20013106-C02-02 MCYT (SPAIN)
References 1. Deb, K., Mohan, M., Mishra, S.: A fast multi-objective evolutionary algorithm for finding well-spread pareto-optimal solutions. Technical Report 2003002 KanGAL (2003) 2. Garulli, A., Kacewicz, B., Vicino, A., Zappa, G.: Error Bounds for Conditional Algorithms in Restricted Complexity Set Membership Identification. IEEE Transaction on Automatic Control, 45(1) (2000) 160-164 3. Goldberg, D.E.: Genetic Algorithms in search, optimization and machine learning. Addison-Wesley (1989) 4. Herrero, J.M., Blasco, X., Salcedo, J.V., Ramos, C.: Membership-Set Estimation with Genetic Algorithms in Nonlinear Models. In Proc. of the XV international Conference on Systems Science, (2004)
Nonlinear Robust Identification
241
5. Blasco, X., Herrero, J.M., Mart´ınez, M., Senent, J.: Nonlinear parametric model identification with Genetic Algorithms. Application to thermal process. Connectionist Models of Neurons, Learning Processes, and Artificial Intelligence (Springer). Lecture notes in computer science 2084 (2001) 6. Johansson, R.; System modeling identification. Prentice Hall (1993) 7. Laumanns, M., Thiele, L., Deb, K., Zitzler, E.: Combining convergence and diversity in evolutionary multi-objective optimization. Evolutionary computation, 10(3) (2002) 8. Ram, B., Gupta, H. Bandyopadhyay, P., Deb, K., Adimurthy, V.: Robust Identification of Aerospace Propulsion Parameters using Optimization Techniques based on Evolutionary Algorithms. Technical Report 2003005 KanGAL (2003) 9. Reinelt, W., Garulli, A., Ljung, L.: Comparing different approaches to model error modelling in robust identification. Automatica, 38(5) (2002) 787-803 10. Pronzalo, L., Walter, E.: Identification of parametric models from experimental data. Springer Verlang (1997)
1
Introduction
When modelling dynamic systems by first principles, the problem of parametric uncertainty always appears. Model parameters are never exactly known and system identification can be used in order to identify their values from measured input-output data. Uncertainty can be caused mainly by measurement noise and model error (e.g. unmodelled dynamics) and it always appears as an error between model and process outputs (identification error for a specific experiment). Well established identification techniques exist for linear models [6], [10], but not for nonlinear ones. Most of these techniques rely on cost function optimization where an error norm is considered. In particular, if the cost function is convex, local optimizers offer the solution, but in the nonconvex case, global optimization becomes necessary. In general, nonlinear models produce nonconvex optimization problems, and in this case, Evolutionary Algorithms [3], [5] offer a good solution. Furthermore, considering several norms of the identification error simultaneously (e.g. ∞ , 1 , 2 , F air, Huber, T urkey), the quality of the estimated models can be improved [8]. Thus, identification of nonlinear processes is stated as a multiobjective optimization problem. To solve this problem, a Multiobjective Evolutionary Algorithm ( −M OEA), based on -dominance concept [1], [7], has been developed. The algorithm converges to a well distributed sample of the Pareto Front. The paper is organized as follows. In section 2, robust identification problem is posed when different error norms are considered simultaneously. Section 3 des´ J. Mira and J.R. Alvarez (Eds.): IWINAC 2005, LNCS 3562, pp. 231–241, 2005. c Springer-Verlag Berlin Heidelberg 2005
232
J.M. Herrero et al.
cribes the −M OEA algorithms to find the Pareto optimal set. Finally, section 4 presents experimental results when −M OEA is applied to the identification of a thermal process represented by a nonlinear model with three unknown parameters.
2
Robust Identification Approach
In this work the following structure is assumed for the nonlinear model: ˙ x(t) = f (x(t), u(t), θ),
ˆ (t) = g(x(t), u(t), θ) y
(1)
where – – – – –
f (.), g(.) are the nonlinear functions of the model, θ ∈ D ⊂ Rq is the vector of unknown model parameters, x(t) ∈ Rn is the vector of model states, u(t) ∈ Rm is the vector of model inputs, ˆ (t) ∈ Rl is the vector of model outputs. y ˆ Let E(θ) = Y − Y(θ), where
– E(θ) is the identification error, – Y are the process output measurements , [y(0), y(T )...y(N T )], when the inputs U = [u(0), u(T )...u(N T )] are applied to the process, ˆ are the simulated1 outputs [ˆ ˆ (T )...ˆ y(0), y y(N T )], when the same inputs – Y: U are applied to the model. Denote E(θ)pi as the pi -norm of the identification error, with i ∈ A := [1, 2 . . . s]. When several norms of the identification error E(θ) are considered simultaneously, model identification can be posed as a multiobjective optimization problem min J(θ) θ∈D
(2)
where J(θ) = {J1 (θ), J2 (θ), . . . , Js (θ)} = {E(θ)p1 , E(θ)p2 , . . . , E(θ)ps } . Consequently, there is not a unique optimal model and a Pareto optimal set ΘP (solutions where no-one dominates others) must be found. Pareto dominance is defined as follows. A model θ1 , with cost function value J(θ1 ) dominates another model θ2 with cost function value J(θ2 ), denoted by J(θ1 ) ≺ J(θ2 ), iff ∀i ∈ A, Ji (θ1 ) ≤ Ji (θ2 ) ∧ ∃k ∈ A : Jk (θ1 ) < Jk (θ2 ) . 1
Model outputs are calculated integrating equations (1). T is the sample time and N is the number of measurements.
Nonlinear Robust Identification
233
Therefore the Pareto optimal set ΘP , is given by ˜ ≺ J(θ)} . ΘP = {θ ∈ D | θ˜ ∈ D : J(θ)
(3)
ΘP is unique and normally includes infinite models. Hence a set Θ∗P , with a finite number of elements from ΘP , will be obtained2 .
3
−M OEA
The variable MOEA ( −M OEA) is a multiobjective evolutionary algorithm based on -dominance concept. Consider the cost function space splitted up in a fixed number of boxes (for each dimension, n boxi cells of i wide)3 . That grid preserves diversity of J(Θ∗P ) since one box can be occupied by just one solution. This fact avoids the algorithm converging to just one point or area inside the cost function space (Fig. 1). The concept of -dominance is defined as follows. For a model θ, boxi (θ) is defined by4 Ji (θ) − Jimin boxi (θ) = ceil . i Let box(J(θ)) = {box1 (θ), . . . , boxs (θ)}. A model θ1 with cost function value J(θ1 ) -dominates the model θ2 with cost function value J(θ2 ), denoted by J(θ1 ) ≺ J(θ2 ), iff box(θ1 ) ≺ box(θ2 ) ∨ (box(θ1 ) = box(θ2 ) ∧ J(θ1 ) ≺ J(θ2 )) . Hence, a set Θ∗P is -Pareto iff Θ∗P ⊆ ΘP ∧ (box(θ1 ) = box(J(θ2 )), ∀θ1 , θ2 ∈ Θ∗P , θ1 = θ2 .
(4)
−M OEA algorithm obtains the -Pareto front J(Θ∗P ), a well-distributed sample of the Pareto front J(ΘP ). The algorithm, which adjusts the width dynamically, is composed of three populations (see Fig. 2). 1. Main population P (t) explores the searching space D during the algorithm iterations (t). Population size is N indP . 2. Archive A(t) stores the solution (Θ∗P ). Its size N indA can be variable and it will never be higher than N ind max A which depends on the number of the boxes (n box). 3. Auxiliary population G(t). Its size is N indG , which must be an even number. The pseudocode of the −M OEA algorithm is given by 2 3 4
Notice that ΘP∗ is not unique. i = (Jimax − Jimin )/n boxi . The function ceil(x) returns the smallest integer greater or equal than x.
234
J.M. Herrero et al. box
max
J2
J(qP)
Grey area is e-dominated by
J(qi)
e2
J(q*P)
qi
n_box1=10
are the points
n_box2=10 min
J2
min
max
e1
J1
J1
Fig. 1. -dominance concept. -Pareto Front J(Θ∗P ) in a bidimensional problem dominate update
e-dominate storeini
A(t)
create
P(t)
G(t)
store e-dominate
Fig. 2. −M OEA algorithm structure
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.
t := 0 A(t) := ∅ P (t) := ini random(D) eval(P (t)) A(t) := storeini (P (t), A(t)) while t < n iterations do G(t) := create(P (t), A(t)) eval(G(t)) A(t + 1) := store(G(t), A(t)) P (t + 1) := update(G(t), P (t)) t := t + 1 end while
Nonlinear Robust Identification
235
Main steps of the algorithm are detailed: Step 3. P (0) is initialized with N indP individuals (models) randomly selected from the searching space D. Step 4 and 8. Function eval calculates cost function value (2) for each individual in P (t) (step 4) and G(t) (step 8). Step 5. Function storeini checks individuals of P (t) that might be included in the archive A(t) as follows: 1. P(t) individuals non-dominated are detected, θN D . 2. Cost function space limits are calculated from J(θN D ). 3. Individuals of θN D are analyzed, one by one, and those not -dominated by individuals of A(t), will be included in A(t). Step 7. Each iteration, function create creates G(t) as follows: 1. Two individuals are randomly selected, one from P (t), θ1 , and another from A(t), θ2 . 2. θ1 and θ2 are crossed over by extended lineal recombination technique creating two new individuals θ1 and θ2 . 3. θ1 and θ2 are included in G(t). This procedure is repeated N indG /2 times until G(t) will be filled up. Step 9. Function store checks, one by one, which individuals of G(t) must be included in A(t) based on their location in the cost function space (see Fig. 3). Thus ∀θG ∈ G(t)
Z2
max
J2
Z3
Z1
min
J2
min
Z4
J1
Z3
max
J1
Fig. 3. Cost function space
236
J.M. Herrero et al.
1. If θG lies on the area Z1 and is not -dominated by any individual of A(t), it will be included in A(t)5 . Individuals of A(t) which are -dominated by θG will be eliminated. 2. If θG lies on the area Z2 then it is not included in the archive, since it is dominated by all individuals of A(t). 3. If θG lies on the area Z3, new cost function limits and i widths are recalculated. One by one the individuals of A(t) are again analyzed using the same procedure as storeini function. Finally, θG is included in A(t). 4. If θG lies on the area Z4, all individuals of A(t) are deleted since all of them are -dominated by θG . θG is included and cost function space limits are J(θG ). Step 10. Function update updates P (t) with individuals from G(t). Every individual θG from G(t) is randomly compared with an individual θP from P (t). If J(θG ) ≺ J(θ) then θG replaces θ, on the other hand, θ is maintained. Finally, individuals from A(t) compound the solution Θ∗P to the multiobjective minimization problem.
4
Robust Identification Example
Consider a scale furnace with a resistance placed inside. A fan continuously introduces air from outside (air circulation) while energy is supplied by an actuator controlled by voltage. Using a data acquisition system, resistance temperature and air temperature are measured when voltage is applied to the process. Fig. 4 shows the input signal applied and the output signal measured for an experiment of length N = 6000. These signals will be used for the robust identification problem. The dynamics of the resistance temperature can be modelled by 1 k1 u(t)2 − k2 ((x1 (t) − Ta (t)) + Of f set , 1000 x˙2 (t) = (1/k3 )(x1 (t) − x2 (t)) ,
x˙1 (t) =
(5)
yˆ(t) = x2 (t) , where – – – – – 5
x˙ 1 (t), x˙ 2 (t) are the model states, u(t) is the input voltage with rank 0 - 100 (%), yˆ(t) is the resistance temperature (o C) (model output), Ta (t) is the temperature inside furnace (o C), θ = [k1 , k2 , k3 ] are the model parameters to be identified,
When the individual is not -dominated and its box is occupied, that individual lying farthest from the box centre will be eliminated.
Nonlinear Robust Identification
237
Output y(ºC)
120 100 80 60 40 20
0
1000
2000
3000
4000
5000
6000
0
1000
2000
3000
4000
5000
6000
0
1000
2000
3000 Time(sec.)
4000
5000
6000
Input u(%)
80 60 40 20
Ta(ºC)
19 18 17 16
Fig. 4. Output process y(t), input process u(t) and disturbance Ta (t). T = 1sec. (sample time), N = 6000 (experiment length)
– Of f set is the correction to ensure zero steady-state error at a particular operating point 6 . To determine Θ∗P , two norms p1 = ∞ and p2 = 1 have been considered, to limit maximum and average identification error respectively. Therefore two cost functions have to be minimized to solve the problem stated at (2) J1 (θ) = E(θ)∞ , E(θ)1 . J2 (θ) = N The parameters of the −M OEA algorithm are set to:
(6) (7)
– The searching space D is k1 ∈ [0.05 . . . 0.12], k2 ∈ [3.0 . . . 8.0] and k3 ∈ [1.0 . . . 25.0]. – N indP = 8000, N indG = 8, and n iterations = 3000. – Cost funcion space parameter n box = [10, 10] (number of divisions for each dimension). At the end of the optimization process, the solution Θ∗P is A(3000), which contains six optimal models. Fig. 5 shows Θ∗P with its projections on the planes (k1 , k3 ), (k2 , k3 ) and (k2 , k1 ) respectively. 6
Before applying input (see Fig. 4) the process is in steady-state. Assuming that 1 (k1 u(0)2 − k2 (x1 (0) − Ta (0)), is not necessary x˙ 1 (0) = x˙ 2 (0) = 0 and Of f set = − 1000 to identify x1 (0) and x2 (0) since x1 (0) = x2 (0) = y(0).
238
J.M. Herrero et al. 8 25
7 6
15 k2
k3
20 10
5
5 8 6 4
0.06
k2
0.08
0.1
4
0.12
3
0.06
k1
20
20
15
15
0.1
0.12
k3
25
k3
25
0.08 k1
10
10
5
5 0.06
0.08 k1
0.1
0.12
3
4
5
6
7
8
k2
Fig. 5. The Pareto optimal set Θ∗P and its projections on (k1 , k3 ), (k2 , k3 ) and (k2 , k1 ) planes respectively
Fig. 6 shows associated Pareto front J(Θ∗P ). Notice that the algorithm has achieved an appropriate distribution of the Pareto front7 due to the gridding of the cost function space. The −M OEA algorithm succeeds in characterizing the front despite its short length and the differences among models. The Pareto optimal set Θ∗P is a sample of non-dominated solutions well distributed over the Pareto front J(ΘP ). Since the final objective is to obtain a unique model (θ∗ ), a way to select it has to be established8 . In this work, the proposal is based on the minimum distance to an ideal value cost function space. The ideal value is obtained from the extremes of the sampled Pareto front J(Θ∗P ): θ1 = arg min J1 (θ) = [0.078, 4.92, 4.77]; J1 (θ1 ) = 0.42, J2 (θ1 ) = 2.9 , θ∈D
θ∞ = arg min J2 (θ) = [0.082, 5.17, 9.44]; J1 (θ∞ ) = 0.480, J2 (θ∞ ) = 1.72 . θ∈D
Hence, Jideal = (J1min , J2min ) = (0.42, 1.72) . 7 8
If a better characterization of Pareto front is required, n box parameter should be increased. In multiobjective literature this is not a trivial issue and there is several methodologies regarding it.
Nonlinear Robust Identification
239
Archive 2.8
2.6
J2
2.4
2.2
J(θ *)
2
1.8
J(θ ideal) 0.425
0.43
0.435
0.44
0.445
0.45
0.455
0.46
0.465
0.47
0.475
J1
Fig. 6. The Pareto front J(Θ∗P ) 120 100
(a)
80 60 40 20
0
1000
2000
3000 Time(sec.)
4000
5000
6000
120 100
(b)
80 60 40 20 3800
3850
3900
3950
4000
4050 4100 Time(sec.)
4150
4200
4250
4300
Fig. 7. (a) For θ∗ model, real output y(t) (solid) and simulated output yˆ(t) (dash). (b) Notice differences between both outputs
240
J.M. Herrero et al.
The proposed solution θ∗ is the nearest model (in the cost function space) to the ideal model (see Fig. 6): θ∗ = [0.08, 5, 7.7]; J1 (θ∗ ) = 0.43, J2 (θ∗ ) = 2 . The estimated nonlinear model θ∗ can be validated by comparing its simulated output against measured output data (see Fig. 7).
5
Conclusions
A multiobjective evolutionary algorithm −M OEA, based on -dominance concept, has been developed for robust identification of nonlinear processes. Consequently, robust identification is posed as a multiobjective optimization problem and −M OEA estimates the nonlinear model set Θ∗P by assuming, simultaneously, the existence of several bounds in identification error. J(Θ∗P ) results in a well distributed sample of the optimal Pareto front J(ΘP ). The algorithm presents the following features: – Assuming parametric uncertainty, all kind of processes can be identified if their outputs can be calculated by model simulation. Differentiability respect to the unknown parameters is not necessary. – The algorithm dynamically adjusts the Pareto front precision without increasing the archive size N ind max A. – The algorithm adapts the extremes of the Pareto front with independence of the parameter n box. Consequently, good distribution and convergence to the Pareto front is achieved.
Acknowledgements This work has been partially financed by AGL 2002-04108-C02-01 and DPI 20013106-C02-02 MCYT (SPAIN)
References 1. Deb, K., Mohan, M., Mishra, S.: A fast multi-objective evolutionary algorithm for finding well-spread pareto-optimal solutions. Technical Report 2003002 KanGAL (2003) 2. Garulli, A., Kacewicz, B., Vicino, A., Zappa, G.: Error Bounds for Conditional Algorithms in Restricted Complexity Set Membership Identification. IEEE Transaction on Automatic Control, 45(1) (2000) 160-164 3. Goldberg, D.E.: Genetic Algorithms in search, optimization and machine learning. Addison-Wesley (1989) 4. Herrero, J.M., Blasco, X., Salcedo, J.V., Ramos, C.: Membership-Set Estimation with Genetic Algorithms in Nonlinear Models. In Proc. of the XV international Conference on Systems Science, (2004)
Nonlinear Robust Identification
241
5. Blasco, X., Herrero, J.M., Mart´ınez, M., Senent, J.: Nonlinear parametric model identification with Genetic Algorithms. Application to thermal process. Connectionist Models of Neurons, Learning Processes, and Artificial Intelligence (Springer). Lecture notes in computer science 2084 (2001) 6. Johansson, R.; System modeling identification. Prentice Hall (1993) 7. Laumanns, M., Thiele, L., Deb, K., Zitzler, E.: Combining convergence and diversity in evolutionary multi-objective optimization. Evolutionary computation, 10(3) (2002) 8. Ram, B., Gupta, H. Bandyopadhyay, P., Deb, K., Adimurthy, V.: Robust Identification of Aerospace Propulsion Parameters using Optimization Techniques based on Evolutionary Algorithms. Technical Report 2003005 KanGAL (2003) 9. Reinelt, W., Garulli, A., Ljung, L.: Comparing different approaches to model error modelling in robust identification. Automatica, 38(5) (2002) 787-803 10. Pronzalo, L., Walter, E.: Identification of parametric models from experimental data. Springer Verlang (1997)
