STOCHASTIC SIMULATION APPROACH FOR NONLINEAR SYSTEMS WITH UNCERTAIN INPUT SIGNALS

The present article is devoted to the problem of the nonlinear technical systems, operating under proba-bilistic uncertainty of the input signals and internal parameters. Numerical statistical experiments in the course of simulating such systems involve excessive overhead costs, levelling of which allows to increase signiﬁcantly the precision of estimates of the research system simulated characteristics. The results presented in the paper can be used in solving problems of applied physics, as well as in the analysis and synthesis of control systems. In particular, they can be used for planning and organizing numerical study with limited computing resources, or in case of excessive duration of a statistical experiment.


Introduction
The method of statistical testing known as the Monte Carlo method is often applied to study characteristics of the complex stochastic systems. The system parameters under study are evaluated statistically, the key role being played by their mathematical expected value and deviation from the true value (accuracy) [Ross, 2012]. The required number of the experiments n needed to obtain a sufficiently reliable result and the accuracy of the numerical study ε have the inverse dependence [Danilov, Ermakov, and Halton, 2000;Graham and Ta-lay, 2013]: where α is the value characterizing the confidence interval, and σ 2 is the variance of the sought value estimate. Obviously, the simulation complexity is also affected by the variance of the parameter being estimated [Asmussen and Glynn, 2007;Law and Kelton, 2000]. Numerous methods are applied for the estimation variance reduction [McLeish, 2005]. Depending on the purposes, they allow to achieve either reduction of the statistical simulation complexity or increase of the parameters estimation accuracy. Among the well-known methods are the following: Importance Sampling (IS) [Kroese, Taimre, and Botev, 2011;Rubinstein and Kroese, 2017], Stratified Sampling (SS) [Keramat and Kielbasa, 1998;Marnay and Strauss, 1991], Multiple Control Variates (MCV) [Glasserman, 2003;Glynn and Szechtman, 2002;Nelson, 1990] and Common Part Variates (CPV) [Emeljanov, Likholet, and Sharov, 2009;Podoplekin and Andriyevskiy, 2005], as a analogue of the latter according to the paper [Vasiliev and Sabinin, 1987]. The Importance Sampling method involves replacing the original response function of the model under study with a new one with the same mathematical expectation [Rubinstein and Kroese, 2017]. To minimize the variance the density of the cumulative distribution function of a new random variable is to be proportional to the product of the state variable of the initial random variable by the density of the random parameter variance [Glasserman, 2003;Kroese, Taimre, and Botev, 2011]. On the other hand [Ermakov, 2009], the new model construction in this case is a task comparable in complexity to the original one.
The Stratified Sampling is formed by splitting the intervals of the accepted values of the components of the random parameters vector V = [V 1 , V 2 , . . . , V m ] into the separate layers. The statistical experiment is carried out separately in each layer. When n is fixed this approach in particular ensures the reduction of the cumulative sampling variance [Rubinstein and Kroese, 2017], which in its turn provides the improved accuracy of estimations. With an increasing number of the random process parameters, as well as of layers for each dimension, the solution of the problem of the multidimensional hypercube optimal stratification turns out to be a rather compute-intensive procedure. In order to solve this problem, Latin Hypercube Sampling (LHS) may be applied [Burhenne, Jacob, and Henze, 2011;Owen, 1997;Veetil, Sylvester, and Blaauw, 2008]. The Common Part Variates and Multiple Control Variates methods appear to be the most flexible, in view of their excellent combinability with the sampling methods and the use of the so-called simplified models concept. Some papers introduce the term "metamodel" for a simplified model [Porta Nova and Wilson, 1993]. A simplified model with statistical characteristics similar to those of the base model enables one to define more exactly the results of the statistical experiment and to adjust them. The simplification is achieved by choosing a model which has a less complex structure and requires fewer resources for conducting a single experiment. Both methods assume different criteria for constructing simplified models. In the case of the Common Part Variates method, a simplified model is to have a geometric proximity with the base model, while in the case of the Multiple Control Variates method the proximity based on the strong correlation relationship [Nelson, 1990].
In this paper the application of the Multiple Control Variates method in combination with the unified method of the system simplified model construction proposed by the authors will be considered.
2 Mathematical Description of the Nonlinear Research System As the research system the normalized nonlinear oscillating circuit model considered in the paper [Nekrasov, 2017] will be used: where k = const is the coefficient of nonlinearity, f (t) is the function that determines the system stochastic process. When solving Cauchy problem for the system (2) will be considered the random process f (t) as the following: Here ω = 5π sec −1 , β and µ -are random input signals, uniformly distributed along the interval (0; 4π).
Obviously, in this case the vector of random parameters and input signals is given as the following: V = [β, µ] , m = 2. For the numerical solution of the differential equation (2) in the paper the Runge-Kutta fourthorder method is used under the following initial conditions: 3 Problem Statement When solving Cauchy problem for the nonlinear system of type (2) a statistical experiment inevitably appears to be obligatory in order to determine the mathematical expectation of value x(t) with the predefined precision ε. The application of stochastic simulation method to estimate the desired parameter with the required accuracy implies the iterative adjustment [Glasserman, 2003;Ross, 2012;Rubinstein and Kroese, 2017] of the experiments required number n in the simulation process in accordance with (1). As previously noted the application of the MCV method oriented to minimize the variance of the required parameter implies constructing the simplified model y(t). The application of this method deals with the problem of developing criteria for the selection of the effective approach to form a simplified model. A large number of papers [Asmussen and Glynn, 2007;Glasserman, 2003;Ross, 2012;Nelson, 1990], are known to describe the criteria for forming the structure of y(t) for the particular classes of the systems under study. However the universal approach is known only for various complexity simplified models in the polynomial form [Emeljanov and Likholet, 2008]. The MCV method is characterized by the following dependence of the statistical characteristics of the initial and simplified models: is the covariance [Bashkirtseva, 2016] betweenX andỸ and Var (Ỹ ) is the variance of theỸ .
The universal approach to the construction of y(t) [Emeljanov and Likholet, 2008] implies the search for the polynomial with the minimum of the value −cσỹ/σx of the following form: The set of polynomial coefficients denoted as C = {C 0 , C 1 , . . .} is the set of parameters to be optimized when forming the structure of universal simplified model y(t). In case of nonlinear systems the authors of the polynomial models in order to further improve their productivity [Emeljanov and Likholet, 2008] also developed an adaptive algorithm for stochastic simulation, which implies, in particular, both the advantages of MCV and SS and the formation of simplified models setŶ = {ŷ 0 (t),ŷ 1 (t), . . .} for each SS layer, where ∀ŷ(t) ∈Ŷ is the polynomial simplified model, developed with the respect of the particular layer sampling. Taking into account (5) as well as the possibility to select the optimal position of the layer edges, in case of the adaptive algorithm the maximum number of the parameters to be optimized is estimated as λ = (2m + |C|) Ŷ . Obviously, at larger λ the adaptive construction of the simplified modelsŷ(t) appears to be a very laborious method. This paper presents an alternative approach of constructing universal simplified models which is to be applied as a part of MCV and CPV methods and which assumes only one parameter to be optimized.

