Interactive Evolutionary Computation in Identification of Dynamical Systems Janos Abonyi, Janos Madar, Lajos Nagy, and Ferenc Szeifert University of Veszprem, Department of Process Engineering, Veszprem, P.O. Box 158, H-8201, Hungary [email protected] Summary. In practical system identification it is often desirable to simultaneously handle several objectives and constraints. In some cases, these objectives and constraints are often non-commensurable and the objective functions are explicitly/mathematically not available. In this paper, Interactive Evolutionary Computation (IEC) is used to effectively handle these identification problems. IEC is an optimization method that adopts evolutionary computation (EC) among system optimization based on subjective human evaluation. The proposed approach has been implemented in MATLAB (EAsy-IEC Toolbox) and applied to the identification of a pilot batch reactor. The results show that IEC is an efficient and comfortable method to incorporate a priori knowledge of the user into a user-guided optimization and identification problems. The developed EASy-IEC Toolbox can be downloaded from the website of the authors: http://www.fmt.vein.hu/softcomp/EAsy.

Key words: System identification, Interactive Evolutionary Computation

1 Introduction The main objective of system identification is to identify a model with good prediction capabilities in the sense that it is able to accurately predict the system’s response to a given class of excitations. Hence, a common identification objective is to minimize the mismatch between model prediction and observed data. However, it is often desirable to introduce additional objectives and constraints into the identification problem. There are several reasons for this [7]: In most practical cases the data sequence used for identification may be incomplete or uncertain, the model structure is not necessarily known and can be overparameterized, there may be some properties that we may want the model to have. Furthermore, there is often additional information available, and it makes sense to use this information to improve the accuracy and validity of the model [1]. Multi-objective system identification tools allow the incorporation of these objectives into the identification procedure. E.g.,

2

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

single constraints or penalties on the parameter space have been suggested, for example to ensure stability [24, 1], fulfillment of balance equations and steady-state data [23], and explicit belief about parameter values [14, 10]. These multi-objective identification problems involve the simultaneous consideration of multiple performance criteria. These objectives are often noncommensurable and are frequently in conflict with one another [7]. Trade-offs exist between some objectives where advancement in one objective will cause deterioration in another. In most cases the cost function is defined prior to the optimization procedure; this requires in-depth information concerning the various trade-offs and valuation of each individual objective. However, in some cases, the objective functions are explicitly/mathematically not available and the identification problem can be decomposed only for large-scale problems. This paper proposes the application of Interactive Evolutionary Computation (IEC) [21] to handle these identification problems, and its application to the identification of tendency models of batch reactors. In the proposed framework, human users evaluate the performances of the individuals (set of solutions) by selecting the best individuals. The human user does not (only) evaluate the mean square fitting error, or time constant of the model but he or she can simultaneously analyze these numerical values with the plotted dynamic and/or steady-state behavior of the model, e.g. step-responses, gains, etc. The proposed approach has been implemented in MATLAB (EAsy-IEC Toolbox) and applied to the identification of a pilot batch reactor. The remaining part of the paper is organized as follow. In Section 2 the algorithm of Evolutionary Strategy is presented, which is tailored taking into account the specialities of IEC. Section 3 shows an application example when the proposed approach is applied to identification of a tendency model of a heating-cooling jacket. Finally, some conclusions are drawn in Section 4.

2 Description of the Algorithm 2.1 Introduction to Evolutionary Algorithms Evolutionary Algorithms (EAs) [18] are stochastic optimization algorithms that mimic the process of natural selection. Evolutionary Algorithms handle several candidate solutions simultaneously. Thus, unlike conventional algorithms, EAs start form a pool of points (candidate solutions), usually referred to as individuals. Every individual represents a point in the search space. Every individual has a fitness value, which is calculated by the evaluation function. The fitness value determines the successfulness of the given individual appropriate to the objectives. The next generation of individuals is generated from the current population using genetic operators: selection, mutation and crossover. Likes in the natural world, the more successful individuals, which have bigger fitness values, have higher probability to transmit their genes into

Interactive Evolutionary Computation

3

the new generation. Because the genetic operators are stochastic, EAs are stochastic optimization methods. In the initialization step, EAs generate individuals randomly. In every evolutionary-loop, all of the individuals are evaluated and their fitness values are calculated based on an objective function. After that the selection operator selects randomly the best individuals to constitute the next generation. The crossover (recombination) operator generates two new individuals by taking two selected individuals and recombining them. The mutation operator changes randomly the individuals. Finally the old individuals are replaced by the new individuals. Generally, EAs are automated algorithms using a given objective function to calculate the fitness values of the individuals and the selection operator performs automatically the selection based on these calculated fitness values (for example using roulette wheel selection). But sometimes instead of a given objective function, a human user can make the valuation of the individuals, which results in an interactive optimization. Interactive Evolutionary Computation (IEC) is an optimization method that adopts evolutionary computation among system optimization based on subjective human evaluation [21]. In contrast to the automated evolutionary optimization algorithms, IEC cannot use many individuals and searching generations because of resulting human fatigue. Hence, there is a need to develop and tune an algorithm which effectively search with a few individuals within a few searching generations. For this purpose, this paper presents an algorithm based on Evolutionary Strategy (ES). In the following the IEC and the ES will be described and then the structure of the developed IEC-Tune Toolbox will be shown. 2.2 Interactive Evolutionary Computation Most optimization techniques which work by improving a single solution step by step are not suited for interactive optimization approach. For example, in simple hill-climbing a user suggestion can either be accepted, which means the loss of results gained so far, or denied, which may lead the algorithm a wrong way. Evolutionary Algorithms and other population based optimization procedures are better suited for interactive optimization, because a suggestion can be put into the population to see whether it is worth further consideration or not. A good suggestion can be expected to survive and lead the optimization process towards a good solution, while a bad suggestion will be eliminated quite soon. This approach become well known as Interactive Evolutionary Computation (IEC) [20, 19]. Simply stated, IEC is a technique from the class of evolutionary algorithms (EAs), in which the fitness function is replaced by a human user. As in interactive evolution, the user selects one or more individual(s) which survive(s) and reproduce(s) (with variation) to constitute a new generation. Hence IEC uses two different spaces for its search: the human user evaluates the output

