Multi-Objective parameter estimation problems: an improved strategy

Abstract

A multi-objective optimization approach has been applied to solve parameter estimation problems. An improved algorithm, based on evolutionary strategies, has been proposed to optimize mathematical model parameters. The algorithm makes use of a new concept of fitness function, which determines the reproduction ratio as a function of the population density, and a new class of operators, which enhance the algorithm performance.

Two processes have been analyzed: a grain cooling process and a grain drying process. In order to estimate the coefficient of heat transfer and the drying rate parameters of these models, minimization of the sum of the least squares of temperature and equilibrium moisture content have been conducted. Experimental data obtained from the soybean cooling in a continuous cross-flow moving bed heat exchanger and the corn drying in a fixed bed dryer have been used to evaluate the estimated parameters. The simulated results demonstrated the algorithm efficiency to perform parameter estimation. The validated model consistently fits the experimental data.

Keywords:

• Parameter estimation
• Multi-objective optimization
• Evolutionary algorithm
• Drying process
• Cooling process

Introduction

Accurate modeling is an important prerequisite for optimization, control and design of engineering processes. Mathematical models are sets of expressions that simulate the behavior of a system, describing quantitatively the mechanisms of the process. The generality of a model depends upon the process complexity as well as the available information concerning the system. In order to validate the model, it is often necessary to estimate the model parameters of nonlinear algebraic or differential equations.

The parameter estimation procedure is based on the minimization of the total sum of every observed error. Such procedure is often conducted using optimization algorithms. Traditionally, the estimation methods are “hill-climbing” methods that converge to only one stationary point. The objective functions based upon a nonlinear model and experimental data frequently present more than one optimum. The hill-climbing algorithms are easily trapped in local optimum, and are unable to exploit the existing optima when the eigenvalue ratio becomes large [1]. Heuristic optimization techniques have been successfully used to solve parameter estimation problems [1–4].

Parameter estimation constitutes a multi-objective optimization problem, as it entails opposing requirements to be satisfied simultaneously. As the targets are often in conflict, a compromise solution has to be sought. The estimation problem has been handled by single objective approaches, in which all targets are combined into a single objective function by using weighting factors. A meaningful optimization approach implies the use of the multi-objective concept, which generates a family of equally good solutions, called Pareto optimal set, instead of a global optimum. The choice of the optimal results requires additional knowledge of the process, in order to select the best compromise solution between the competing goals.

This contribution is focused on the development of a multi-objective optimization algorithm to solve parameter estimation problems. An improved approach based on evolutionary strategies is proposed. A new class of operators consisting of a Pareto-set filter, an elitism operator and a niche test operator is introduced. These operators enhance the performance of the algorithm and are considerably simple to implement. A new ranking strategy is proposed to treat multidimensional problems. The new strategy does not require any previous knowledge of the relative importance of individual objectives. A fitness function based on each rank population size and rank level is also provided to determine the reproduction ratio.

In order to evaluate the algorithm performance, two different processes have been studied: a soybean grain cooling process and a corn grain drying process. The mathematical model of each process was employed to generate the dependent variable values under the same conditions of the experiments. The minimization of every discrepancy between the measured and the predicted variables guides the search for the best-fit parameters. The heat and moisture transfer coefficients of these processes are estimated by the temperature and equilibrium moisture content least squares. Experimental data obtained from a continuous cross-flow moving bed heat exchanger and a fixed bed dryer have been used to evaluate the estimated parameters.

The simulated results demonstrated the ability of the method to handle multi-objective optimization problems. The proposed algorithm can successfully seek tradeoff surface regions and find the Pareto optimal set. It has been shown that the validated model fits the experimental data satisfactorily, indicating to be a consistent strategy to solve parameter estimation problems.

Multi-Objective Evolutionary Algorithm

Evolutionary algorithms simulate a natural evolution process: the fittest species survive and propagate while the less successful tend to disappear. The optimization procedure consists of a search for nondominated solutions. The concept of nondominance refers to the solutions for which no objective can be improved without worsening at least one of the other objectives. The progress strategy is guided by the fitness evaluation, and consists of performing the population with genetic operators to generate the next population. Different adaptations of the evolutionary techniques are presented in the literature [5–8]. A detailed background on the evolutionary theory is reported in [9, 10].

