One dynamical input reconstruction problem: tuning of solving algorithm via numerical experiments

1.   Introduction

The problems of reconstructing unknown input parameters of controlled systems based on inaccurate/incomplete information on the phase state arise in many scientific studies and applications (in flight mechanics and guidance theory, in financial and actuarial mathematics, for designing engineering and production processes, in ecology and medicine, etc.) and attract a great attention in recent time. The first publications concerning this field go back to the 60's of the previous century ^[1,2,3], when some criteria of the unique solvability of inverse problems and of the continuous "input/output" dependence for systems described by ordinary differential equations were obtained under the conditions of sufficiently smooth inputs. The key definition of an invertible dynamical system was introduced in ^[2].

Input reconstruction problems can be treated as the special case of identification problems. They fall into the range of inverse problems of dynamics of controlled systems and, as a rule, are ill-posed (due to the nonuniqueness of the desired characteristic, as a consequence of inaccuracies, and/or due to the discontinuity of the inverse operator) and require the application of regularizing procedures.

The huge amount of works is devoted to different statements (for example, to operational statements), where solving approaches are based on a posteriori character of desired regularizing algorithms, which can use the whole history of output measurements. We list only some of monographs containing detailed reviews and bibliography. In the classical monograph ^[4], the emphasis is on different classes of recurrent methods for the identification of nonstationary objects, including the theoretical analysis and experimental verification. The monograph ^[5] contains the introduction to identification theory. Algorithms for estimating time-varying inputs acted upon linear systems are considered in the monograph ^[6]. The monographs ^[7,8] are devoted to the foundations of the theory of inverse and ill-posed problems.

One specific solution approach named the method of dynamic inversion is extremely important in cases when there is the necessity to reconstruct some characteristics dynamically and preferably in real time. It was proposed and developed by Kryazhimskii, Osipov, and their colleagues, see ^[9,10,11,12]. The approach is based on the combination of principles of the theory of positional control ^[13,14] and ideas of the theory of ill-posed problems ^[7,8]. A reconstruction problem is reduced to a feedback control problem for an auxiliary dynamical system called a model. The goal is to find an appropriate control law for this model in such a way that the adaptation of the model controls to the results of current observations provides an approximation (in a desired sense) to the unknown inputs. We want to have solving algorithms both dynamical and stable. As the dynamic property, we understand the property that the algorithm uses current input and output values in real time for decision making during the process. In its turn, the stability property means that the algorithm produces an approximation that is as precise as one likes if the observations are sufficiently accurate.

First, the method of dynamic inversion was realized for systems described by ordinary differential equations (ODEs) ^[9,10,11], then, for functional differential equations ^[15,16], equations and variational inequalities with distributed parameters ^[12,17]. The case of partially observed systems was also intensively investigated ^[11,18].

As to the application of the theory of dynamic inversion to stochastic objects, the problem of positional simulation of an unknown stochastic control in a system described by an ODE was considered for the first time in ^[19]. The present paper continues the investigations (in the framework of the theory above) of the problems of reconstructing unknown deterministic characteristics of a linear/quasi-linear stochastic differential equation (SDE), see ^[20,21,22]. For a linear SDE, the problem of dynamical reconstruction of a disturbance entering into the Ito integral and characterizing the amplitude of a random noise, on the basis of the measurements of realizations of the whole phase vector, was considered in ^[20]. Then, the case was supplemented by the possibility of measuring a part of coordinates of the stochastic process ^[21]. The inverse problem for a quasi-linear system with diffusion depending on the phase state in the statement assuming the simultaneous reconstruction of disturbances both in the deterministic and stochastic terms was discussed in ^[22]. Namely the latter problem is under further investigation in the present paper. As is shown in ^[22], the problem for a SDE can be reduced to an inverse problem for the system of ODEs describing the mathematical expectation and covariance matrix of the process. The applicability of the algorithms for reconstructing unknown parameters developed earlier for systems of ODEs ^[9,10,11] is substantiated, an appropriate finite-step software-oriented modification is proposed, and an estimate for its convergence rate with respect to the number of measurable realizations is obtained, see ^[22].