4

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

of the target system according to the distance between the target goal and the system output in psychological space [20], on the other hand, the EA searches in the parameter space. These techniques were applied already in computer graphics [9], animation, creating forms, textures, and motion [2]. Potential applications of interactive evolution include artificial life design, e.g., development of components of biological nature [17, 6] and engineering construction design [13]. Although, IEC has been applied in knowledge acquisition and data mining, according to our best knowledge, the applicability of this approach in system identification has not been studied yet. However, in recent years, the application of Evolutionary Computing to black-box and grey-box model identification has received considerable interest since the seminal paper by Kristinsson and Dumont [8]. Subsequent evolutionary system identification applications are applied to build dominant structural identification with local parametric tuning without the need of a differentiable performance index in the presence of noisy data. As the overview of Fleming [5] illustrates evolutionary techniques provide an excellent fitting performance and is capable of accommodating multiple objectives such as to examine the relationships between model complexity and fitting accuracy during the model building process. These approaches can be extended to the identification of gray-box models, to use of a priori knowledge [22]. 2.3 Evolution Strategy In contrast to automated evolutionary optimization algorithms, the IEC cannot use many individuals and searching generations because the human fatigue. Hence, there is a need to develop an algorithm which can effectively search with a few individuals within a few searching generations. For this purpose, this paper presents an algorithm based on Evolutionary Strategy. Evolution Strategy (ES) was developed by [15], with selection, mutation, and a population of size one. Schwefel [16] introduced recombination and populations with more than one individual, and provided a nice comparison of ES with more traditional optimization techniques. Evolutionary Strategy are typically applied to real-valued parameter optimization problems. The characteristic feature of ES lies in the self-adaptation of the standard deviation of the Gaussian distribution used in the mutation [3]. The basic idea is to add these adaptation parameters to the genotype and have them undergo evolution themselves. The main elements of the algorithm used in this paper are the following: Representation An ES-individual aj = (xj , σ j ) consists of two components: the object variables xj = [xj,1 , . . . , xj,n ] and the strategy variables σ j = [σj,1 , . . . , σj,n ]. The strategy parameters allow the algorithm to adapt to the topology of the objective function by controlling the mutation step size.

Interactive Evolutionary Computation

5

Mutation As small changes occur frequently but large ones only rarely in the nature, the mutation operator adds normal distributed random numbers zj,i ∼ N (0, σj,i ) to the individuals: xj,i = xj,i + zj,i . (1) Before the object variables are changed, the standard deviations are mutated using a multiplicative normally distributed process: (t)

(t−1)

σj,i = σj,i

exp(τ 0 N (0, 1) + τ Ni (0, 1)),

(2)

with exp (τ 0 N (0, 1)) as a global factor which allows an overall change of the mutability and exp (τ Ni (0, 1)) allowing for individual changes of the mean step sizes σj,i . The τ 0 and τ parameters can be interpreted in the sense of global learning rates. Schwefel suggests to set them as [16]: 1 1 τ0 = √ , τ = p √ . 2n 2 n

(3)

