Multi-objective Dynamical System Parameters and Initial Value

Identification Approach in Chemical Disintegration Reaction

Modelling

Ivan Ryzhikov, Christina Brester and Eugene Semenkin

Institute of Computer Sciences and Telecommunication, Siberian State Aerospace University, Krasnoyarsk, Russia

Keywords: Multi-output System, Linear Differential Equation, Multi-objective Optimization, Parameters Identification,

Initial Value Estimation.

Abstract: A multi-criteria multi-output dynamical system identification problem is considered. The inverse

mathematical problem of estimating the parameters of a system of differential equations and its initial point

using the measured data is provided for the hexadecane disintegration reaction. The aim of modelling is to

approximate the dynamical behaviour of hexadecane and the concentrations of its products, which

according to chemical kinetics are determined by a differential equation. Since the dynamical model

adequacy is based on the error between its output and the sample data and the output itself depends on the

initial point values, the inverse mathematical modelling problem is the simultaneous estimation of the model

parameters and the initial point. At the same time, the initial point is unknown and the sample data is noisy,

and for this reason, the inverse mathematical modelling problem is reduced to a two-objective optimization

problem. The reduced problem is a sample of black-box optimization problems; it is complex, multimodal

and requires a reliable technique to solve it. This is why a specific heterogeneous multi-objective genetic

algorithm with the island meta-heuristic is used and its efficiency in solving this problem is proved by the

investigation results.

1 INTRODUCTION

The inverse mathematical modelling problem for

dynamical systems is a problem of current

importance. It has many applications in different

scientific fields, including chemistry, robotics,

mechanics, information science, biology and

financial mathematics. There are many different

problem statements and most of them are solved by

finding the solution of the reduced optimization

problem. In general, these optimization problems are

so-called black-box optimization problems and, in

this case, the model determination using a set of

parameters and the extremum seeking algorithm are

the most important objects to develop and improve.

In this study, we consider a multi-output system

without any control input, but with the vector on the

right-hand side for the process of hexadecane

disintegration. The proposed approach is applicable

to a linear dynamical system with many control

inputs, but in this study, the investigated systems

have no control inputs and so these inputs were set

as unit-step functions. Our aim is to build a

mathematical model of the hexadecane and the

functions of its disintegration product

concentrations, which would give the information

about the amount of acids at different time points in

the process. Two different process types were

considered: diffusion and static, and for each type

different disintegration reaction products were

considered.

The problem of estimating the parameters of the

differential equation using evolutionary, stochastic

and other heuristic approaches can be found in many

studies. Commonly, evolution-based and nature-

inspired algorithms are reliable solvers of reduced

optimization problems of different natures. Single-

objective and multi-objective optimization problems

appear in estimating the parameters of dynamical

systems in many different studies. Methods of

solving the parameter identification problem with

heuristic algorithms of these classes can be found in

many different studies.

There is a significant difference between solving

one-criterion and multi-objective optimization

Ryzhikov, I., Brester, C. and Semenkin, E.

Multi-objective Dynamical System Parameters and Initial Value Identiﬁcation Approach in Chemical Disintegration Reaction Modelling.

DOI: 10.5220/0006431504970504

In Proceedings of the 14th International Conference on Informatics in Control, Automation and Robotics (ICINCO 2017) - Volume 1, pages 497-504

ISBN: 978-989-758-263-9

497

problems and a significant difference in the solution

of the optimization problem. For the multi-objective

optimization problem, an approximation of the

Pareto front is required, or, in other words, the set of

non-dominated solutions. The reason why we need

to reduce the inverse mathematical modelling

problem to a multi-objective optimization problem is

the impossibility of estimating the system

parameters without its initial value for dynamical

systems.

The current study is based on the works

described in (Ryzhikov and Semenkin, 2017) and

(Ryzhikov et al., 2016), where the initial point and

the system coefficients were estimated by solving

the optimization problem, but for a single criterion.

The approach proposed in this study is an extension

of previous ones and allows different initial value

estimations to be used. Because of the sample size,

which is small, the initial value was compared with

the first measurements in the sample, but in a

general case, it can be compared with some average

value, its statistical estimation or some expected

value given by an expert, which makes the proposed

approach flexible and useful.

It is widely known that each particular problem

requires specific problem-oriented optimization

algorithm modifications or meta-heuristics, which

provide an improvement in the algorithm

performance. In this study, we propose a dynamical

system and initial value vector determination, and

form an unconstrained extremum problem on a real

value vector field.

In our consideration of the multi-objective

optimization problem, there is a need for a reliable