The Proposed Algorithm

The multi-objective optimization algorithm developed is based on the concept of non-dominance, proposed by Goldberg [9]. The ranking procedure is extended to treat multidimensional problems. New operators are introduced to enhance the algorithm performance: (a) a niche test operator, which prevents genetic drift and maintains a uniformly distributed population along the optimal set; (b) a Pareto-set filter, which avoids missing optimal points during the evolutionary process and (c) an elitism operator, which insures the propagation of the best result of each individual objective function. These operators reduce the necessary number of generations and are computationally feasible for even complicated problems [11]. A new concept of fitness function is adopted in order to provide the evolutionary process with robust and stable convergence. The computational code developed operates in a continuous variable space, which is computationally fast and stable in converging to the Pareto optimal set.

The proposed GA procedure works through the following steps:

1.	creation of a random initial population;
2.	evaluation of the individuals;
3.	ranking of the individuals, calculation of the fitness, registration of the best individuals;
4.	registration of all nondominated individuals in the filter operator;
5.	selection of pairs of individuals as parents;
6.	crossover of the parents to generate the children;
7.	replacement of the individuals using the niche test operator;
8.	genetic mutation;
9.	replacement of the individuals using the elitism operator.

The improvements made to the original method are discussed in the next sections.

Ranking Procedure

The ranking procedure consists of the classification of the individuals into categories according to the concept of dominance. First, all nondominated individuals of the population are identified and assigned rank 1. These individuals are virtually removed from the population and a new evaluation is conducted on the remaining individuals. The next set of nondominated points are identified and assigned rank 2. This procedure continues until all the individuals are classified. Figure 1 illustrates the ranking procedure for a maximization problem.

FIGURE 1 Population ranking for a two-dimensional maximization problem.

Mathematically, the classification method proposed by the authors is conducted in four basic steps:

1.	the points are sorted according to the evaluation of an objective function, randomly chosen as a reference function;
2.	all points that produce opposite effects on the function values are selected as potential candidates to the rank, i.e., the sequence of points that simultaneously produces an increase on the reference function and a decrease on the other functions, or vice versa, will be selected;
3.	repeat Steps (a) and (b) until all the objective functions have been chosen as the reference function;
4.	all points selected in the Step (b) will be assigned to the rank.

These steps are repeated for each rank, until all points are classified.

Fitness Function

After being ranked, all individuals in the population are evaluated by a fitness function. The fitness value represents a measure of each individual performance, which will be used on the selection procedure. Every individual belonging to the same rank class is considered equivalent and has the same fitness value. Thus, it has the same probability of being selected for reproduction. The fitness function, F _k , is determined to each individual of the same rank k by the following equations, as proposed by Cheng and Li [5]:

where N _r is the highest rank of the population, P _s is the population size and P _sk is the population size of rank k. According to this fitness definition, the larger the population size at a rank is, the smaller the fitness of a point. Hence, the reproduction ratio of individuals at each rank depends on both the rank level and population size.

Niche Test Operator

The niche test operator consists of a replacement procedure of individuals. The test is applied to determine which individuals will go to the next population. Each child generated by the crossover procedure has its fitness calculated in the domain of the parental population. If the children’ superior fitness exceeds the best fitness of the parents, the children will replace the parents. Otherwise, both parents go to the next generation, if the number of children per crossover is set on two, or just the best parent will propagate, if this number is set on one. This operator helps to attenuate the phenomenon known as genetic drift, which makes a population become clustered at certain regions. The genetic drift results from stochastic errors associated with the genetic operators and the use of a population of finite size. The niche test operator is also adopted in order to maintain the appropriate diversity of the population [5].

Elitism Operator

The evolutionary process, as mentioned before, is guided by the ranking of the population. In the ranking procedure, points are categorized into groups according to their non-dominance levels. The closer to the Pareto set a group of points is, the higher its probability to propagate to the next generation. The concept of non-dominance refers to a characteristic behavior presented by a sequence of points sorted according to their objective function values. Hence, the classification of a group of points depends on their best objective function values.