To apply practically algorithms of such kind is technically rather difficult due to the necessity of adjusting the algorithm to a specific dynamical system. The matter is that up to now there is no a universal procedure for fitting parameters of algorithms of dynamical reconstruction even in the case of ODEs, see the discussions on numerical examples in ^[10,11,16]. In fact, this important applied aspect of the theory of dynamic inversion remains vague. In the paper, we make a first attempt to fill this gap. The idea is in the automatic tuning of algorithm's parameters for the given system using some model disturbances, which are typical for the problem, i.e., which possess the features (a priori known) inherent in real disturbances (for example, subject to the same geometrical/integral constraints, the same bounds on variation). After performing this empirical procedure, we hope that the set of parameters obtained will be suitable for the reconstruction of other disturbances acting on the system. For the problem considered, such a process assumes numerous program runs and the simulation of a large amount of SDE's trajectories. Thus, the novelty of the present paper consists namely in considering the applied aspects of solving reconstruction problems. To optimize the process of the adaptation of the algorithm to a system, the parallelization of time-taking calculations is involved.

The paper is organized as follows. In Section 1, we give the statement of the dynamical input reconstruction problem for a quasi-linear system of SDEs. Section 2 is devoted to the description of the procedure of reducing the problem for the system of SDEs to the problem for a system of ODEs. A constructive algorithm for solving the latter problem and estimates of its convergence rate are presented in Section 3. The first three sections contain the basic facts from ^[22], which are necessary to formulate the main results of the present paper. Different aspects of the automatic tuning of algorithm's parameters in order to get best approximations in a specific dynamical system and the necessity of parallelizing calculations are discussed in Section 4. In Section 5, a model example illustrating and testing the method proposed is given. In Conclusion, we resume the results obtained and outline some directions for perspective investigations.

2.   Reconstruction problem statement

Simple linearized models can be useful when studying the price dynamics on a commodity market under the influence of random factors, the change of the size of a biological population in a stochastic medium or the processes of chaotic motion of one-type particles. In this context, we consider a quasi-linear system of SDEs with diffusion depending on the phase state:

Here, $t\in T = [0, \vartheta]$, $x = (x_1, x_2, \ldots, x_n) \in \mathbb{R}^{n}$ is a column vector, $\xi \in \mathbb{R}$; $x_{0}$ is a known deterministic or random (normally distributed) vector of initial conditions; $\omega \in \Omega$, $(\Omega, F, P)$ is a probability space; $\xi(t, \omega)$ is a standard scalar Wiener process (i.e., a process starting from zero with zero mathematical expectation and dispersion equal to $t$); $A(t) = \{a_{ij}(t)\}$, $B(t) = \{b_{ij}(t)\}$, and $f(t) = \{f_{i}(t)\}$ are continuous matrix functions of dimension $n\times n$, $n\times r$, and $n\times 1$, respectively. Two external disturbances act on the system: a vector $u_1(t) = (u_{11}(t), u_{12}(t), \ldots, u_{1r}(t)) \in \mathbb{R}^{r}$ and a diagonal matrix $U_2(t) = \{u_{21}(t), u_{22}(t), \ldots, u_{2n}(t) \} \in \mathbb{R}^{n\times n}$. The action $u_1$ enters into the deterministic term and influences the mathematical expectation of the desired process. Since $U_2 x d\xi = (u_{21} x_1 d\xi, u_{22} x_2 d\xi, \ldots, u_{2n} x_n d\xi)$, we can assume that the vector $u_2 = (u_{21}, u_{22}, \ldots, u_{2n})$ regulates the amplitude of random noises. Let the vectors $u_1$ and $u_2$ take values from given convex compact sets $S_{u1}$ and $S_{u2}$, and let their elements belong to the spaces $L_{2}(T; \mathbb{R}^{r})$ and $L_{2}(T; \mathbb{R}^{n})$, respectively.