Recombination Recombination in ESs can be either sexual, where only two parents are involved in the creation of an offspring, or global, where up to the whole population contributes to a new offspring. Sexual recombination of just two individuals is often called local while the contribution of all individuals is called global recombination. Traditional recombination operators are discrete recombination, intermediate recombination, and geometric recombination, all existing in a sexual and global form. When F and M denote two randomly selected individuals from the µ parent population, the following operators can be defined: xF,i no recombination xF,i or xM,i discrete 0 xi = (4) (x + x ) /2 intermediate F,i M,i Pµ global avarage k=1 xK,i /µ σF,i no recombination σ or x discrete F,i M,i + σ ) /2 intermediate σi0 = (σ . (5) M,i pF,i (σ σ ) geometric Pµ F,i M,i global avarage k=1 σK,i /µ Throughout this work discrete recombination of the object variables and global intermediate combination of the strategy parameters were used.

6

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

Selection and Replacement In the applied (µ + λ) Evolution Strategy, at a given generation (t), there are µ parents, and λ offspring are generated by recombination and mutation. Then, the µ+λ members of then next population are sorted according to their objective function values given by the user, and he or she selects the best µ to become the parents of the next generation. 2.4 IEC Integration to MATLAB IEC systems seek to interface EA to human preference in order to create systems capable of evolving artifacts that require a human expertise that has not yet succumb to computation. Hence, the interfacing of human ability with machine computation requires resolving difficult issues [12]. Further, progress in the design of interactive evolutionary systems allows a glimpse into how very human abilities such as intuition, projection, and holistic perception interplay with the mechanics of machine computation. This section reports on one such interactive evolutionary system that seeks to combine human perception with the previously presented evolutionary strategy to evolve the parameters of MATLAB and Simulink models (see Figure 1). The EAsy-IEC Tune Toolbox is downloadable from the website of the the author: www.fmt.vein.hu/softcomp/EAsy, is designed to be applicable for different types of optimization problems (e.g., system identification, controller tuning, data mining). The number of displayed individuals, µ + λ is set to 9, which can be displayed spatially. The number of searching generations is limited to 20 generations at the most, due to the fatigue of human operators. In the proposed framework human users evaluate the performances of the individuals (set of solutions). Not the genotypes and phenotypes are applied directly, but rather the output of the target system realized by the individuals. For example, the user does not (only) evaluate the mean square fitting error, or time constant of the model but he or she can simultaneously analyze these numerical values with the plotted dynamic and/or steady state behavior of the model, e.g.

Fig. 1. Scheme of the Proposed IEC System Integrated to MATLAB

Interactive Evolutionary Computation

7

step-responses, gains, etc. An example of such evaluation figure is given in Figure 4. Due to small number of individuals, the human evaluation was restricted to selecting the individuals which will survive and constitute the next generation.

3 Application Example 3.1 Process Description In recent years our department focuses to the development of control strategies of multipurpose production plants. Hence, a prototype of pharmaceutical process systems was designed and installed in our laboratory. The central element of this process unit is a 50 liter stirred reactor with heating-cooling jacket (see fig. 2).

Fig. 2. The stirred tank reactor with heating-cooling jacket

This paper focuses to the identification of the tendency model of the heating-cooling jacket. Through the jacket the direct heating-cooling system allows cooling with chilled water or heating with steam in steam or hot water mode. In the hot water heating mode the water is circulated by a pump, the steam is introduced through a mixer, while the excess water is removed through an overflow. 3.2 Tendency Model of the Heating-Cooling Jacket The jacket and the circulating liquid can be described using the common lumped parameter enthalpy or heat balances given in the chemical engineering literature. In the model, zero-volume distributors and mixers and are applied and the overflow as well as the feed of the steam and the fresh cooling water

8

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

are taken into account. The obtained simplified first principle model can be regarded as the tendency model [4] of the most important phenomena and its scheme is given on Figure 3.

Fig. 3. Tendency model of the jacket of the stirred ractor

V 1 and V 2 are the valve positions (0%-100%), T H is the temperature of the cooling water coming from the environment, T F is the equivalent steam temperature calculated from the boiling temperature, the latent heat and specific heat. The valve characteristics are given in form of second order polynomials (V ALV E1 and V ALV E2). The two first order transfer functions 1 1 ( T 1s+1 and T 2s+1 ) are obtained from the lumped parameter model of the jacket and the model of the thermometer at the jacket inlet. The K1, K2, and K3 are gains and the DT is the dead-time. The identification problem is identifying these parameters of the model (gains, time constants, gains, dead-time and parameters of second order polynomials). For this purpose open-loop experiments were conducted on the whole operating range of the composite manipulated variable u (because of the split-range control: V 1 and V 2 are not opened at the same time, u = 100% means that V 1 is opened fully, u = 0% means that V 2 is opened fully). 3.3 Results Based on the collected input-output data, the identification problem can be formulated as an optimization problem and can be solved using common nonlinear optimization methods (e.g. SQP). To formulate an optimization problem we need a cost-function which is minimized under the optimization procedure. Usually for this purpose the mean square model error is used, n P 2 M SE = (y(k) − y˜(k)) /n, where y is the observed output, y˜ is the prek=1

dicted output of the model, n is the number of observation. Beside the MSE, other objectives should be considered. E.g. the steadystate prediction error and the characteristics of the valves. But it is difficult to formulate these objectives and it is difficult to select the weights of the

Interactive Evolutionary Computation

9

objectives in the cost-function. Hence it is problematical to introduce these additional objectives into the identification problem. The IEC handles this problem, because the human user can consider and balance these objectives. For comparison, beside the proposed IEC optimization scheme, the direct minimization of the MSE prediction error is considered. For this purpose, Sequential Quadratic Programming (SQP) was used. The SQP method is a standard general purpose algorithm for solving NLP optimization problems. It is a gradient based method which calculates the search direction vector from the quadratic approximation of the objective function around the actual local point [11]. The SQP based identification method minimizes automatically the cost-function subject to the parameters of the model. The cost-function was evaluated by simulating the model and calculating the MSE. In contrast to SQP based method, the user evaluates the model in the IEC based identification method. Firstly, the identified model was simulated and then the human user selected the best solutions based on the dynamic and the near steady-state behavior of the system, and the identified characteristic of the valves. Beside these plots, some numerical values were also given to the user to glimpse. Such an evaluation display is depicted in Figure 4. 100

100

100

100

100

100

100

100

50

50

50

50

50

50

50

50

50

Temp. [C]

100

0 0

Temp. [C]

100

0 0 0 0 0 0 0 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] 100 100 100 100 100 100 100 100

50

0

F1%, F2%

100

50

50

50

50

50

50

50

0 0 0 0 0 0 0 0 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] 100 100 100 100 100 100 100 100

50

0

50

50

50

50

50

50

50

50

50

0 0 0 0 0 0 0 0 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% ISE: 15.1 K1: 0.055 K2: 0.357 K3: 0.025 T0: 0.559 Th: 0.422

ISE: 9.3 K1: 0.055 K2: 0.382 K3: 0.024 T0: 0.560 Th: 0.424

ISE: 38.8 K1: 0.052 K2: 0.382 K3: 0.024 T0: 0.536 Th: 0.424

ISE: 101.6 K1: 0.038 K2: 0.375 K3: 0.016 T0: 0.392 Th: 0.417

ISE: 339.2 K1: 0.065 K2: 0.387 K3: 0.030 T0: 0.687 Th: 0.428

ISE: 42.4 K1: 0.043 K2: 0.378 K3: 0.019 T0: 0.438 Th: 0.420

ISE: 213.2 K1: 0.024 K2: 0.344 K3: 0.010 T0: 0.229 Th: 0.411

ISE: 409.8 K1: 0.013 K2: 0.340 K3: 0.010 T0: 0.143 Th: 0.404

ISE: 3948.8 K1: 0.071 K2: 0.390 K3: 0.033 T0: 0.719 Th: 0.433

Fig. 4. MATLAB Evaluation figure in the IEC-Tune toolbox

From an arbitrary initial value, the application of IEC resulted in 5.0 MSE prediction error after 14 generations. This corresponds to 146 evaluations of the model. If the same numbers of function calls are applied in SQP, the obtain model results in 26.0 MSE. The drawback of IEC is that is slower than the SQP method due to the time-requirement of the human evaluation. Hence while the IEC procedure consumed about 8.5 minutes, the SQP consumed only 2.5 minutes. (See the Table 1. and Figures 5, 6.)

10

J. Abonyi, J. Madar, L. Nagy and F. Szeifert Table 1. Comparative results Evaulation number 146 149 506

