Automatic linear differential equation identification in analytical form Текст научной статьи по специальности «Математика»

2. Mironovskiy S. N., Repetskiy O. V., Ryzhikov I. N.

Effekt rasstroyki lopatok rabochikh koles tyrbomashin (effect of mistuning blades on the rotors of turbomachines). Vestnik stipendiatov DAAD, 2012, no. 2, p. 75-81.

3. Gisli Ottarsson, Mathew P. Castanier, Christophe Pierre. A reduced-order modeling technique for mistuned bladed disks. Journal AIAA, 1994-1640, p. 2552-2562.

4. Jan Ronnie Bladh. Efficient predictions of the vibratory response of mistuned bladed disks by reduced order modeling. Dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan, 2001.

5. Repetskiy O. V. Komp'juternyj analiz dinamiki i prochnosti turbomashin (Computer analysis of dynamics and strength turbomachines). Irkutsk, IrGTU Publ., 1999.

6. Zenkevich O. Me tod konechnyh jelementov v tehnike (The Finite Element Method in engineering science). Moscow, Mir Publ., 1975.

7. Bladh J. R., Castanier M. P., Pierre C. Component-Mode-Based Reduced order modeling techniques for Mistuned bladed Disks - Part I: Theoretical Models. Journal of Engineering for Gas turbines and Power, 2001, no. 8, p. 89-99.

8. Boriskin O. F. Kulibaba V. V. Repetskiy O. V. Konechnojelementnyj analiz kolebaniy mashin (Finite element analysis of vibration machines). Irkutsk, IrGTU Publ., 1989.

9. Hurty W. C. Dynamic analysis of structural systems using component modes. Journal AIAA, 1965, no. 4, p. 678-685.

© Репецкий О. В., До Мань Тунг, 2014

УДК 005; 519.7; 303.732

AUTOMATIC LINEAR DIFFERENTIAL EQUATION IDENTIFICATION IN ANALYTICAL FORM*

I. S. Ryzhikov

Siberian State Aerospace University named after academician M. F. Reshetnev 31, Krasnoyarsky Rabochy Av., Krasnoyarsk, 660014, Russian Federation Е-mail: ryzhikov-88@yandex.ru

In this study a multi-agent evolutionary strategies algorithms system is considered in application to linear dynamic system identification problem solving. The proposed approach is based on previous results of designing the universal simultaneous parameters and structure identification technique that uses the observations of system output and input. Every agent is model-based and has an aim to find an extremum for agent's criterion. The solution of the problem is symbolic and can be automatically found for the small samples and noised data.

Keywords: evolutionary strategies, multi-agent system, identification, structure and parameters, differential equation.

АВТОМАТИЧЕСКАЯ ИДЕНТИФИКАЦИЯ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

В АНАЛИТИЧЕСКОМ ВИДЕ

И. С. Рыжиков

Сибирский государственный аэрокосмический университет имени академика М. Ф. Решетнева Российская Федерация, 660014, г. Красноярск, просп. им. газ. «Красноярский рабочий», 31

Е-mail: ryzhikov-88@yandex.ru

Рассматривается применение многоагентной системы алгоритмов эволюционных стратегий для решения задач идентификации линейных динамических систем. Предлагаемый подход основан на результатах предыдущих работ по разработке и совершенствованию универсального метода, позволяющего автоматически решать задачу идентификации одновременно параметров и структуры объекта по данным наблюдений его входа и выхода. Каждый агент основан на модели и имеет цель найти оптимум функционала; функционалы отличаются для всех агентов. Решение ищется в аналитическом виде автоматически, в том числе и для случаев малой выборки.

Ключевые слова: эволюционные стратегии, многоагентные системы, идентификация, структура и параметры, дифференциальное уравнение.

* The study was supported by The Ministry of education and science of Russian Federation, project № 14.B37.21.1521. The second International Workshop on Mathematical Models and its Applications (IWMMA 2013).

The problem of making a model for a dynamic plant or process is not new and has a lot of definitions. For some problem definitions there are quite a lot of techniques, so problems can be called classic today. For instance, there are techniques that are based on the mathematical statistics maximum likelihood method [1] or techniques that are based on autocorrelation function and nonparametric kernel statistics [2]. Still we can face different systems or definitions that require new approaches to reach a success, since the nature of the task makes it useless or impossible to apply common techniques.