A solution of system (2.1) is defined as a stochastic process satisfying (for any $t$ with probability 1) the corresponding integral identity containing the stochastic Ito integral on the right-hand side. As is known, under the above assumptions, there exists a unique solution, which is a normal Markov process with continuous realizations, see ^[23], Theorem 5.2.1.

The problem under discussion consists in the following. At discrete, frequent enough, times $\tau _{i} \in T$, $\tau _{i} = i\delta$, $\delta = \vartheta /l$, $i \in [0:l],$ the information on some number $N$ of realizations of the stochastic process $x(\tau _i)$ is received. We assume that $l = l(N)$ and there exist an estimate $m_{i}^{N}$ of the mathematical expectation $m(t) = Mx(t)$ and an estimate $D_{i}^{N}$ of the covariance matrix $D(t) = M(x(t)-m(t))(x(t)-m(t))'$ (the prime means transposition) such that

and $h(N), \, g(N) \rightarrow 0$ as $N\rightarrow \infty$. By $\| \cdot \|$, we denote the Euclidean norm in the corresponding space. The standard statistical procedures ^[24], Ch. 22, of constructing the estimates $m_{i}^{N}$ and $D_{i}^{N}$ admit modifications providing the validity of (2.2).

It is required to design an algorithm for the dynamical reconstruction of the unknown disturbances $u_1(\cdot)$ and $u_2(\cdot)$ generating the stochastic process $x(\cdot)$ from the discrete information on its realizations. The probability of an arbitrarily small deviation of approximations from the desired inputs in the metric of the spaces $L_{2}(T; \mathbb{R}^{r})$ and $L_{2}(T; \mathbb{R}^{n})$ should be close to 1 for sufficiently large $N$ and the time discretization step $\delta = \delta(N) = \vartheta /l(N)$ concordant with $N$ in an appropriate way.

The specific properties of quasi-linear system (2.1) admit the reduction (by analogy with the method of moments, see ^[25]) of the problem formulated for a SDE to a problem for the system of ODEs describing the mathematical expectation and covariance matrix of the desired process. This allows us to organize a procedure of the simultaneous reconstruction of the disturbances both in the deterministic and stochastic terms of the right-hand side in the case when simultaneous measurements of a sufficiently large amount of trajectories (for example, of the motion of one-type particles) are possible.

3.   Reduction of the problem

Following ^[20,22], we reduce the reconstruction problem for a quasi-linear SDE to a problem for a system of ODEs. We need

Lemma. The standard estimates $m_{i}^{N}$ of the mathematical expectation $m(\tau _i)$ and $D_{i}^{N}$ of the covariance matrix $D(\tau _i)$ constructed from $N$ $(N > 1)$ realizations $x^{1}(\tau_{i}), x^{2}(\tau_{i}), \ldots, x^{N}(\tau_{i})$ of the random variables $x(\tau_{i})$ by the rules ^[24], Ch. 22:

provide the validness of relation (2.2).

Note that estimate (3.2) can be rewritten to optimize calculations (in particular, to parallelize them) in the form:

The proof of the lemma is presented in ^[22], where the following parameters of estimate (2.2) are obtained:

Here, $0 < \epsilon < 1/2$, $0 < \alpha < 1$ are parameters allowing us to regulate the orders of smallness in the approximating and probabilistic parts of estimate (2.2). Obviously, $\delta(N) \rightarrow 0$, as well as $h(N), \, g(N) \rightarrow 0$, as required in relation (2.2), for any admissible values $\epsilon$ and $\alpha$. By $C_i$ we denote auxiliary constants, which depend on the structure of the system and on the specified constraints on disturbances (and independent of the functions to be reconstructed) and can be written explicitly.