Time (sec) 508 150 501

MSE 5.0 26.0 1.5

Temp. [C]

100

50

0

0

10

20

30 Time [sec]

40

50

60

Temp. [C]

100

50

0

80

100

120

140 Time [sec]

40

50 V1%, V2%

160

180

200

F1%, F2%

100

50

0

0

10

20

30

60

70

80

90

100

Fig. 5. Results using IEC after 146 model evaluations

Temp. [C]

100

50

0

0

10

20

30 Time [sec]

40

50

60

Temp. [C]

100

50

0

80

100

120

140 Time [sec]

40

50 V %, V %

160

180

200

100 F1%, F2%

Method IEC SQP SQP

50

0

0

10

20

30

1

60

70

80

90

100

2

Fig. 6. Results using SQP after 149 model evaluations

Interactive Evolutionary Computation

11

Our experiences and these results show that the proposed IEC approach is extremely useful to find good initial conditions to SQP that is sensitive to stacking in local minima.

4 Conclusions This paper demonstrated how interactive evolution algorithms can be applied to solve multi-objective identification problems. A MATLAB toolbox has been developed to provide a human-machine interface of the IEC. With the use of this tool, the user can easily select the individuals which then are reproduced and mutated to constitute the next generation. The proposed approach applied to the identification of a pilot batch reactor. The results show that IEC is an efficient and comfortable method to incorporate the priori knowledge and preferences of the user into the identification problem. Our experiences show, IEC is extremely useful to find good initial conditions to other nonlinear optimization algorithms, like SQP.

Acknowledgement The authors would like to acknowledge the support of the Cooperative Research Center (VIKKK) (KKK-I-7), the Hungarian Ministry of Education (FKFP-0073/2001), and the Hungarian Science Foundation (OTKA TO37600). Janos Abonyi is grateful for the financial support of the Janos Bolyai Research Fellowship of the Hungarian Academy of Science.

References 1. J. Abonyi, R. Babuska, H.B. Verbruggen, and F. Szeifert. Constraint parameter estimation in fuzzy modelling. In FUZZ-IEEE’99 Conference, pages 951–956. Seoul, Korea, 1999. 2. K. Aoki, H. Takagi, and N. Fujimura. Interactive ga-based design support system for lighting design in computer graphics. In Int’l Conf. on Soft Computing (IIZUKA’96), Iizuka, Fukuoka, Japan, pages 533–536. World Scientic, 1996. 3. A.E. Eiben and M. Schoenauer. Evolutionary computing. Information Processing Letters, 82, 2002. 4. C. Filippi-Bossy, C. Georgakis, J. Bordet, J. Villermaux, and S. MarchalBrassely. Batch reactor optimization by use of tendency models. Comp. & Chem. Eng, 13:35–47, 1989. 5. P.J. Fleming and R.C. Purshouse. Evolutionary algorithms in control system engineering: a survey. Control Engineering Practice, 10:1223–1241, 2002. 6. H. Furuta, K. Maeda, and E. Watanabe. Application of genetic algorithm to aesthetic design of bridge structures. Microcomputers in Civil Engineering, 10(6):415–421, 1995.

12

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

7. T.A. Johansen. Multi-objective identification of fir models. In IFAC SYSID 2000. Santa Barbara, USA, 2000. 8. K. Kristinsson and G.A. Dumont. System identication and control using genetic algorithms. IEEE Transactions on Systems, Man, and Cybernetics, 22(5):1033– 1046, 1992. 9. J. McCormack. Interactive evolution of l-system grammars for computer graphics modelling. In D.G. Green and T. Bossomaier, editors, Complex Systems: from Biology to Computation, pages 118–130. IOS Press, Amsterdam, Netherlands, 1993. 10. C. Moons and B.De Moor. Parameter identification of induction motor drives. Automatica, 31:1137–1147, 1995. 11. M. Papadrakakis N.D. Lagaros and G. Kokossalakis. Structural optimization using evolutionary algorithms. Computers and Structures, 80, 2002. 12. R.W. Noel and S. Acchione-Noel. Objet trouv, holism, and morphogenesis in interactive evolution. CogSci, 2002, Fairfax Virginia, Aug., 2002. 13. M. Ohsaki and H. Takagi. Application of interactive evolutionary computation to optimal tuning of digital hearing aids. In Int’l Conf. on Soft Computing (IIZUKA’98), Iizuka, Fukuoka, Japan. World Scientic, 1998. 14. A. Pages, H. Pinguad, M. meyer, and X. Joulia. From process simulation to general estimation. Computers and Chemical Engineering, 20:1413–1424, 1996. 15. I. Rechenberg. Case studies in evolutionary experimentation and computation. Comput. Methods Appl. Mech. Engrg., 186:125–140, 2000. 16. H.P. Schwefel. Numerical Optimization of Computer Models. John Wiley and Sons, New York, 2nd edition, 1995. 17. J.R. Smith. Designing biomorphs with an interactive genetic algorithm. In 4th Int’l Conf. on Genetic Algorithms (ICGA’91), San Diego, CA, USA, pages 535–538. Morgan Kaufmann Publisher, 1991. 18. W.M. Spears, K.A. De Jong, T. Back, D.B. Fogel, and H. Garis. An overview of evolutionary computation. In Proceedings of European Conference on Machine Learning. 1993. 19. H. Takagi. Interactive evolutionary computation - cooperation of computational intelligence and human kansei. In 5th Intl Conf. on Soft Computing (IIZUKA98), pages 41–50. World Scientic, Iizuka, Fukuoka, Japan, 1998. 20. H. Takagi. Interactive evolutionary computation: System optimization based on human subjective evaluation. In IEEE Intl Conf. on Intelligent Engineering Systems (INES98), pages 1–6. Vienna, Austria, 1998. 21. H. Takagi. Interactive evolutionary computation: Fusion of the capabilities of ec optimization and human evaluation. IEEE, 89(9), 2001. 22. K.C. Tan and Y. Li. Grey-box model identification vie evolutionary computing. Control Engineering Practice, 10:673–684, 2002. 23. M.L. Thompson and M.A. Kramer. Modeling chemical processes using prior knowledge and neural networks. AIChE J., 40:1328–1340, 1994. 24. H.J.A.F. Tulleken. Grey-box modelling and identification using physical knowledge and bayesian techniques. Automatica, 29:285–308, 1993.