The elitism procedure consists of the propagation of the best solutions to the next generation. For this purpose, it maintains an elitism file where the best solution of each individual objective function is registered. At each generation, this file is updated: if a better solution was generated it replaces the one stored. The individuals selected to the next population by the crossover and mutation procedures are submitted to the elitism operator. The points registered in the elitism file will randomly replace some of the individuals selected as candidates to the next generation. If the individual chosen to be replaced is better than the one registered in the elitism file, a new candidate is selected.

Since the elitism operator stores the best solutions obtained till the current iteration, it guarantees the propagation of the best solutions during the evolutionary search procedure. As a consequence, it increases the convergence of the optimization process as well as the robustness of the algorithm.

Pareto-set Filter Operator

The Pareto-set filter operator works during the whole optimization process. All points assigned rank 1 at each generation are registered in a Pareto-set file. This file is dynamically updated by using the filter operator. In this filter, the nondominated solutions of the current population are compared with those already stored in the file, from the previous generations. A new evaluation is conducted in all points in the filter, according to the following rules:

1.	all points in the filter reassigned rank 1, i.e., points identified as nondominated, are recorded in the Pareto-set file. The dominated ones are discarded;
2.	if the number of points in the file is inferior to the population size, the new nondominated points are stored. Otherwise, if the file is full, the most similar points in the Pareto-set file are replaced.

The concept used to measure the similarity between two points in rule (b) is based on the euclidean distance between the values of the objective functions. Points with minimum distance relative to the others are removed from the file. Such procedure maintains an even distribution of the points in the file, which helps to provide the optimal set with the best solutions. At the end of the optimization process, the file itself comprises the Pareto optimal set and constitutes the result of the optimization.

The proposed algorithm is detailed in a previous work [11]. The computational load associated to the algorithm relies on the number of objective functions of the problem. At each generation, the objective functions are evaluated just once for each point. Therefore, the total number of the objective evaluations is a function of the size of the population, the number of generations and the number of objective functions optimized.

Parameter Estimation Problem

Grain Cooling Process

The first case study is the estimation of the overall heat transfer coefficient in a grain cooling process. The process is conducted in a continuous cross-flow moving bed heat exchanger, fed with warm soybean grains and dry air. The scheme of the process unit is summarized in Fig. 2.

FIGURE 2 Scheme of the grain cooling process.

The model formulation considers the convection on the particle surface as the dominant mechanism in the gas–solid heat transfer [12]. Some assumptions were made to simplify the mathematical formulation of the problem:

1.	gas and solid velocity profiles are linear;
2.	the temperature profiles in both phase feeding regions are considered linear;
3.	heat losses through the equipment walls are negligible;
4.	the heat transfer perpendicular to the flow is negligible comparing to the heat transfer in the flow direction;
5.	the bed properties are homogeneous;
6.	solid and gas flows are unidirectional;
7.	physical and chemical properties are constant throughout the bed.

A hyperbolic system of differential equations are provided to describe the steady-state heat transfer mechanism as follows: Gas phase:

Solid phase:

where S = (6(1 - ϵ))/(ϕ Dp) Subject to the initial conditions

The physical chemistry properties of the soybean and the air are presented in Table I. The operational conditions of the problem are shown in Table II.

TABLE I Solid and gas physical chemistry properties [12]

ε (−)	3.778 × 10^−1
ϕ (−)	0.9
Dp (cm)	0.6
ρ (g/cm^3)	1.17
Cp s (cal/g C)	0.53
Cp g (cal/g C)	0.25

TABLE II Operational conditions [12]

T g 0 (°C)	39.8
T s 0 (°C)	60.0
G g (g s/cm^2)	8.63 × 10^−2
G s (g s/cm^2)	4.38 × 10^−2

After integration, Eqs. (3) and (4) can be expressed in dimensionless form as:

where I ₀ is the modified Bessel function of the first kind of order zero, and

The optimization problem involves the minimization of the least squares of the measured temperature and the one predicted by the model. The heat transfer coefficient is the decision variable. The parameter estimation problem is formulated as follows:

Experimental results obtained in a pilot scale heat exchanger have been used to estimate the heat transfer coefficient.

In order to investigate the multi-objective behavior of the formulated problem, the decision variable is plotted in Fig. 3 as a function of the objective functions. The y-axis was rescaled to enhance the visualization of the multi-objective effect.

FIGURE 3 Heat transfer coefficient of the grain cooling process as a function of the objective functions F ₁ and F ₂.

The objective functions are affected in opposing ways by changes in the decision variable. Therefore, the estimation parameter of the soybean cooling process constitutes a multi-objective optimization problem, and is suitable to evaluate the algorithm performance.

Grain Drying Process

The second problem deals with the estimation of the heat transfer coefficient and the drying rate parameters of a corn drying process, conducted in a fixed bed dryer. The scheme of the process unit is shown in Fig. 4.

FIGURE 4 Scheme of the grain drying process.

The drying model constitutes a set of differential equations describing the temperature and moisture content profiles in the gas and solid phases [13]. An algebraic–differential approach has been adopted to express the energy balances in terms of the enthalpies. The hypothesis considered in the model formulation involve [14]:

1.	a uniform plug flow of air through the bed;
2.	adiabatic walls, the heat losses are negligible;
3.	the temperature gradients within the individual particles are negligible;
4.	the heat capacities are constant;
5.	the packing density of the grain in the bed is uniform;
6.	the evaporating water extracts the heat of desorption from the grain and enters the air stream as water vapor at the grain temperature.

The differential equations describing the heat and mass transfer mechanisms are presented: Solid phase mass balance

Solid phase energy balance

Gas phase mass balance

Gas phase energy balance

The enthalpies are related to the temperatures as:

The physical properties of the corn grains and the air are summarized in Table III. The operational conditions of the dryer as well as the parameters of the process used in the simulations are presented in Table IV.

TABLE III Solid and gas physical chemistry properties [15]

ε (−)	0.36
ϕ (−)	0.75
ρ s (kg/m^3)	1.23 × 10^3
ρ g (kg/m^3)	1.17
Dp (cm)	0.0078
Cp s 0 (kJ/kg K)	1.90
Cp v 0 (kJ/kg K)	1.89
Cp g 0 (kJ/kg K)	1.01
Cp l 0 (kJ/kg K)	4.186
λ TR	2.4584 × 10^3

TABLE IV Operational conditions of the dryer [15]

T g 0 (K)	363.15
T s 0 (K)	298.15
G g (kg/m^2 s)	0.55
Y s 0 (−)	0.218
Y g 0 (−)	0.016
t (s)	21600
d 1	−3.82 × 10^−5
d 2	1.8
d 3	410.0

The drying rate expression for corn grains is described by Mancini [16] as

where b ₁ and b ₂ are the parameters to be estimated in order to determine the drying rate. The grain moisture content in equilibrium, Y _se is expressed by the sorption isotherms for corn grains (T _s ≥ 40°C), described by Henderson's modified isotherm [15,17]:

The partial differential and algebraic equation system is discretized in the spatial domain into 50 points, using the implicit Euler form. The resulting system consists of 300 equations, which were integrated using the solver DASSL [18].

The formulation of the multi-objective parameter optimization problem comprises the minimization of the following objective functions:

Three decision variables are selected: the heat transfer coefficient, h, and the drying rate parameters of Eq. (19), b ₁ and b ₂. The experimental data obtained by Calçada [15] in a soybean drying process were used to estimate the parameters above. Additional parameters related to the process are presented in the same work.

Results and Discussion

The proposed algorithm has been applied to estimate the model parameters for the two cases formulated. The optimizations have been carried out using the algorithm parameters summarized in Table V, established according to a sensitivity study.

TABLE V Algorithm parameters

Population size	20 individuals
Crossover probability	75%
Mutation probability	5%
Number of children/crossover	1

Grain Cooling Process

Figure 5 shows the Pareto set obtained for the soybean cooling process. The number of generations required to optimize the heat transfer coefficient was 32 iterations. Table VI presents the numerical results of the optimization process.

TABLE VI Optimization results for the soybean cooling process