We briefly describe how to derive the system of ODEs for the mathematical expectation and covariance matrix of the desired process (the detailed procedure is given in ^[20,22]). Introduce the notation: $m_{0} = Mx_{0}$, $D_{0} = M(x_{0}-m_{0})(x_{0}-m_{0})'$. Since the original system is quasi-linear and the mathematical expectation of an Ito integral is zero, the value $m(t)$ does not depend on $u_2(t)$; its dynamics is described by the equation

The covariance matrix $D(t)$ does not depend on $u_1(t)$ explicitly; its dynamics is derived according to the scheme of the method of moments ^[25]:

This matrix equation is rewritten in the form of a vector equation, which is more traditional for the problems under consideration. By virtue of the symmetry of the matrix $D(t)$, its dimension is defined as ${n_{d} = (n^2+n)/2}$. We introduce the vector $d(t) = \{d_{s} (t)\}, \, s\in [1:n_{d}]$, consisting of successively written and enumerated elements of the matrix $D(t)$ taken line by line starting with the element located at the main diagonal. The coordinates of this vector are found from the elements of the matrix $D(t) = \{d_{ij} (t) \}, \, i, j\in [1:n]$:

Then, the matrix equation is written in the form

Here, the matrix ${\bar{A}(t):\, T\rightarrow R^{n_{d}\times n_{d}}}$ and the diagonal matrix ${\bar{B}(d(t), m(t)):\, T\rightarrow R^{n_{d}\times n_{d}}}$ can be written explicitly ^[22], the vector $u_3(t) = \{u_{3s} (t)\}$, $s\in [1:n_{d}]$ is found from $u_2(t)$ by the formulas $u_{3s} = u_{2i} u_{2j}$, $i\leq j$, $s = (n-i/2)(i-1)+j$, it takes values from some convex compact set $S_{u3} \in \mathbb{R}^{n_d}$ and has elements belonging to the space $L_{2}(T; \mathbb{R}^{n_d})$. The initial state $d_{0}$ is obtained from $D_{0}$, whereas the measurements $d_{i}^{N}$, from $D_{i}^{N}$.

Thus, we have the system of equations of dimension $n+n_d$:

Now, we can reformulate the original problem. Taking into account (2.2), we assume that during the process, at the discrete times $\tau _{i} \in T$ the inaccurate information $m_{i}^{N}, \, d_{i}^{N}$ on the phase state of system (3.5), (3.6) is received:

where $h(N), \, g(N) \rightarrow 0$ as $N\rightarrow \infty$. It is required to design an algorithm of dynamical reconstruction of the unknown disturbances $u_1(\cdot)$ and $u_3(\cdot)$ from information (3.7). The probability of an arbitrarily small deviation of approximations from the desired inputs in the metric of the spaces $L_{2}(T; \mathbb{R}^{r})$ and $L_{2}(T; \mathbb{R}^{n_d})$, respectively, should be close to 1 for sufficiently large $N$ and the time step of measurements $\delta = \delta(N) = \vartheta /l(N)$ concordant with $N$ in an appropriate way. Note that, under additional assumptions, the function $u_2(\cdot)$ is also reconstructed by means of the coordinates of the vector $u_3(\cdot)$, which are equal to $u^2_{2i}$, $i\in [1:n]$.

System (3.5), (3.6) is nonlinear in phase variables but is linear in control. The problem formulated for this system corresponds to the problem considered, for example, in ^[9,10,11]. In the present paper, for the modified version of the finite-step solution algorithm proposed earlier for ODEs, we use the fact that it admits a constructive concordance of its parameters with the number of measurable realizations of the original stochastic process.

4.   Reconstruction algorithm

Let us present the solving algorithm for system (3.5), (3.6). At the initial time $\tau _{0} = 0$, we fix a value $N$; then, determine the values ${l^{N} = l(N)}$, ${h^{N} = h(N)}$, and $g^{N} = g(N)$ (see (3.4)) and construct the uniform partition of the interval $T$ with the step $\delta^{N} = \vartheta /l^N$: $\tau_{i} \in T, \, \tau_{i} = i\delta^{N}, \, i\in [0:l^N]$. We introduce a discrete model, which dynamics and initial state are defined by the relations