multi-objective optimization tool. The multi-criteria

cooperative heterogeneous genetic algorithm with

the island meta-heuristic has been chosen, because

this algorithm has proved its performance in the

solving of many multi-criteria optimization

problems of different natures. Experimental results,

which are presented in various figures, show that

this algorithm is applicable to the considered

problem and achieves a very good result.

2 MULTI-OUTPUT DYNAMICAL

SYSTEM IDENTIFICATION

PROBLEM

In this chapter, we consider the problem statement

and related real problem. The hexadecane

disintegration reaction consists in the forming of

different components: spirits, carbonyl components,

lactones and acids for the diffusion reaction and

spirits and carbonyl compounds for the static

reaction. We assume that the concentrations of

products and the hexadecane itself are related and

affect one another. The influence of each

disintegration reaction component concentration on

other component concentrations can be

mathematically determined by the addition of some

coefficient to the particular equation. The process of

how the concentration changes is dynamical and

linear, and in this case it is also assumed in chemical

kinetics theory that the concentrations can be

determined with a linear differential equation of the

first order.

Let a set





,, 1,

Yt i s , be a data sample, where

YR is the dynamical system output

measurements at the time point

t ,

is the size of

the sample and

is the system output dimension.

As was mentioned above, the solution of the inverse

mathematical problem is a differential equation and

the initial point, so we assume that the real system

output

















...

txt xt is determined by

the Cauchy problem:













11,11 1, 1nn

taxt axtb





 ,

…,













,1 1 ,nn nnnn

taxt axtb













0...

xx .

(1)

where

,1,

aij n are the system coefficients

and





0x is the initial point. The system parameters

and the initial point coordinates, generally, are

unknown and to be estimated.

We use the following notation for the

concentration functions of the chemical reaction

products: hexadecane, spirits, carbonyl compounds,

lactones and acids and similar notation for the static

chemical reaction.

We also know that the output data

,1,

Yi s is

distorted by the additive noise

:()0,()ED



 :

() , 1,

iii

YXt i s.

(2)

where the





t function is a solution of the Cauchy

problem (1).

Let the model be determined by the given

problem: a differential equation in a matrix form and

the initial value vector (we know the number of

system outputs n )

ICINCO 2017 - 14th International Conference on Informatics in Control, Automation and Robotics

498



ˆˆ

Xt B

 ,

















ˆˆ

00...0

Xx x

(3)

where

is a system matrix, B

is a vector of