Key words: System identification, Interactive Evolutionary Computation

1 Introduction The main objective of system identification is to identify a model with good prediction capabilities in the sense that it is able to accurately predict the system’s response to a given class of excitations. Hence, a common identification objective is to minimize the mismatch between model prediction and observed data. However, it is often desirable to introduce additional objectives and constraints into the identification problem. There are several reasons for this [7]: In most practical cases the data sequence used for identification may be incomplete or uncertain, the model structure is not necessarily known and can be overparameterized, there may be some properties that we may want the model to have. Furthermore, there is often additional information available, and it makes sense to use this information to improve the accuracy and validity of the model [1]. Multi-objective system identification tools allow the incorporation of these objectives into the identification procedure. E.g.,

2

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

single constraints or penalties on the parameter space have been suggested, for example to ensure stability [24, 1], fulfillment of balance equations and steady-state data [23], and explicit belief about parameter values [14, 10]. These multi-objective identification problems involve the simultaneous consideration of multiple performance criteria. These objectives are often noncommensurable and are frequently in conflict with one another [7]. Trade-offs exist between some objectives where advancement in one objective will cause deterioration in another. In most cases the cost function is defined prior to the optimization procedure; this requires in-depth information concerning the various trade-offs and valuation of each individual objective. However, in some cases, the objective functions are explicitly/mathematically not available and the identification problem can be decomposed only for large-scale problems. This paper proposes the application of Interactive Evolutionary Computation (IEC) [21] to handle these identification problems, and its application to the identification of tendency models of batch reactors. In the proposed framework, human users evaluate the performances of the individuals (set of solutions) by selecting the best individuals. The human user does not (only) evaluate the mean square fitting error, or time constant of the model but he or she can simultaneously analyze these numerical values with the plotted dynamic and/or steady-state behavior of the model, e.g. step-responses, gains, etc. The proposed approach has been implemented in MATLAB (EAsy-IEC Toolbox) and applied to the identification of a pilot batch reactor. The remaining part of the paper is organized as follow. In Section 2 the algorithm of Evolutionary Strategy is presented, which is tailored taking into account the specialities of IEC. Section 3 shows an application example when the proposed approach is applied to identification of a tendency model of a heating-cooling jacket. Finally, some conclusions are drawn in Section 4.

2 Description of the Algorithm 2.1 Introduction to Evolutionary Algorithms Evolutionary Algorithms (EAs) [18] are stochastic optimization algorithms that mimic the process of natural selection. Evolutionary Algorithms handle several candidate solutions simultaneously. Thus, unlike conventional algorithms, EAs start form a pool of points (candidate solutions), usually referred to as individuals. Every individual represents a point in the search space. Every individual has a fitness value, which is calculated by the evaluation function. The fitness value determines the successfulness of the given individual appropriate to the objectives. The next generation of individuals is generated from the current population using genetic operators: selection, mutation and crossover. Likes in the natural world, the more successful individuals, which have bigger fitness values, have higher probability to transmit their genes into

Interactive Evolutionary Computation

3

the new generation. Because the genetic operators are stochastic, EAs are stochastic optimization methods. In the initialization step, EAs generate individuals randomly. In every evolutionary-loop, all of the individuals are evaluated and their fitness values are calculated based on an objective function. After that the selection operator selects randomly the best individuals to constitute the next generation. The crossover (recombination) operator generates two new individuals by taking two selected individuals and recombining them. The mutation operator changes randomly the individuals. Finally the old individuals are replaced by the new individuals. Generally, EAs are automated algorithms using a given objective function to calculate the fitness values of the individuals and the selection operator performs automatically the selection based on these calculated fitness values (for example using roulette wheel selection). But sometimes instead of a given objective function, a human user can make the valuation of the individuals, which results in an interactive optimization. Interactive Evolutionary Computation (IEC) is an optimization method that adopts evolutionary computation among system optimization based on subjective human evaluation [21]. In contrast to the automated evolutionary optimization algorithms, IEC cannot use many individuals and searching generations because of resulting human fatigue. Hence, there is a need to develop and tune an algorithm which effectively search with a few individuals within a few searching generations. For this purpose, this paper presents an algorithm based on Evolutionary Strategy (ES). In the following the IEC and the ES will be described and then the structure of the developed IEC-Tune Toolbox will be shown. 2.2 Interactive Evolutionary Computation Most optimization techniques which work by improving a single solution step by step are not suited for interactive optimization approach. For example, in simple hill-climbing a user suggestion can either be accepted, which means the loss of results gained so far, or denied, which may lead the algorithm a wrong way. Evolutionary Algorithms and other population based optimization procedures are better suited for interactive optimization, because a suggestion can be put into the population to see whether it is worth further consideration or not. A good suggestion can be expected to survive and lead the optimization process towards a good solution, while a bad suggestion will be eliminated quite soon. This approach become well known as Interactive Evolutionary Computation (IEC) [20, 19]. Simply stated, IEC is a technique from the class of evolutionary algorithms (EAs), in which the fitness function is replaced by a human user. As in interactive evolution, the user selects one or more individual(s) which survive(s) and reproduce(s) (with variation) to constitute a new generation. Hence IEC uses two different spaces for its search: the human user evaluates the output

4

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