Here, $i\in [0:l^N-1]$, $w_m$, $w_d$ are model variables approximating the phase state of system (3.5), (3.6); $v_{1i}^N$, $v_{3i}^N$ are control actions calculated at the time $\tau_{i}$ by rules that are specified below.

The work of the algorithm is decomposed into $l^N$ identical steps. At the $i$th step performed on $(\tau_{i}, \, \tau_{i+1}]$, the input data for calculations are the estimates $m_{i}^{N}$, $d_{i}^{N}$, and the model state $w_m(\tau_{i})$, $w_d(\tau_{i})$ obtained by this moment. The model controls are found as follows: $v_{1i}^{N}$ and $v_{3i}^{N}$ are unique solutions of the corresponding extremal problems

where $\langle \cdot, \, \cdot \rangle$ is the scalar product in the corresponding Euclidean space, $\alpha_{m}^{N}$, $\alpha_{d}^{N}$ are regularization parameters. After the calculation of controls (4.3) and (4.4), the model state $w_m(\tau_{i+1})$, $w_d(\tau_{i+1})$ is recomputed by formulas (4.1) and (4.2). The process stops at the terminal time $\vartheta$.

Let $U_{*} = U_{*}(m(\cdot), d(\cdot))$ be the set of all disturbances $u(\cdot) = (u_1(\cdot), u_3(\cdot)) \in L_2(T; \mathbb{R}^{r+n_d})$ generating the pair $(m(\cdot), d(\cdot))$. An element of the set $U_{*}$ that has the minimal $L_2(T; \mathbb{R}^{r+n_d})$-norm is denoted by $u_{*}(\cdot)$. Its existence and uniqueness follow from the convexity of $U_{*}$ and the strict convexity of the norm in $L_2(T; \mathbb{R}^{r+n_d})$, see ^[9,11]. The function $u^{N}(t) = (v_{1}^{N}(t), v_{3}^{N}(t))$, $v_{1}^{N}(t) = v_{1i}^{N}$, $v_{3}^{N}(t) = v_{3i}^{N}$, $t\in [\tau _{i}, \tau _{i+1})$, $i\in [0:l^N-1]$ is found from relations (4.3) and (4.4). Below is the main result of ^[22] concerning the convergence of the algorithm.

Theorem. Let the following conditions of concordance of the parameters hold for $N\rightarrow \infty$:

Then, the sequence $\{u^{N} (\cdot)\}$ is compact in $L_2(T; \mathbb{R}^{r+n_d})$ and we have the convergence as $N\rightarrow \infty$:

Under the assumption that the variation of real disturbances is bounded, the following estimate for the accuracy of the algorithm with respect to the number of measurable realizations of the process is valid:

where $D_1$ and $D_2$ are some positive constants.

As is shown in ^[22], to obtain estimate (4.7), it is necessary to set

where $K_1$, $K_2$, and $K_3$ are some positive constants.

Note that the constants in (4.7) and (4.8) depend on the parameters of the given system, on structural a priori constraints, in particular, on the variation of real disturbances, but are independent of the functions under reconstruction. The choice of the power exponents of the value $1/N$ in the approximating and probabilistic parts of estimates like (4.7) admits certain arbitrariness in the range of change of $\epsilon$ and $\alpha$ from relations (3.4). In the theorem, we orientate to the coincidence of the orders of smallness for the parts, it is provided by relations (4.8). Evidently, other variants are possible. In the case when the vector $u_2 = (u_{21}, u_{22}, \ldots, u_{2n})$ is unique and its coordinates are nonnegative (such constraints seem natural to characterize a noise with zero mean value; below, we assume that these conditions are fulfilled), algorithm (4.1)–(4.5) reconstructs the function $u_2(\cdot)$ itself in the metric of the space $L_2(T; \mathbb{R}^{n})$ with an accuracy estimate of form (4.7).