coefficients and

)0(

X is an initial point.

It can be seen that the solving of the inverse

mathematical modelling problem is related to the

parameters and initial point estimation problems on

the basis of the sample data, which consists of noisy

output measurements. Moreover, the sample size for

some particular problems is small and the sample

data is flat. These data characteristics make many

estimation approaches useless, and require specific

tools, which would provide the simultaneous

estimation of the parameters and the initial point.

Since there is a cross influence between the

parameters and the initial point, but at the same time

the values of these variables are estimated and based

on the sample data, we require not just a criterion,

which minimizes the error between the data and the

model output, but also a number of criteria.

One of the proposed criteria is the distance

between the sample data and the model output

 



ˆˆ

,, 0

min

iji

AB X

YXt









(4)

where





t is the solution of the Cauchy problem

(3) for matrix

, vector

and ,1,

Dj n is a

proposed norming value



max min

ji i

DY Y.

The second criterion is the distance between the

initial point estimation and the sample



0min

CYX  .

(5)

The concentrations according to the physical

meaning of this term cannot be negative; to achieve

this, we add the penalty function into the criterion

(4) and criterion (5) with some penalty coefficients.

The penalty function is defined by the expression



0, 0











(6)

Finally, we receive the following criteria for the

multi-objective optimization problem (4), (5) and (6)



ˆˆ

,, 0

min

Pji

AB X

CCc PXt



 



(7)



0min

C С cPX



 



(8)

where

cc

are the penalty function

coefficients. In the current study, these coefficients

were equal to

, 1000

cc The current model

determination leads to an optimization problem in

the

nn n n

RRR



 field. By determining the problem

in this way, the model can be represented with a

vector from this field. As can be seen, the problem

dimension adds extra complexity to the

identification problem.

In this study, the Cauchy problem (3) is solved

numerically using the Runge-Kutta integration

scheme, which makes the designed approach flexible

and useful for nonlinear or nonstationary systems

with a known symbolic form but unknown

parameters or different input functions.

Now we consider the multi-objective

optimization tool, which has been designed to find a

solution for complex black-box optimization

problems. The proposed algorithm consists of some

stochastic and evolution search operators and meta-

heuristics. The fitness function is performed on the

basis of criteria



ˆˆ

1arg

fit a

Ca C







ˆˆ

1arg

fit a

Ca C





(9)

where

arg( )

I

is a specific transformation of an

alternative

to the arguments of a criterion

. This

transformation decodes the individual vector

coordinates into matrix coefficients, right-side

vector coefficients and into the initial point in the

case of the

C criterion. For the

C criterion, it

converts the solution vector into the initial point.

3 MULTI-OBJECTIVE

COOPERATIVE GENETIC

ALGORITHM WITH THE

ISLAND META-HEURISTIC

While solving multi-objective optimization

problems, we tend to reach a trade-off between

competing criteria. The Pareto-dominance idea

(Goldberg, 1989) underlies the way we may

Multi-objective Dynamical System Parameters and Initial Value Identiﬁcation Approach in Chemical Disintegration Reaction Modelling

499

compare alternative solutions. As a result, we

expect to obtain a set of non-dominated points which

cannot be preferred to one another based on all the

criteria considered.

Population-based algorithms (in particular,

genetic algorithms) operate with a number of

candidate-solutions at each generation, and

therefore, it was decided to use them as an effective

tool to find Pareto set and front approximations.

However, there are some open questions researchers

usually encounter when they apply multi-objective

genetic algorithms (MOGAs) in real problems.

Whether choosing one of the existing MOGAs or

designing a new one, researches should elaborate

three main concepts which are incorporated into the

scheme of any MOGA and opt for the most

appropriate implementation of each concept for the

problem being solved. Firstly, various fitness

assignment strategies might be proposed (Zitzler,

2004): the dominance rank (the number of points by

which the candidate-solution is dominated), the

dominance depth (a population is divided into

several fronts or niches and it is determined which

front an individual belongs to), or the dominance

count (the number of points dominated by the

candidate-solution) might be used to assign a fitness

function.

Secondly, to keep variety within the Pareto set

and front approximations, diversity preservation

techniques are applied. In (Silverman, 1986) several

variants of these techniques are presented: kernel

methods assess the density with a Kernel function

which takes the distance to another point as an

argument; nearest neighbour techniques are based on

estimating the distance between a given point and its

k-th nearest neighbour; and histograms, using a

hypergrid to calculate neighbourhoods, relate to

another class of density estimators. In most cases,

these approaches calculate the distance between

points in the objective space.

Moreover, to avoid the loss of effective

individuals during the algorithm execution due to

stochastic effects, the idea of elitism has been

suggested. There are two basic ways to implement it.

The first way is to merge the parent population with

the offspring and then to employ environmental

selection taking into account the fitness values of

individuals from the mating pool. Another variant is

based on the usage of an additional set called an

archive for copying promising solutions at each

generation.

Thinking about these issues, in this study we

decided to apply a cooperative MOGA (Brester et

al., 2015) which combines three algorithms based on

different heuristics. The use of the cooperative

MOGA allows us to eliminate the choice of the most

effective algorithm and avoid multiple experiments

with many different MOGAs. The cooperative

MOGA uses an island model (Whitley et al., 1997)

and includes NSGA-II (Non-Sorting Genetic

Algorithm II) (Deb et al., 2002), PICEA-g

(Preference-Inspired Co-Evolutionary Algorithm

with goal vectors) (Wang, 2013), and SPEA2

(Strength Pareto Evolutionary Algorithm 2) (Zitzler

et al., 2002) as its islands work in a parallel way

(Figure 1).

Figure 1: The Island Model Implemented.

The initial number of individuals M is spread

across L subpopulations: M

=M/L, i=1,…,L and the

same number L of threads is initialized. Thus, the

fitness function evaluation for different

subpopulations is implemented in parallel threads.

At each T-th generation algorithms exchange the

best solutions (migration). There are two parameters:

migration size, the number of candidates for

migration, and migration interval, the number of

generations between migrations. Moreover, it is

necessary to define the island model topology, in

other words, the scheme of migration. The fully

connected topology is used, meaning that each

algorithm shares its best solutions with all other

algorithms included in the island model. The multi-

agent model is expected to preserve a higher level of

genetic diversity. The benefits of the particular

algorithm are advantageous in different stages of

optimization.

The cooperative MOGA was investigated on the

set of complex benchmark problems CEC 2009

(Zhang, 2008) and proved its effectiveness. It also

was applied in a wide range of practical problems:

emotion recognition from speech (for feature

selection and neural-network design) (Brester et al.,

2016), the prediction of cardiovascular diseases (for

feature selection) (Brester et al., 2016), and

spacecraft control (for the choice of control contour

variant by solving an optimization problem)

(Semenkina et al., 2016).

ICINCO 2017 - 14th International Conference on Informatics in Control, Automation and Robotics

500

4 HEXADECANE

DISINTEGRATION REACTION

PRODUCT CONCENTRATION

MODELLING

To solve the problem we chose the following

algorithm resources: 3 different populations with

200 individuals and 500 generations, the migration

size was equal to 25, the migration interval was 50.

The number 300000 limited the total amount of the

fitness function evaluations. After 25 algorithm runs

we selected some solutions to be demonstrated.

From the whole Pareto front estimation, we selected

the model with the highest first criterion value

second criterion value

and the one model with

the values in between

. All the models are

presented in the same Figures: the first model is a

long dashed curve, the second is a short dashed

curve and the last one is a dotted curve. The

measurements are marked with grey crosses on the

plot.

Figure 2: Hexadecane concentration: model outputs and

sample measurements.

Figure 3: Spirit concentration: model outputs and sample

measurements.

Figure 4: Concentration of carbonyl compounds: model

outputs and sample measurements.

Figure 5: Lactone concentration: model outputs and

sample measurements.

In Figure 2, Figure 3, Figure 4, Figure 5 and

Figure 6 the hexadecane concentration, spirits,

carbonyl compounds, lactones and acids are

presented, respectively.

If we compare the initial value and the measured

initial value, it can be seen that the model that fits

the data sample better has the highest inaccuracy in

the initial value. At the same time, the model with

the highest value of the second criterion does not fit

the sample data as well as other models.

Figure 6: Acid concentration model outputs and their

measurements.

X ,

Multi-objective Dynamical System Parameters and Initial Value Identiﬁcation Approach in Chemical Disintegration Reaction Modelling

501

In Figure 7 the Pareto front estimation obtained

by multiple algorithm runs is given. In the same

Figure, the solution found in one algorithm run is

marked with grey points.

Figure 7: Pareto front estimations for each run (black)

and a single run Pareto front estimation (grey).

We see that the Pareto front estimation is a very

difficult problem, which is why there is a very

complicated relation between the closeness to the

initial point estimation and how well it fits the

sample data. This means that there is a necessity to

solve the inverse mathematical problem as an

optimization problem with two criteria.

The same problem was solved for another

chemical experiment with different reaction

characteristics. In this case it was necessary to build

a model for only three outputs: hexadecane

concentration, spirit concentration, and carbonyl

compound concentration. For this problem we also

took three different solutions from the Pareto front

estimation in a similar way to how it was performed

above.

Figure 8: Hexadecane concentration: model outputs and

sample measurements.

Figure 9: Spirit concentration: model outputs and sample

measurements.

Figure 10: Concentration of carbonyl compounds: model

outputs and sample measurements.

As can be seen for the second problem, there is

an abnormal measurement. However, the model

“ignores” this measurement and fits the sample data.

We can conclude that the second problem is easier,

since the models do not differ as much as they did in

the case above. Also, the Pareto front estimation

proves this hypothesis; the estimation is given in

Figure 11.

Figure 11: Pareto front estimations for each run (black)

and a single run Pareto front estimation (grey).

X ,

ICINCO 2017 - 14th International Conference on Informatics in Control, Automation and Robotics

502

The last Figure gives us the same information as

Figure 7 but for the current problem. It is easier to

estimate the Pareto front since many points are

localized near the same curve.

4 CONCLUSION

The experimental results of this study prove that the

proposed approach is useful in solving inverse

mathematical problems for dynamical systems in

cases when the initial point and the noisy sample

data are unknown. Using this approach, many

models of hexadecane disintegration reaction

product concentrations were build. It was

demonstrated that these models fit the observation

data well and behave normally.

In this paper, the multi-output dynamical system

identification problem was solved by means of the

multi-objective heterogeneous genetic algorithm

with the island meta-heuristic. The results prove the

high efficiency of the algorithm used and the

applicability of the proposed approach, which allows

us to solve the inverse mathematical modelling

problem in the case of having no information about

the initial point value and satisfy the trajectory

constraints.

It can be seen that the model output fits the

sample data well and represents the physical

properties of the process. The multi-objective

problem reduction allows us to receive the Pareto

front estimation on the basis of estimations of the

initial point and system parameters, so the expert can

vary the degree of belief in the initial point values

and choose the mathematical model that would

satisfy his modelling needs. Moreover, the proposed

two-criterion approach allows mathematical models

to be found, the parameters and initial value

characteristics of which can contradict. As can be

seen in Figures 7 and 11, different problems have

different Pareto fronts, but the criteria have a

complex relation and so they cannot be represented

as a single one.

The considered sample data has a small size,

which makes it impossible to apply statistical

methods for the initial value estimation or apply

some other identification techniques based on

approximating the model output as a static function.

This is the reason why the differential equation

based models are the most important part of

modelling the dynamical processes and so it is

important to develop the algorithms for the equation

parameter identification.

Further work is related to the inverse

mathematical problem solving for multi-output

dynamical systems of higher order and control

inputs. Another direction is the developing of

heuristic optimization tools for the single and multi-

criteria problems of dynamical system identification,

and designing problem-oriented modifications.

ACKNOWLEDGEMENTS

This research is supported by the Russian

Foundation for Basic Research within project No 16-

01-00767.

REFERENCES

Brester, C., Kauhanen, J., Tuomainen, T.P., Semenkin, E.,

Kolehmainen, M., 2016. Comparison of two-criterion

evolutionary filtering techniques in cardiovascular

predictive modelling. Proceedings of the 13th

International Conference on Informatics in Control,

Automation and Robotics (ICINCO’2016), Lisbon,

Portugal, vol. 1: pp. 140–145.

Brester, C., Semenkin, E., 2015. Cooperative

Multiobjective Genetic Algorithm with Parallel

Implementation. Advances in Swarm and

Computational Intelligence, LNCS 9140: pp. 471–

478.

Brester, C., Semenkin, E., Sidorov, M., 2016. Multi-

objective heuristic feature selection for speech-based

multilingual emotion recognition. Journal of Artificial

Intelligence and Soft Computing Research, vol. 6, no.

4: pp. 243-253.

Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A

fast and elitist multiobjective genetic algorithm:

NSGA-II. IEEE Transactions on Evolutionary

Computation 6 (2): pp. 182-197.

Goldberg, D., 1989. Genetic algorithms in search,

optimization, and machine learning. Addison-wesley.

Ryzhikov, I., Semenkin, E., 2017. Restart Operator Meta-

heuristics for a Problem-Oriented Evolutionary

Strategies Algorithm in Inverse Mathematical MISO

Modelling Problem Solving. IOP Conference Series:

Materials Science and Engineering, vol. 173.

Ryzhikov, I., Semenkin, E., Panfilov, I., 2016.

Evolutionary Optimization Algorithms for Differential

Equation Parameters, Initial Value and Order

Identification. ICINCO (1) 2016: pp. 168-176.

Semenkina, M., Akhmedova, Sh., Brester, C., Semenkin,

E., 2016. Choice of spacecraft control contour variant

with self-configuring stochastic algorithms of multi-

criteria optimization. Proceedings of the 13th

International Conference on Informatics in Control,

Automation and Robotics (ICINCO’2016), Lisbon,

Portugal, vol. 1: pp. 281–286.

Multi-objective Dynamical System Parameters and Initial Value Identiﬁcation Approach in Chemical Disintegration Reaction Modelling

503

Silverman, B., 1986. Density estimation for statistics and

data analysis. Chapman and Hall, London.

Wang, R., 2013. Preference-Inspired Co-evolutionary

Algorithms. A thesis submitted in partial fulfillment

for the degree of the Doctor of Philosophy, University

of Sheffield: p. 231.

Whitley, D., Rana, S., and Heckendorn, R., 1997. Island

model genetic algorithms and linearly separable

problems. Proceedings of AISB Workshop on

Evolutionary Computation, vol.1305 of LNCS: pp.

109-125.

Zhang, Q., Zhou, A., Zhao, S., Suganthan, P. N., Liu, W.,

Tiwari, S., 2008. Multi-objective optimization test

instances for the CEC 2009 special session and

competition. University of Essex and Nanyang

Technological University, Tech. Rep. CES–487, 2008.

Zitzler, E., Laumanns, M., Bleuler, S., 2004. A Tutorial on

Evolutionary Multiobjective Optimization. In:

Gandibleux X., (eds.): Metaheuristics for

Multiobjective Optimisation. Lecture Notes in

Economics and Mathematical Systems, Vol. 535,

Springer.

Zitzler, E., Laumanns, M., Thiele, L., 2002. SPEA2:

Improving the Strength Pareto Evolutionary Algorithm

for Multiobjective Optimization. Evolutionary

Methods for Design Optimisation and Control with

Application to Industrial Problems EUROGEN 2001

3242 (103): pp. 95-100.

ICINCO 2017 - 14th International Conference on Informatics in Control, Automation and Robotics

504