of the target system according to the distance between the target goal and the system output in psychological space [20], on the other hand, the EA searches in the parameter space. These techniques were applied already in computer graphics [9], animation, creating forms, textures, and motion [2]. Potential applications of interactive evolution include artificial life design, e.g., development of components of biological nature [17, 6] and engineering construction design [13]. Although, IEC has been applied in knowledge acquisition and data mining, according to our best knowledge, the applicability of this approach in system identification has not been studied yet. However, in recent years, the application of Evolutionary Computing to black-box and grey-box model identification has received considerable interest since the seminal paper by Kristinsson and Dumont [8]. Subsequent evolutionary system identification applications are applied to build dominant structural identification with local parametric tuning without the need of a differentiable performance index in the presence of noisy data. As the overview of Fleming [5] illustrates evolutionary techniques provide an excellent fitting performance and is capable of accommodating multiple objectives such as to examine the relationships between model complexity and fitting accuracy during the model building process. These approaches can be extended to the identification of gray-box models, to use of a priori knowledge [22]. 2.3 Evolution Strategy In contrast to automated evolutionary optimization algorithms, the IEC cannot use many individuals and searching generations because the human fatigue. Hence, there is a need to develop an algorithm which can effectively search with a few individuals within a few searching generations. For this purpose, this paper presents an algorithm based on Evolutionary Strategy. Evolution Strategy (ES) was developed by [15], with selection, mutation, and a population of size one. Schwefel [16] introduced recombination and populations with more than one individual, and provided a nice comparison of ES with more traditional optimization techniques. Evolutionary Strategy are typically applied to real-valued parameter optimization problems. The characteristic feature of ES lies in the self-adaptation of the standard deviation of the Gaussian distribution used in the mutation [3]. The basic idea is to add these adaptation parameters to the genotype and have them undergo evolution themselves. The main elements of the algorithm used in this paper are the following: Representation An ES-individual aj = (xj , σ j ) consists of two components: the object variables xj = [xj,1 , . . . , xj,n ] and the strategy variables σ j = [σj,1 , . . . , σj,n ]. The strategy parameters allow the algorithm to adapt to the topology of the objective function by controlling the mutation step size.

Interactive Evolutionary Computation

5

Mutation As small changes occur frequently but large ones only rarely in the nature, the mutation operator adds normal distributed random numbers zj,i ∼ N (0, σj,i ) to the individuals: xj,i = xj,i + zj,i . (1) Before the object variables are changed, the standard deviations are mutated using a multiplicative normally distributed process: (t)

(t−1)

σj,i = σj,i

exp(τ 0 N (0, 1) + τ Ni (0, 1)),

(2)

with exp (τ 0 N (0, 1)) as a global factor which allows an overall change of the mutability and exp (τ Ni (0, 1)) allowing for individual changes of the mean step sizes σj,i . The τ 0 and τ parameters can be interpreted in the sense of global learning rates. Schwefel suggests to set them as [16]: 1 1 τ0 = √ , τ = p √ . 2n 2 n

(3)