One of these problems is so-called identification of the black box, when the structure is unknown and there are only observations of inputs and outputs. It is still possible to make a model in a form of stochastic differential equation [3], if the order of the equation is known and there is an acceptable observation of the reaction on the unit-step function. Also, it is possible to estimate the transient function via fuzzy output estimation [4], artificial neural network modelling [5; 6], genetic programming [7], or another function estimation technique. Every of suggested methods have their own advantages and disadvantages. But all of them allow to find or the estimation for the observed time segment only, or the estimation, which symbolic representation makes it useless for analytic methods of control or analysis. Some of them are dependent on the size of the sample. In this study we consider a special case of the problem, when the observations are noised, the size of the sample can be small and the aim is not just to estimate the dynamic model of the system but estimate it in symbolic form.

In the article [8] the linear dynamic system identification with the 2nd order ordinary differential equation via the genetic algorithm was examined. Although, the genetic algorithm is a well-known and reliable global optimization technique, it is more preferable to apply the seeking not only on the given compact and avoid quantizing of the arguments. That is why it was suggested to implement the optimization technique for the real valued arguments.

In the previous work the approach of reducing the initial identification problem to extremum problem on the real vector field was suggested and the modified and hybrid evolutionary strategies algorithm was designed and examined. The efficient algorithm settings were found to be promising in application to the identification task.

Linear differential equations models are useful in many different fields, so the linear dynamics identification problem is significant and takes place in different problem definitions. The filtering problem and articulatory identification [7; 9], are related to the stochastic ODE identification. The proposed approach can be extended to stochastic ODE or Bessel equations identification. Hence, it is also can be applicable to Markov processes [10; 11]. That is why the linear differential equation identification can be useful is some fields related to speech recognition problem, gesture recognition problem and, probably, many more.

Let one have the sample of one input, one output

system measurements {yi, ui, ti}, i = 1, s , where s is its

size, yi e R is dynamic system output measurement at

time ti, and ui = u(ti) is control measurement. It is also

known, that the dynamic system is linear and stationary, so it can be described with the ordinary differential equation (ODE):

ak • x(k) + ak-1 • x(k-1) +... + a0 • x = b • u(t) , x(0) = x0 . (1)

Here x0 is supposed to be known. In the case of the

transition observation, we can put forward a hypothesis about initial point: the system output is known at time t = 0 and the derivative values can be set to zero, because usually the system observation starts with it's being in the steady state. In general, the initial point can be approximated or even being the part of optimization problem. Using the sample data we need to identify parameters and the system order m , which is assumed to be limited, so m < M, M e N . M is a parameter that is set by the researcher. This value limits the structure of the differential equation, i. e., it limits the ODE order. It is also assumed that there is an additive noise Q: E(Q) = 0, D(Q) , that affects the output measurements:

yt = x(tt) + Q . (2)

Since the identification problem definition requires the model in symbolic form, the system structure is needed. For the ODE structure can be defined with the order of the differential equation and its parameters. That is why we put forward a hypothesis about ODE order, that it is less than chosen maximum order for the equation M . So now, the proposed approach is uderparametrized.

Without loss of the generality, let the leading coefficient of ODE be the constant equal to 1, so that

x( k)

+ • x(k-1) +... + ^ • x = b • u(t),

or

v( k ).

+ ak • x(k 1)1 +... + <5j • x = b • u(t).

(3)

(4)

Then one can be searching for the solution of the identification task as a vector of parameters for linear differential equation with the order m

?Cm) .

xv"; + aim • x(m u +... + ail • x = a0 • u(t) , x(0) = x0 , (5) where the vector of equation parameters a = (0, ..., 0, aim, ..., (5j, a0)T eRn, n <M +1, delivers an extremum to the functional

I^n) = £|y; - x(ti)|

^ min . (6)

aeRn, n<M+1

In general case, the solution x(t) is evaluated with a numerical integration technique, because the control function has no analytical from, rather is given algorithmically. We prefer the criterion (6) instead of quadratic criteria because of its robustness.