F 1 (−)	F 2 (−)	h (cal/cm^2 s°C)
2.784 × 10^−02	1.533 × 10^−02	3.322 × 10^−03
2.765 × 10^−02	1.533 × 10^−02	3.325 × 10^−03
2.746 × 10^−02	1.534 × 10^−02	3.328 × 10^−03
2.726 × 10^−02	1.535 × 10^−02	3.332 × 10^−03
2.709 × 10^−02	1.537 × 10^−02	3.335 × 10^−03
2.692 × 10^−02	1.540 × 10^−02	3.338 × 10^−03
2.672 × 10^−02	1.544 × 10^−02	3.343 × 10^−03
2.651 × 10^−02	1.550 × 10^−02	3.347 × 10^−03
2.630 × 10^−02	1.558 × 10^−02	3.352 × 10^−03
2.611 × 10^−02	1.567 × 10^−02	3.357 × 10^−03
2.598 × 10^−02	1.575 × 10^−02	3.361 × 10^−03
2.578 × 10^−02	1.590 × 10^−02	3.368 × 10^−03
2.569 × 10^−02	1.600 × 10^−02	3.372 × 10^−03
2.558 × 10^−02	1.613 × 10^−02	3.376 × 10^−03
2.549 × 10^−02	1.628 × 10^−02	3.381 × 10^−03
2.542 × 10^−02	1.642 × 10^−02	3.385 × 10^−03
2.537 × 10^−02	1.656 × 10^−02	3.390 × 10^−03
2.533 × 10^−02	1.673 × 10^−02	3.394 × 10^−03
2.530 × 10^−02	1.696 × 10^−02	3.400 × 10^−03
2.529 × 10^−02	1.717 × 10^−02	3.405 × 10^−03

FIGURE 5 Pareto optimal set for the soybean cooling process.

The optimal values of the heat transfer coefficient corresponding to both extremes of the Pareto set were used to verify the parameter estimation.

Figures 6 and 7 show that the model simulation using any of the estimated parameters perfectly fits the experimental data. In this case, the determination of the final solution should be based on additional information of the process, such as accuracy and reliability of the measurements.

FIGURE 6 Air temperature profiles as a function of the position in the heat exchanger.

FIGURE 7 Soybean temperature profiles as a function of the position in the heat exchanger.

Grain Drying Process

The Pareto set obtained for the corn drying process is depicted in Fig. 8. A total of 64 iterations were required to generate the optimal solution set. The results of the optimization as well as the estimated parameters are presented in Table VII.

TABLE VII Optimization results for the corn drying process

F 1 (−)	F 2 (−)	h (kJ/m^3 s)	b 1 (−)	b 2 (−)
5.151 × 10^−05	1.489 × 10^−02	63.230	1.442	6.296 × 10^−02
5.084 × 10^−05	1.490 × 10^−02	63.212	1.438	6.333 × 10^−02
4.903 × 10^−05	1.490 × 10^−02	63.257	1.402	5.939 × 10^−02
4.727 × 10^−05	1.497 × 10^−02	63.694	1.388	5.999 × 10^−02
4.657 × 10^−05	1.511 × 10^−02	63.392	1.385	6.097 × 10^−02
4.603 × 10^−05	1.540 × 10^−02	63.387	1.385	6.250 × 10^−02
4.562 × 10^−05	1.551 × 10^−02	63.357	1.380	6.249 × 10^−02
4.524 × 10^−05	1.577 × 10^−02	63.165	1.378	6.325 × 10^−02
4.497 × 10^−05	1.610 × 10^−02	63.269	1.379	6.427 × 10^−02
4.453 × 10^−05	1.641 × 10^−02	63.304	1.374	6.460 × 10^−02
4.428 × 10^−05	1.714 × 10^−02	63.203	1.376	6.637 × 10^−02
4.405 × 10^−05	1.754 × 10^−02	62.562	1.374	6.691 × 10^−02
4.395 × 10^−05	1.793 × 10^−02	63.140	1.374	6.757 × 10^−02
4.387 × 10^−05	1.817 × 10^−02	63.122	1.373	6.791 × 10^−02
4.355 × 10^−05	1.865 × 10^−02	63.209	1.366	6.801 × 10^−02
4.339 × 10^−05	1.903 × 10^−02	63.246	1.362	6.824 × 10^−02
4.329 × 10^−05	1.994 × 10^−02	63.168	1.361	6.931 × 10^−02
4.323 × 10^−05	2.026 × 10^−02	58.926	1.362	6.978 × 10^−02
4.317 × 10^−05	2.099 × 10^−02	63.162	1.356	7.015 × 10^−02
4.310 × 10^−05	2.147 × 10^−02	63.189	1.353	7.034 × 10^−02