Recombination Recombination in ESs can be either sexual, where only two parents are involved in the creation of an offspring, or global, where up to the whole population contributes to a new offspring. Sexual recombination of just two individuals is often called local while the contribution of all individuals is called global recombination. Traditional recombination operators are discrete recombination, intermediate recombination, and geometric recombination, all existing in a sexual and global form. When F and M denote two randomly selected individuals from the µ parent population, the following operators can be defined: xF,i no recombination xF,i or xM,i discrete 0 xi = (4) (x + x ) /2 intermediate F,i M,i Pµ global avarage k=1 xK,i /µ σF,i no recombination σ or x discrete F,i M,i + σ ) /2 intermediate σi0 = (σ . (5) M,i pF,i (σ σ ) geometric Pµ F,i M,i global avarage k=1 σK,i /µ Throughout this work discrete recombination of the object variables and global intermediate combination of the strategy parameters were used.

6

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

Selection and Replacement In the applied (µ + λ) Evolution Strategy, at a given generation (t), there are µ parents, and λ offspring are generated by recombination and mutation. Then, the µ+λ members of then next population are sorted according to their objective function values given by the user, and he or she selects the best µ to become the parents of the next generation. 2.4 IEC Integration to MATLAB IEC systems seek to interface EA to human preference in order to create systems capable of evolving artifacts that require a human expertise that has not yet succumb to computation. Hence, the interfacing of human ability with machine computation requires resolving difficult issues [12]. Further, progress in the design of interactive evolutionary systems allows a glimpse into how very human abilities such as intuition, projection, and holistic perception interplay with the mechanics of machine computation. This section reports on one such interactive evolutionary system that seeks to combine human perception with the previously presented evolutionary strategy to evolve the parameters of MATLAB and Simulink models (see Figure 1). The EAsy-IEC Tune Toolbox is downloadable from the website of the the author: www.fmt.vein.hu/softcomp/EAsy, is designed to be applicable for different types of optimization problems (e.g., system identification, controller tuning, data mining). The number of displayed individuals, µ + λ is set to 9, which can be displayed spatially. The number of searching generations is limited to 20 generations at the most, due to the fatigue of human operators. In the proposed framework human users evaluate the performances of the individuals (set of solutions). Not the genotypes and phenotypes are applied directly, but rather the output of the target system realized by the individuals. For example, the user does not (only) evaluate the mean square fitting error, or time constant of the model but he or she can simultaneously analyze these numerical values with the plotted dynamic and/or steady state behavior of the model, e.g.

Fig. 1. Scheme of the Proposed IEC System Integrated to MATLAB

Interactive Evolutionary Computation

7

step-responses, gains, etc. An example of such evaluation figure is given in Figure 4. Due to small number of individuals, the human evaluation was restricted to selecting the individuals which will survive and constitute the next generation.

3 Application Example 3.1 Process Description In recent years our department focuses to the development of control strategies of multipurpose production plants. Hence, a prototype of pharmaceutical process systems was designed and installed in our laboratory. The central element of this process unit is a 50 liter stirred reactor with heating-cooling jacket (see fig. 2).

Fig. 2. The stirred tank reactor with heating-cooling jacket

This paper focuses to the identification of the tendency model of the heating-cooling jacket. Through the jacket the direct heating-cooling system allows cooling with chilled water or heating with steam in steam or hot water mode. In the hot water heating mode the water is circulated by a pump, the steam is introduced through a mixer, while the excess water is removed through an overflow. 3.2 Tendency Model of the Heating-Cooling Jacket The jacket and the circulating liquid can be described using the common lumped parameter enthalpy or heat balances given in the chemical engineering literature. In the model, zero-volume distributors and mixers and are applied and the overflow as well as the feed of the steam and the fresh cooling water

8

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

are taken into account. The obtained simplified first principle model can be regarded as the tendency model [4] of the most important phenomena and its scheme is given on Figure 3.

Fig. 3. Tendency model of the jacket of the stirred ractor

V 1 and V 2 are the valve positions (0%-100%), T H is the temperature of the cooling water coming from the environment, T F is the equivalent steam temperature calculated from the boiling temperature, the latent heat and specific heat. The valve characteristics are given in form of second order polynomials (V ALV E1 and V ALV E2). The two first order transfer functions 1 1 ( T 1s+1 and T 2s+1 ) are obtained from the lumped parameter model of the jacket and the model of the thermometer at the jacket inlet. The K1, K2, and K3 are gains and the DT is the dead-time. The identification problem is identifying these parameters of the model (gains, time constants, gains, dead-time and parameters of second order polynomials). For this purpose open-loop experiments were conducted on the whole operating range of the composite manipulated variable u (because of the split-range control: V 1 and V 2 are not opened at the same time, u = 100% means that V 1 is opened fully, u = 0% means that V 2 is opened fully). 3.3 Results Based on the collected input-output data, the identification problem can be formulated as an optimization problem and can be solved using common nonlinear optimization methods (e.g. SQP). To formulate an optimization problem we need a cost-function which is minimized under the optimization procedure. Usually for this purpose the mean square model error is used, n P 2 M SE = (y(k) − y˜(k)) /n, where y is the observed output, y˜ is the prek=1

dicted output of the model, n is the number of observation. Beside the MSE, other objectives should be considered. E.g. the steadystate prediction error and the characteristics of the valves. But it is difficult to formulate these objectives and it is difficult to select the weights of the

Interactive Evolutionary Computation

9

objectives in the cost-function. Hence it is problematical to introduce these additional objectives into the identification problem. The IEC handles this problem, because the human user can consider and balance these objectives. For comparison, beside the proposed IEC optimization scheme, the direct minimization of the MSE prediction error is considered. For this purpose, Sequential Quadratic Programming (SQP) was used. The SQP method is a standard general purpose algorithm for solving NLP optimization problems. It is a gradient based method which calculates the search direction vector from the quadratic approximation of the objective function around the actual local point [11]. The SQP based identification method minimizes automatically the cost-function subject to the parameters of the model. The cost-function was evaluated by simulating the model and calculating the MSE. In contrast to SQP based method, the user evaluates the model in the IEC based identification method. Firstly, the identified model was simulated and then the human user selected the best solutions based on the dynamic and the near steady-state behavior of the system, and the identified characteristic of the valves. Beside these plots, some numerical values were also given to the user to glimpse. Such an evaluation display is depicted in Figure 4. 100

100

100

100

100

100

100

100

50

50

50

50

50

50

50

50

50

Temp. [C]

100

0 0

Temp. [C]

100

0 0 0 0 0 0 0 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] 100 100 100 100 100 100 100 100

50

0

F1%, F2%

100

50

50

50

50

50

50

50

0 0 0 0 0 0 0 0 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 100 150 200 Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] Time [sec] 100 100 100 100 100 100 100 100

50

0

50

50

50

50

50

50

50

50

50

0 0 0 0 0 0 0 0 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 0 50 100 V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% V1%, V2% ISE: 15.1 K1: 0.055 K2: 0.357 K3: 0.025 T0: 0.559 Th: 0.422

ISE: 9.3 K1: 0.055 K2: 0.382 K3: 0.024 T0: 0.560 Th: 0.424

ISE: 38.8 K1: 0.052 K2: 0.382 K3: 0.024 T0: 0.536 Th: 0.424

ISE: 101.6 K1: 0.038 K2: 0.375 K3: 0.016 T0: 0.392 Th: 0.417

ISE: 339.2 K1: 0.065 K2: 0.387 K3: 0.030 T0: 0.687 Th: 0.428

ISE: 42.4 K1: 0.043 K2: 0.378 K3: 0.019 T0: 0.438 Th: 0.420

ISE: 213.2 K1: 0.024 K2: 0.344 K3: 0.010 T0: 0.229 Th: 0.411

ISE: 409.8 K1: 0.013 K2: 0.340 K3: 0.010 T0: 0.143 Th: 0.404

ISE: 3948.8 K1: 0.071 K2: 0.390 K3: 0.033 T0: 0.719 Th: 0.433

Fig. 4. MATLAB Evaluation figure in the IEC-Tune toolbox

From an arbitrary initial value, the application of IEC resulted in 5.0 MSE prediction error after 14 generations. This corresponds to 146 evaluations of the model. If the same numbers of function calls are applied in SQP, the obtain model results in 26.0 MSE. The drawback of IEC is that is slower than the SQP method due to the time-requirement of the human evaluation. Hence while the IEC procedure consumed about 8.5 minutes, the SQP consumed only 2.5 minutes. (See the Table 1. and Figures 5, 6.)

10

J. Abonyi, J. Madar, L. Nagy and F. Szeifert Table 1. Comparative results Evaulation number 146 149 506

Time (sec) 508 150 501

MSE 5.0 26.0 1.5

Temp. [C]

100

50

0

0

10

20

30 Time [sec]

40

50

60

Temp. [C]

100

50

0

80

100

120

140 Time [sec]

40

50 V1%, V2%

160

180

200

F1%, F2%

100

50

0

0

10

20

30

60

70

80

90

100

Fig. 5. Results using IEC after 146 model evaluations

Temp. [C]

100

50

0

0

10

20