5.   Tuning of algorithm's parameters

The number $N$ of measurable trajectories of SDE (2.1) determines both the accuracy $h^{N}$ of input estimates (2.2), (3.7) and the parameters of the algorithm $\delta^{N}$, $\alpha_{m}^{N}$, and $\alpha_{d}^{N}$ providing convergence (4.6), (4.7). Relations necessary for calculations are not specified explicitly but are given as the order of convergence (as $1/N \rightarrow 0$) with accuracy to constants (4.8), for which upper bounds can be written while their values need more precise definition. Note that there is an empirically proven sensitivity of the algorithm to the constraints for the functions forming system (2.1), to the sets $S_{u1}$ and $S_{u2}$ of admissible values of the disturbances, to their variation, etc. It is necessary to find constants in relations (4.8) providing the successful reconstruction of any possible disturbances (satisfying a priori constraints) in a specific system of form (2.1). Up to now, there is no a universal procedure for fitting parameters like $\delta^{N}$, $\alpha_{m}^{N}$, and $\alpha_{d}^{N}$ of reconstruction algorithms ^[10,11,16]. To get over this difficulty, we suggest to use some model disturbances, which are typical for the problem, i.e., which possess the features (a priori known) inherent in real disturbances (subject to the same geometrical/integral constraints, the same bounds on variation, and so on). This is what the fine tuning of the algorithm for. At the current stage, we assume that "optimal" relations between parameters can be determined empirically through a solution of the following extremal problem:

Here, $L_2(T) = L_2(T; \mathbb{R}^{r+n_d})$, $u^{N}_{K_{1}, K_{2}, K_{3}}(\cdot)$ is an algorithm's output depending on the constants from (4.8); $I^{N}$ is the set of chosen values of $N$, $\beta^{N}$ are weighting coefficients, $\sum_{N \in I^{N}} \beta^{N} = 1$; $S_{K}$ is some set of admissible values of the vector $(K_{1}, \, K_{2}, \, K_{3})$, as a rule, it is a parallelepiped in the positive octant. We fix all the constraints imposed on the admissible inputs, fix some test function $u(\cdot)$ satisfying these constraints, introduce a uniform grid with respect to $K_{1}$, $K_{2}$, $K_{3}$ and then solve problem (5.1) by exhaustive search. Such a procedure requires numerous runs of the calculative module with modeling a large amount of independent trajectories of SDE (2.1) to form estimates (2.2) for different $N$. We hope for the applicability of this empirical optimization procedure providing a solution of (5.1) to the reconstruction of different disturbances satisfying given constraints in the specific system, since the test function has all main features, which are inherent in real inputs and influencing the constants in relations (4.8). Numerous experiments (one of them is described in the next section) demonstrate results in favor of our hypothesis. To simulate the dynamics of equation (2.1), we use the Euler method with the substitution of a sequence of random impulses for the Wiener process (the detailed description of the approximation scheme can be found, for example, in ^[26]). The mean square accuracy order of this method for a quasi-linear SDE is equal to $O(\delta_s)$, where $\delta_s$ is a simulation step. This step is, as a rule, rather small comparing with the step $\delta^{N}$ of receiving the information on the realizations and, correspondingly, of reconstructing the unknown function $u(t)$. Thus, two time grids are involved into the scheme of the tuning of algorithm's parameters for a specific system. The simulations with the smaller step take place only at the stage of tuning; when solving a real problem, the input information should come from "outer space". Note that the probabilistic character of convergence estimate (4.7) assumes the algorithm's output averaging over a series of runs for a fixed set of parameters.