i=1

Previous study was focused on designing a special technique that allows automatically estimate the order and the parameters with an algorithm, which argument is a vector with fixed length. The advantage of that technique is an ability of simultaneous problems solving. Anyway, the main disadvantage was the fact that proposed reduction required a special scheme for integration of ODE. And since the order reduction is connected with setting the first vector's arguments to 0, the optimization algorithm should deal with this fact. Actually, it means that for stochastic optimization algorithms on the real vector field some special modifications are needed to be designed and examined, because of the solution sensitiveness to small disturbances that takes place.

The current study is based on another principle of order and parameters estimation, but still the algorithm provides automatic simultaneous seeking. Let every optimization problem for different values of variable n be solving with a distinct algorithm.

As a basis of the optimization technique the evolutionary strategies algorithm was chosen. The basic idea of this algorithm is described in [12]. This heuristic evolution-based technique was chosen, because the identification problem leads to solving the multimodal optimization task with complex criterion, that can be evaluated only numerically. The system structure and its parameters are defined with an integer and a vector. The criteria (6) for these vectors are complex and sensitive to its components, which are changing by stochastic search operators. This is why we have to develop the specific modification for the global optimization technique. The proposed approach is based on hybridized and modified evolutionary strategies algorithm (HMES).

Let every individual be represented with tuple