30 Time [sec]

40

50

60

Temp. [C]

100

50

0

80

100

120

140 Time [sec]

40

50 V %, V %

160

180

200

100 F1%, F2%

Method IEC SQP SQP

50

0

0

10

20

30

1

60

70

80

90

100

2

Fig. 6. Results using SQP after 149 model evaluations

Interactive Evolutionary Computation

11

Our experiences and these results show that the proposed IEC approach is extremely useful to find good initial conditions to SQP that is sensitive to stacking in local minima.

4 Conclusions This paper demonstrated how interactive evolution algorithms can be applied to solve multi-objective identification problems. A MATLAB toolbox has been developed to provide a human-machine interface of the IEC. With the use of this tool, the user can easily select the individuals which then are reproduced and mutated to constitute the next generation. The proposed approach applied to the identification of a pilot batch reactor. The results show that IEC is an efficient and comfortable method to incorporate the priori knowledge and preferences of the user into the identification problem. Our experiences show, IEC is extremely useful to find good initial conditions to other nonlinear optimization algorithms, like SQP.

Acknowledgement The authors would like to acknowledge the support of the Cooperative Research Center (VIKKK) (KKK-I-7), the Hungarian Ministry of Education (FKFP-0073/2001), and the Hungarian Science Foundation (OTKA TO37600). Janos Abonyi is grateful for the financial support of the Janos Bolyai Research Fellowship of the Hungarian Academy of Science.

References 1. J. Abonyi, R. Babuska, H.B. Verbruggen, and F. Szeifert. Constraint parameter estimation in fuzzy modelling. In FUZZ-IEEE’99 Conference, pages 951–956. Seoul, Korea, 1999. 2. K. Aoki, H. Takagi, and N. Fujimura. Interactive ga-based design support system for lighting design in computer graphics. In Int’l Conf. on Soft Computing (IIZUKA’96), Iizuka, Fukuoka, Japan, pages 533–536. World Scientic, 1996. 3. A.E. Eiben and M. Schoenauer. Evolutionary computing. Information Processing Letters, 82, 2002. 4. C. Filippi-Bossy, C. Georgakis, J. Bordet, J. Villermaux, and S. MarchalBrassely. Batch reactor optimization by use of tendency models. Comp. & Chem. Eng, 13:35–47, 1989. 5. P.J. Fleming and R.C. Purshouse. Evolutionary algorithms in control system engineering: a survey. Control Engineering Practice, 10:1223–1241, 2002. 6. H. Furuta, K. Maeda, and E. Watanabe. Application of genetic algorithm to aesthetic design of bridge structures. Microcomputers in Civil Engineering, 10(6):415–421, 1995.

12

J. Abonyi, J. Madar, L. Nagy and F. Szeifert

7. T.A. Johansen. Multi-objective identification of fir models. In IFAC SYSID 2000. Santa Barbara, USA, 2000. 8. K. Kristinsson and G.A. Dumont. System identication and control using genetic algorithms. IEEE Transactions on Systems, Man, and Cybernetics, 22(5):1033– 1046, 1992. 9. J. McCormack. Interactive evolution of l-system grammars for computer graphics modelling. In D.G. Green and T. Bossomaier, editors, Complex Systems: from Biology to Computation, pages 118–130. IOS Press, Amsterdam, Netherlands, 1993. 10. C. Moons and B.De Moor. Parameter identification of induction motor drives. Automatica, 31:1137–1147, 1995. 11. M. Papadrakakis N.D. Lagaros and G. Kokossalakis. Structural optimization using evolutionary algorithms. Computers and Structures, 80, 2002. 12. R.W. Noel and S. Acchione-Noel. Objet trouv, holism, and morphogenesis in interactive evolution. CogSci, 2002, Fairfax Virginia, Aug., 2002. 13. M. Ohsaki and H. Takagi. Application of interactive evolutionary computation to optimal tuning of digital hearing aids. In Int’l Conf. on Soft Computing (IIZUKA’98), Iizuka, Fukuoka, Japan. World Scientic, 1998. 14. A. Pages, H. Pinguad, M. meyer, and X. Joulia. From process simulation to general estimation. Computers and Chemical Engineering, 20:1413–1424, 1996. 15. I. Rechenberg. Case studies in evolutionary experimentation and computation. Comput. Methods Appl. Mech. Engrg., 186:125–140, 2000. 16. H.P. Schwefel. Numerical Optimization of Computer Models. John Wiley and Sons, New York, 2nd edition, 1995. 17. J.R. Smith. Designing biomorphs with an interactive genetic algorithm. In 4th Int’l Conf. on Genetic Algorithms (ICGA’91), San Diego, CA, USA, pages 535–538. Morgan Kaufmann Publisher, 1991. 18. W.M. Spears, K.A. De Jong, T. Back, D.B. Fogel, and H. Garis. An overview of evolutionary computation. In Proceedings of European Conference on Machine Learning. 1993. 19. H. Takagi. Interactive evolutionary computation - cooperation of computational intelligence and human kansei. In 5th Intl Conf. on Soft Computing (IIZUKA98), pages 41–50. World Scientic, Iizuka, Fukuoka, Japan, 1998. 20. H. Takagi. Interactive evolutionary computation: System optimization based on human subjective evaluation. In IEEE Intl Conf. on Intelligent Engineering Systems (INES98), pages 1–6. Vienna, Austria, 1998. 21. H. Takagi. Interactive evolutionary computation: Fusion of the capabilities of ec optimization and human evaluation. IEEE, 89(9), 2001. 22. K.C. Tan and Y. Li. Grey-box model identification vie evolutionary computing. Control Engineering Practice, 10:673–684, 2002. 23. M.L. Thompson and M.A. Kramer. Modeling chemical processes using prior knowledge and neural networks. AIChE J., 40:1328–1340, 1994. 24. H.J.A.F. Tulleken. Grey-box modelling and identification using physical knowledge and bayesian techniques. Automatica, 29:285–308, 1993.