Simplified Models Construction Technique
According to [Ross, 2012] MCV implies the following dependence: The substitution of the expressions (6) and (4) into (1) allows to get the estimate of the experiments required number n for the MCV method: The statistical experiment with the iterative adjustment n applying (7) allows to estimate E [X] with the predefined precision ε. The proposed method of the simplified model structure construction implies achieving a controlled correlation coefficientρ xy when the difference between x(t) and y(t) is sufficient for practical application (7) and (4). This requires obeying the following conditions: 1. organization of a statistical experiment in three stages: initialization, post-initialization, iterative; 2. formation of the simplified model structure in the form m-dimentional "heightmap".
The initialization stage is characterized by the random sample N 0 ∈ R m which has the size of n 0 and on the basis of which the "heightmap" is constructed. The sample N 0 is the base one for the simplified models of the considered type. The second stage implies the primary estimation of the desired statistical characteristics of the original model in accordance with (4). The final stage implies an iterative adjustment of the experiments required number (7) and a repeated evaluation of the characteristics. The experiment is terminated with the convergence of the required number of experiments n to the value of theñ and with the achievement of the required accuracy ε. When constructing simplified model y(t) in the form of the "heightmap" the base model x(t) is considered as a function of the point defined at m-dimensional Euclidean space that corresponds to the y( , moreover if ∀p, q ∈ {A, B, C} the following equation is correct: Let A ∈ N 0 , B ∈ N 0 , C / ∈ N 0 and ∀D ∈ N 0 : R DC > R AC , R DC > R BC , then if R AC < R BC and φ = R AC /R BC for simplified model y(t) of the considered form the following equation is also correct: With regard to (9) the task of using y(t) in the form of "heightmap" for ∀C at the post-initialization and iterative stages is reduced to the k-Nearest Neighbor Search (k-NNS) problem in the m-dimensional Euclidean space, which has an effective solution, for example, with usage of the k-d trees [Bentley, 1975].
It should be separately noted that the only parameter, the structure of the y(t) andρ xy in this case depend on is n 0 .

Simulation and Experimental Results
Let's consider the results of the stochastic simulation of the (2) with the presence of the random process (3) in the nonlinear system under study. It also makes sense to experimentally verify the relationship between the parameters n 0 andρ xy .  Figure 1 shows the base model x(β, µ, t) surface visualization where t = const. The model surface is formed as a "heightmap" with the gradient color description of the function resulting value. Further on this illustration will be allowed to assess visually the geometric proximity of the base and simplified models. The original presumption concerning the simplified models built in accordance with (9), is based on an assumption that geometrically similar models are highly likely to have strongly prounounced correlative relationships.
Let's consider the solution of the Cauchy problem for the system (2) applying statistical simulation with the estimation precision ε = 3 × 10 −3 . When comparing the results of the Direct Monte Carlo Simulation and the MCV in combination with the simplified models in the form of the "heightmaps", one can get convinced of achieving the predefined precision, and also can estimateñ for both experiments (valuesñ base andñ M CV respectively). In this case, the following evaluation of the efficiency of the statistical simulation complexity reduction with implementing the proposed approaches can be applied: (10) Figure 2 shows the dynamics of value η with different values of n 0 for adaptive simplified models in the form ofŶ and simplified models y (β, µ, t) in the form of "heightmaps". The fluctuations on chart η for the simplified model Y indicate that in case of multiparametric optimization the global extremum or the local global extremum are hard to be achieved. The revision of theŶ parameters optimization method in some cases may significantly improve the extreme points reachability. Nevertheless, in view of the results presented in this paper this issue has not been studied. The surfaces obtained as the result of solving the above formulated problem with the application of the simplified model in the form of the "heightmap" are shown in Figure 3. For the indicated values of n 0 the following values of the correlation coefficient are achieved:ρ xy = 0.58,ρ xy = 0.83 andρ xy = 0.90 respectively. Obviously, with the increase of the "heightmap" base sample size the visual (geometric) proximity of the simplified and original models significantly increases. This circumstance can be used when projecting the obtained results on the CPV method. The dependence shown in Figure 4 demonstrates the  ability to control characteristicρ xy by determining the optimal size of sample N 0 when using simplified models in the form of "heightmaps" as a part of MCV. For a simplified model ofŶ form, the valueρ xy was not estimated, since in accordance with the mathematical apparatus of the adaptive method [Emeljanov and Likholet, 2008] forŷ(β, µ, t) the degree of correlation with the original model is estimated solely within the confines of a single SS layer and does not spread over the entire range of the accepted values of the random vector V components. The obtained results allows to state that the universal approach for simplified model formation proposed in this paper for the MCV and CPV methods is not inferior in efficiency to the approach proposed in the paper [Emeljanov and Likholet, 2008]. Along with this, it allows to significantly simplify the method of constructing the simplified model structure preserving its key statistical characteristics.

Conclusions
The article demonstrates the universal simplified system construction approach to be used as a part of such methods of stochastic simulation complexity reduction as CPV and MCV. The presented results are comparable in efficiency with the results obtained when the adaptive approach is applied. On the other hand, the key feature of the construction method for a simplified model in the form of a "heightmap" is its single parameter to be optimized. In case of the adaptive approach the number of the parameters to be optimized appears to be considerably larger. This may adversely affect the simplified model construction complexity. As mentioned earlier, the proposed approach can be used to solve the problems of applied physics and the nonlinear stochastic control theory in conditions of limited computing resources. In turn, the fundamental dependence of the simplified model on the single parameter to be optimized makes the practical application of the considered proposals quite attractive.