Hi = (op1, sp1, fitness(op1 , i = 1, Nj , where opi e R, j = 1, k is the set of objective parameters of the differential equation; spj e R+, j = 1, k is the set of

snoffspring = srpffspring + z APi - PPi "Tz

• N (0,1)|;

strategic parameters; N j is the

fitness(x): Rk ^ (0,1], fitness(x) = -

population size; 1

is

the

1 +1 (x) fitness function.

Selection types were borrowed from the genetic algorithm: proportional, rank-based and tournament-based. The algorithms in this study are based on produces one offspring from two parents and every next population have the same size as previous. Recombination types initially were chosen from the intermediate, weighed, randomly weighed and discrete crossovers. The mutation of every offspring's gene happens with the chosen probability pm . If we have the random value z = {0,1}, P(z = 1) = pm , which is generated for every current objective gene and its strategic parameter then

= opoffsrring ■

:• N (0, spfspnng ),

where N (m, a ) is the normally distributed random value

with the expected value m and the variance a2. The other way of implementing the mutation operand is with using the exponential adaptation of strategic parameters with an extra parameter of adaptation T :

op°ifspring = opfspring + z . N(0,spfsprmg ),

pffsp^g = pjfipmg el-z-N (0,1)

spi — ¿pi c

The stochastic evolutionary algorithm was hybridized with stochastic coordinate-wise optimization technique. Let N1 randomly chosen individuals for N 2 randomly chosen objective chromosomes are changing with N3 iterations of local search. The local optimization step is hl . After every circle and iteration the initial point is

being replaced with best solution. If all others solutions are worse, it does not change.

During the development of the ODE identification system that was based on modified evolutionary strategies algorithms parallelization some investigations were done. First, the sample for the ODE with chosen order was generated, and then every client program was searching for the best parameters of model. Every client was working with its own order. After the system stops evaluations, it brings as many models with different orders of ODE, as many clients were launched. Here we present some results of relation between the maximum fitness function values for the best solution found and the order of the system. On the fig. 1 diagrams are demonstrating this relation for the samples generated by systems with different orders of the ODE: a - second order; b - fourth order, c - sixth order, d - eighth order. So, the every sample was identified with ODE, orders of models varied from 1st to 9th. The algorithm is stochastic, that is why the results were averaged for 20 runs in every case. The settings of algorithms were examined, so every client had the same and the most efficient settings due to efficiency estimation done.

On these figures one can see that for the every system, the model with the same order fits best than other models. Moreover, increasing the difference in orders between the model and the system leads to decreasing the fitness value. The relation has a form of unimodal mapping with the pike in point with correct order. That is why it was suggested to develop an algorithm that would reallocate resources for all HMES solving optimization problems related to ODE parameters identification for different orders of equation.

We skipped the varying of all the HMES settings in the current research. First of all, the aim of this study was to develop the identification technique that is based on a different approach for estimation of the order. We do believe that the paradigm would be successful with any settings that are efficient enough to reach the adequate estimation. Second of all, in previous studies [13; 14] and [15] a lot of settings investigations were done, also for the algorithm that is described above. Hence, the most

efficient settings were chosen. In a further work there would be another complete examination, since there is no doubt that with the new paradigm another algorithm's settings could improve the efficiency.

As it was mentioned earlier, every HMES is applied for the different ODE order. For instance, let the parameter be set: m = 9, so there would be 9 different

extremum problems with criterion (6) for n = 2,10 . Let the every HMES be the agent with its own aim - to find best solution on the given space and computational resources. Working in cooperation, these agents are trying to figure out, which one of them is the most efficient, and then to reallocate all the resources between agents. In the current work we consider a simple cooperation between the agents; reallocation and cooperation both can be called momentary and proportional-based. It is also possible to implement the dynamic reallocation and based on set of different characteristics.

Let besti, i = 1, m be the fitness of the best solution found by the i -th algorithm. Similarly, let sizet be the size of the population, N1i, N2i - parameters for local optimization, pm i - mutation probabilities, ki - number

of parameters. To prevent collapse of any agent, the minimum population size Min _ size > 0 is needed. So, if the sum of all the individuals in populations is

Sum _ size = X s

then the free resources are

Free = Sum _ size - Min _ size .

In the simplest case, when the agents are model-based reaction agents the reallocation can be implemented in following form:

sizei = round (fw (i, best, aw ) • Free) , where fw : N x R x Rw ^ [0;1] is increasing weight function and V best* e best*, besti e best : best* > besti ^ fw (i,best , aw) > fw (i,best, aw); aw - are the set of its parameters; round: N ^ R turns its argument to nearest integer. We considered these weight functions:

besti

fW(i, best) =-

- • Free,

2 e

fw (i, best, aw ) = —

X bestj

j=1

aw •bestj

- • Free, aw > 0 .

(7)

(8)

^ eaw ■besti j=1

Of course, using the best solution of every HMES is not the only way of how to design the scheme. It is also possible to use, for example, the average of fitness function.

To provide the investigation of the scheme testing samples were generated: 10 for every order from first to ninth. Parameters of the systems were randomly

generated: a'k = U(0,5), bk = U(0,5), i = 1,9, k = 1,i, where U (0,5) is uniformly distributed random value. The time of the process was set to 5. The control function was the step function, u (t) = 1. For sample forming s = 50 points were taken randomly. For every task 10 runs of the algorithm were executed. The goal of examination is to investigate the approach, determine its characteristics.

After examining the HMES algorithm in previous works, it was suggested to use the following settings: 50 individuals for 50 populations, local optimization parameters are N1 = 10, N2 = 10 and N3 = 1 with hj = 0,1, the tournament selection with the tournament size equals 2, the discrete crossover and the mutation with

the probability pm = —.

k

One can be interested in if there is resources reallocating does it need special regularization, since the different agents works with different dimension and so, requires more or less resources. In this work two different investigations were made. The first thing is that there is no relation between efficiency increasing and increasing the number of individuals or the number of local optimization steps. For every sample and every algorithm 10 more tests with 10 launches in every of them were

made: 1 - N1,j = 10 + j-10, j = _____, 2 - N{J = 50 +10 • j .

On the fig. 2, a, b similar diagram to one that was shown on fig. 1 is presented for the investigations of increasing local search resources and individuals, respectively. The sample was generated with the system of the 2nd order.

On the fig. 2, both 2, a and 2, b, there are different polylines; each is the average of the solution fitness function for the given experiment for different orders of model. The first interesting fact is that with increasing number of local search steps, we slightly improved the efficiency. In a test with varying the population size the variation of different runs for the same order have larger variance in a scew zone. Varying the local optimization steps leads the increasing of variance in a pike and shifts the best solution to the next to proper order solution more often than in test with population size.

Anyway, the statistics shows that varying both the number of steps and the size of population gives no sufficient improvements. Moreover, sometimes it even leads to small distortions in identifying the proper order. That is why, reallocating of the resources in the new heuristic algorithm is made due to the common sense.

Now let us consider the main strategy of reallocating and examine how the varying of the weight function parameters would react on its efficiency. The reallocating strategy in the current study was implemented as it is shown below:

- run the same number of circles every HMES;

- reallocate the resources for the number of individuals;

- set the number of individuals that are going to be improved to round (sizet • 0,2);

- check the stop criterion.

i=i

As the stop criterion it is possible to use the achievement of desired accuracy of the model or the exceeding of evaluations.

Here we investigate the influence of the parameter aw

value on the sizes of populations. First, we compare the distribution of the resources by the end of algorithm evaluation for different parameter values. For instance, the system of the 2nd order was chosen, the results are the average for 20 launches. The final sizes of populations in dependence of parameter's value are shown on fig. 3. Another diagram on fig. 4 is to show the dependence of final solution fitness for every order of the parameter value. The linear polyline is for weight function (7).

One can see that with increasing of the weight coefficient the pike of resources distributions are becoming

narrower. It means that there can be lack of resources to improve the solutions for other orders of ODE and in the same time there can be no improvement of solution that is already found. That solution, actually, is forming a pike.

On this figure once can easily see the problem that was noticed above. The large values of weight lead to a gap in identification efficiency for orders that are far from the pike. The next fig. 5 shows the character of the best individual fitness function value, average fitness for the population and population size. To fit everything on the same scales, we normed the population size, its max value was 178. And best fitness value were normed by its min and max values.

123456789 10 The order of the ODE model

1,1 1Я '

1Л л

« <U 1

.e — 0 9

S > 0 8

2 ,2 о 7

% = 0,6

£ 0,5

123456789 10 The order of the ODE model

123456789 10 The order of the ODE model

c

— —i

£ - 09

4- (5

<u >

во с

2 .2 0 7

Q) *

0.5

123456789 10 The order of the ODE model

d

Fig. 1. Fitness function of the best solution found average value in relation with the order of ODE model for different cases: samples generated by 2nd order ODEs (a); samples generated by 4th order ODEs (b); samples generated by 6th order ODEs (c); samples generated by 8th order ODEs (d)

123456789 10

123456789 10

The order of ODE

The order of ODE 6

Fig. 2. Varying the number of local optimization steps (a), and size of population (b), for different orders of model

а

0.9

0.9

0.8 -

0.8

а

Fig. 3. Distribution of populations' sizes in relation to weight parameter value

Fig. 4. The best solution fitness function value in relation to weight parameter value

Fig. 5. The characteristics of the population

All of the settings that were listed above and the weight function (8) with parameter aw = 19, which was determined to be the most reliable the proposed approach during the test demonstrated its high efficiency with less computational efforts as the HMES algorithm with rounding of its parameters for some problems. In the further work the aim of investigation would be to compare both algorithms and to design and examine the new approaches to reallocate resources.

The proposed approach can be used for automatic identification of the linear dynamic system by its observations. In this paper we presented a scheme of interaction between different optimization algorithms. Further work despite would it be the HMES with rounding or the co-evolution-like scheme is focusing also on implementing the cooperation of three different but dependent optimization problems: identification of the structure and parameters, identification of the input, identification of the starting point.

But first of all, after examining different ways of interaction between agents and reallocating techniques the approach we looking forward to implement identification system that would allow solving MIMO dynamic problems with unknown order of every system state.

References

1. Astreon K. J., Bohlin T. Numerical Identification of Linear Dynamic Systems from Normal Operating Records. Technical paper. IBM Nordic Laboratory, Lidingo, Sweden: 38 p., 1967.

2. Medvedev A. V. Identification and control for linear dynamic system of unknown order. Optimization Techniques IFIP Technical Conference. Berlin -Heidelberg - New York: Springer-Verlag, p. 48-56, 1975.

3. Zoteev V. Parametrical identification of linear dynamical system on the basis of stochastic difference equations. Matem. Mod. 2008. Vol. 20, no. 9, p. 120-128.

4. Passino K. M., Yurkovich S. Fuzzy control. Addison Wesley Longman, Inc. 1998. 502 p.

5. Tan Y., Dong R., Chien H., He H. Neural networks based identification of hysteresis in human meridian systems. Int. Journal of Applied Mathematics and Computer Science. 2012. Vol. 22, no. 3, p. 685-694.

6. Jovanovic O., Identification of Dynamic System Using Neural Network // Architecture and Civil Engineering. 1997. Vol. 1, no. 4, p. 525-532.

7. Cao H., Kang L., Chen Y., Yu J. Evolutionary modeling of systems of ordinary differential equations with genetic programming. Genetic Programming and Evolvable Machines1 (40) p. 309-337, 2000.

8. Parmar G., Prasad R., Mukherjee S. Order reduction of linear dynamic systems using stability equation method and GA. International Journal of computer and Inforna-tion Engeneering 1:1, 2007.

9. Reimer M., Rudzicz F. Identifying articulatory goals from kinematic data using principal differential analysis. Proceedings of Interspeech 2010, Makuhari Japan. 2010. P. 1608-1611.

10. Mineiro P., Movell J. R., Williams R. J. Modeling path distributions using partially observable diffusion networks: a Monte-Carlo approach. Technica lReport CogSci.UCSD-99.01, Department of Cognitive Science, UCSD, SanDiego, 1999.

11. Saerens M. Viterbi algorithm for acoustic vectors generated by a linear stochastic differential equation. Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Detroit, 1995, p. 233-236.

12. Schwefel Hans-Paul Evolution and Optimum Seeking : New York: Wiley & Sons., 1995.

13. Ryzhikov I., Semenkin E. Evolutionary Strategies Algorithm Based Approaches for the Linear Dynamic System Identification. Adaptive and Natural Computing Algorithms, Springer: p. 477-483, 2013.

14. Ryzhikov I., Semenkin E. Modified Evolutionary Strategies Algorithm in Linear Dynamic System Identification. ICINCO (1): p. 618-621, 2012

15. Ryzhikov I., Semenkin E. The Application of Evolutionary Algorithm for the Linear Dynamic System Modelling. SIMULTECH: p. 234-237, 2012.

© Pbi^HKOB H. C., 2014

УДК 629.19.22

ИССЛЕДОВАНИЯ ФЕРРОМАГНИТНЫХ ПОРОШКОВ И ПОКРЫТИЙ НА ИХ ОСНОВЕ

В. Н. Саунин, С. В. Телегин

Сибирский государственный аэрокосмический университет имени академика М. Ф. Решетнева Российская Федерация, 660014, г. Красноярск, просп. им. газ. «Красноярский рабочий», 31

E-mail: telegin@sibsau.ru

Приведены исследования химического состава, морфологии и структуры порошков магнитомягких материалов из кристаллических сплавов 79НМ, 10СЮ и аморфного сплава 71КНСР для их применения в качестве исходных порошков при плазменном напылении ферромагнитных сердечников и магнитных экранов с высокой проницаемостью. Проведен анализ структуры и магнитных характеристик плазмонапыленных покрытий. Установлено, что порошки 79НМ, 10СЮ из-за наличия оксидных пленок по границам зерен имеют магнитную проницаемость около 100 Гс/Э, а в покрытиях она находилась в области 100-300 Гс/Э при напряженности магнитного поля от 0,5 до 3,0 Э. Порошки сплава 71КНСР отличались от кристаллических отсутствием оксидных образований на поверхности частиц и однородностью структуры и имели магнитную проницаемость до 9000 Гс/Э. Максимальная магнитная проницаемость напыленных образцов из порошков сплава 71КНСР при толщине покрытий от 200 до 600 мкм составила 2000-6000 Гс/Э при напряженности магнитного поля 0,10,35 Э.

Ключевые слова: магнитомягкие материалы, аморфный сплав, магнитная проницаемость, плазмотрон, напыление.

RESEARCH OF FERROMAGNETIC POWDER AND THE COATINGS ON THEIR BASIS

V. N. Saunin, S. V. Telegin

Siberian State Aerospace University named after academician M. F. Reshetnev 31, Krasnoyarsky Rabochy Av., Krasnoyarsk, 660014, Russian Federation E-mail: telegin@sibsau.ru