The numerical experiments showed that the tuning procedure is rather time-taking (since it requires a large amount of program runs on the grid with respect to several parameters) when sequential schemes are used but it admits a natural (and effective) parallelization. The latter is based on the standard scheme "master-worker" ^[27,28] with a unique loading MPI-module, almost uniform computational load for all processor cores, and the absence of the informational exchange between the workers. Actually, the simulation of a large amount of independent SDE trajectories with a small time step is the only time-taking part of the procedure; therefore, it is shared uniformly among all the cores. At every moment, which is an entry point for coming input information on estimates (2.2), (3.1), (3.3), all the modules compute the same size portions of trajectories $x^{r}(\tau_{i})$ and average both their values itself and their squares, according to formulas (3.1), (3.3). Then, the workers send their contributions into the estimates to the master, who adds own portion and, averaging the sums obtained over the number of cores (the legitimacy of this operation is provided by the same size of portions), finally forms (3.1) and (3.3). All the calculations in the framework of reconstruction algorithm (4.1)–(4.4) is performed by the master. The flowchart of the main fragment of the parallel procedure is presented in Figure 1.

6.   Illustrative example

Consider a quasi-linear SDE describing a stochastic process, which can be treated as a "spoiled" mean-reverting Ornstein–Uhlenbeck process ^[23]:

Here, $u(t) = (u_1(t), \, u_2(t))$, $u_1(t) = (u_{11}(t), \, u_{12}(t))$ and $u_2(t) = (u_{21}(t), \, u_{22}(t))$ are unknown disturbances to be reconstructed, $u_{1i}(t) \in [-2, 2]$, $u_{2i}(t) \in [0, 2]$, $i = 1, 2$; $\xi(t)$ is a standard scalar Wiener process. Equations like (6.1) are used, in particular, in some simplest models describing the dynamics of relatively stable biological populations ^[23]. In this case, the value $x(t) = (x_1(t), \, x_2(t))$ represents the current size (in arbitrary units) of two interacting biological species; the structure of the deterministic part of the equation defines some sizes (for example, averages for some previous interval), which are the aim of "subconscious" return of the populations. Such a return can be prevented by the neighbors as well as by the action of external factors via the disturbance $u_1(t)$ influencing the mathematical expectation of the process and $u_2(t)$ characterizing the amplitude of random noises. Realizations of $x(t)$ corresponding to different colonies are measurable at discrete times. Let us write for (6.1) the system of ODEs of form (3.5), (3.6) describing the mathematical expectation and symmetric covariance matrix:

Here, $u_{31} = u_{21}^2$, $u_{32} = u_{21}u_{22}$, $u_{33} = u_{22}^2$.

The main goal of numerical experiments was the tuning of the algorithm for specific dynamical system (6.1)–(6.2), i.e., the derivation of "universal" relations between parameters $\delta^{N}$, $\alpha_{m}^{N}$, $\alpha_{d}^{N}$, and $N$ via $K_{1}$, $K_{2}$, $K_{3}$ (see (4.8)) to reconstruct any disturbance satisfying the given constraints. In the first series, to solve problem (5.1), we chose the functions

and parameters $K_i \in [0.1, 5]$, $i = 1, 2, 3$, with $10-20$ grid nodes, $I^{N} = \{10^{3}, 10^{6}\}$, $\beta^{N} = \{1/2, 1/2\}$. Here, we took into account a priori information and the constraints imposed on real disturbances. Optimal (in the sense of (5.1)) values of the constants from (4.8) were obtained: $K^{*}_{1} = 0.6$, $K^{*}_{2} = 0.35$, $K^{*}_{3} = 3.5$. These values determine the parameters $\delta^{N}$, $\alpha_{m}^{N}$, and $\alpha_{d}^{N}$ of the algorithm. The results of reconstructing the functions $u_1(t)$ and $u_2(t)$ for found $K_i$ and different $N$ are presented in Figures 2 and 3. These results are principally in agreement with the main assertion of the paper, thus confirming the convergence like (4.6): the more is $N$, the less is the reconstruction accuracy.