FIGURE 8 Pareto set for the corn drying process.

Simulations of the model using the best set of results for each objective function are illustrated in Figs. 9and 10. The influence of the multi-objective behavior can be observed in both cases by the difference on the corn moisture content profiles.

FIGURE 9 Simulation of the moisture content profile as a function of time using the best set of results for F ₁ and F ₂.

FIGURE 10 Simulation of the moisture content profile at different heights using the best set of results of F ₁ and F ₂.

In Fig. 9, the parameter estimation based on the best result of the objective function F ₁ provides the best description of the moisture content variation as a function of the time. In Fig. 10, the moisture content profile at different heights is better simulated using the best set of results for F ₂. Similar behaviors can be detected in the air temperature profiles illustrated in Figs. 11 and 12, but in a slight way. The selection of the best compromise solution will depend on the judgment of the decision maker, based on additional knowledge of the system.

FIGURE 11 Simulation of the air temperature profile as a function of time using the best set of results for F ₁ and F ₂.

FIGURE 12 Simulation of the air temperature profile at different heights using the best set of results of F ₁ and F ₂.

Finally, for the sake of comparison, two different correlations proposed in the literature for this process by Mancini [16] and Lewis [19] are presented. Figures 13–16 show the simulation of the process using the estimated parameters and the two correlations, for operation times of two and six hours.

FIGURE 13 Variation of the corn moisture content during the drying process (z = 50 cm).

FIGURE 14 Corn moisture profiles at different bed heights (t = 120 min).

FIGURE 15 Air temperature profiles during the drying process (z = 50 cm).

FIGURE 16 Air temperature profiles at different bed heights (t = 120 min).

It can be observed that the estimated results are the ones that best fit the experimental data in all tested cases.

Conclusions

Our focus in this contribution has been to evaluate the performance of an evolutionary algorithm in conducting multi-objective parameter estimation. The algorithm has been applied to a grain cooling process and a grain drying process in order to determine the heat and moisture transfer coefficients.

The model was validated against the experimental data. Simulations were carried out to verify the estimated parameters. The model predictions have been shown to match the measured values within experimental error.

Nomenclature

Table

Cp	=	specific heat (kJ/kgK or cal/g °C)
Dp	=	particle diameter (m or cm)
fm	=	drying rate
F	=	objective function
F k	=	fitness function
G	=	mass velocity of drying gas (kg/m^2s or gs/cm^2)
h	=	heat transfer coefficient (kJ/m^3s or cal/cm^2 s °C)
H	=	specific enthalpy (kJ/kg)
N	=	number of experimental data
S	=	solid superficial area/bed volume (m^−1)
t	=	time (s)
T	=	temperature (K)
UR	=	relative humidity
x	=	gas flow direction
X	=	dimensionless parameter
Y	=	moisture content of the grain (−)
z	=	bed height coordinate
Z	=	dimensionless parameter

Greek Symbols

Table

ε	=	void fraction in bed of solids
ϕ	=	particle shape factor
λ	=	latent heat of evaporation (kJ/kg)
ρ	=	density (kg/m^3 or g/cm^3)

Subscripts and Superscripts

Table

exp	=	experimental
g	=	gas
gt	=	total gas
l	=	liquid
R	=	reference
s	=	solid
se	=	solid in equilibrium
v	=	vapor

Acknowledgments

We would like to thank the CNPq – Conselho Nacional de Desenvolvimento Científico e Tecnológico – for the scholarship provided.

Additional information

Notes on contributors

Evaristo C. Biscaia

†E-mail: [email protected]

Notes

^†E-mail: [email protected]