Автоматическая идентификация линейных дифференциальных уравнений в аналитическом виде

Полный текст

He 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, tt}, i = 1, s , where s is its size, yt e R is dynamic system output measurement at time ti, and ut = u(t{) 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 £,: E(^) = 0, D(Q <да, that affects the output measurements: yi = x(t,) + . (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); ог yi k ). + ak • x(k :) +... + <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 xs"'> + am • x m ) +... + a1 • x = a0 • u(t) , x(0) = x0 , (5) where the vector of equation parameters a = (0, ., 0, am, ..., <5j, a0)T eRn, и <M +1, delivers an extremum to the functional 1 (Хn) = XIу, _ x(ti )l ^ 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 67 Математика, механика, информатика 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 H = (op1, sp1, fitness(opi )^, 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 strategic parameters; N j is the population size; fitness(x): Rk ^ (0,1], fitness(x) =-1- is the 1 + J (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 op<o°fspring = opoffspring + z • N(0,spfspring) , spofvnrn = \spfspnng + z • N(0,1)|, where N (m, a2) 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 т : opfsPring = opfspring + z • N (0, spoffspring ), spoffspring = sr,offsVring eTzN(0,1) Pi ~ 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 ht. 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 68 Вестник СибГАУ. № 1(53). 2014 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 bestt, i = 1, m be the fitness of the best solution found by the i -th algorithm. Similarly, let sizei be the size of the population, NM, N2,i - 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 m Sum _ size = ^ sizet, then the free resources are i=1 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*, bestt e best : bestt > bestt ^ fw (i,best , aw) > f (i,best, aw); aw - are the set of its parameters; round: N ^ R turns its argument to nearest integer. We considered these weight functions: fW (i, best) = ^bet- Free , (7) X bestj j=1 eaw ■besti fW (i, best, aw ) =--Free, aw > 0 . (8) J w 47 ’ w' m w ^ 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 = 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 = I_10, 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. 69 Математика, механика, информатика 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 (Л ’ <л л g <u 1 £-09 (0 ' « 51 08 ад с ' 2 .2 о 7 a) -и ' £ | 0,6 <U 4- 0,5 123456789 10 The order of the ODE model 123456789 10 The order of the ODE model c 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) The order of ODE The order of ODE б Fig. 2. Varying the number of local optimization steps (a), and size of population (b), for different orders of model а а 70 Вестник СибГАУ. № 1(53). 2014 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 coevolution-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.

Об авторах

Иван Сергеевич Рыжиков

Сибирский государственный аэрокосмический университет имени академика М. Ф. Решетнева

Email: ryzhikov-88@yandex.ru
аспирант Института информатики и телекоммуникаций

Список литературы

Дополнительные файлы