In addition, it should be noted that the output accuracy obtained in the variant from Figure 3 seems satisfactory ($\| u^{N}(\cdot) - u(\cdot) \|_{L_2(T)} = 0.118$, whereas the measuring accuracy $h ^{N}$ of estimates (2.2), (3.7) does not exceed $0.001$ ($h ^{N}$ is close to the order $1/N^{1/2}$, see (3.4), $N = 10^{6}$)). The component $u_1(\cdot)$ is reconstructed considerably better than $u_2(\cdot)$; the matter is in the differences in accuracies of the statistical estimates $m_{i}^{N}$ of the mathematical expectation and $D_{i}^{N}$ of the covariance matrix.

Then, the parameter relations obtained were tested on other disturbances meeting the given geometrical constraints. Let us give the example (see Figure 4) for $u_{2}(t) = (\sqrt{2}(1-2|t-0.5|), 1+ 0.5\sin{30t})$ where all the values of algorithm's parameters were the same as in the variant from Figure 3. The reconstruction accuracy is comparable with the results of the latter variant. Note that the second coordinate of the disturbance $u_{2}$ is reconstructed worse than the first one; it can be explained by the larger amplitude of the noises in the second equation of system (6.1). Additional experiments demonstrated similar results.

During the experiments, we observed quite different interdependencies between the accuracy and parameters of the algorithm: from rather trivial, a monotone function, see Figure 5 (left), up to a function having a minimum inside the interval, see Figure 5 (right). Both variants are not stable; this fact complicates the optimization of the exhaustive search applied when solving problem (5.1).

To test the parallelization quality characteristics, we chose the variant shown in Figure 3 with the simulation of 128000 trajectories of equation (6.1). The calculations were performed at the Krasovskii Institute of Mathematics and Mechanics of UB RAS by means of a hybrid machine of cluster type named "Uran" with the peak performance about 215 TFlops. The results of the test are presented in Table 1. We measured (in seconds) the computation time for one fixed set of parameters $(K_{1}, \, K_{2}, \, K_{3})$ in problem (5.1), whereas the number of such variants in the example described above was up to 8000. Here, $T_p$ is the performance time on $p$ cores (then, $T_1$ is the performance time for the sequential algorithm), $S_{p} = T_{1}/T_{p}$ and $E_{p} = S_{p}/p$ are speedup and efficiency, respectively ^[27,28]. Recall that the speedup is equal to the number of engaged cores for an "ideal" parallel algorithm; in this case, we have the unit efficiency.

It follows from Table 1 that the parallelization efficiency is rather high and does not fall below a reasonable level if the number of engaged cores increases. For example, using 128 cores, we get the acceptable computation time for solving optimization problem (5.1). It is important since we need numerous program runs both for the tuning of the algorithm and for its testing in the probabilistic sense.

7.   Conclusion

In the present paper, the inverse problem for a quasi-linear diffusion SDE consisting in the dynamical reconstruction of two unknown nonrandom disturbances in the deterministic and stochastic terms of the equation is investigated. As input information, measurements of several realizations of the stochastic process at discrete times are used. The problem is reduced to the inverse problem for the system of two ODEs describing the mathematical expectation and covariance matrix of the original process. To solve the latter problem, the feasible reconstruction algorithm is designed; the estimate for its convergence rate with respect to the number of measured realizations is obtained. The main result of the paper is in an attempt to develop an empirical procedure for the automatic adaptation of the algorithm to a specific dynamical system and constraints imposed on its structure and on possible unknown disturbances. This implies numerous program runs on a parameter grid and the simulation of a large amount of SDE's trajectories to obtain necessary statistical estimates for the moments. Toward this aim, the parallelization of time-taking calculations is applied. As perspective directions of investigations, we mark out the theoretical substantiation of some procedure (described above or similar) for the automatic tuning of a reconstruction algorithm to a system via finding an extremum of an appropriate quality criterion over a set of admissible parameters and the approbation of the method on other dynamical input reconstruction algorithms.

Acknowledgments

The authors would like to thank the generous support from the Krasovskii Institute of Mathematics and Mechanics, Yekaterinburg.

Conflict of interest

The authors declare no conflicts of interest in this paper.