Dynamic System Identification via Randomized Stochastic Optimization Under Unknown-but-Bounded Noise

1. Introduction

An exact solution to any problem is possible with an accurate formulation of the problem, but connections and relationships in the real world are so complex and diverse that it is almost impossible to describe many phenomena strictly in mathematical language. A typical approach in theory is to select a mathematical model close to real processes and include various noises related, on the one hand, to the roughness of the mathematical model and, on the other hand, characterizing uncontrolled external disturbances affecting the considered object or system. For all mathematical models, the result of the experiment is a mathematical object: a number, a set of numbers, a curve, etc. From a mathematical point of view, a significant range of applied problems aims to restore characteristics from experimental data (parameters) of the object. At the same time, real systems are rarely described thoroughly by limited mathematical models. When choosing a model to solve a real problem, it is common to consider so-called systematic error (model error), which can be quantified by the distance from the real operator to the selected model. Another type of error that an experimenter may encounter is associated with measurement errors. Such errors are called statistical errors (random errors). The process of selecting characteristics (parameters) of a model from a given class of models to describe the results in a best way is one of the general definitions of estimation. In practice the estimation process can often be related to some quantitative characteristic of the quality of estimation and, it is natural while choosing the estimates, to try to minimize the negative impact of errors, both statistical and, if possible, systematic [2]. Examples include:

• formation of multiscale vortex structures in turbulent fluid flows and plastic flow of solids under pulsed load;
• clustering in a stream of concentrated dispersive mixtures;
• propagation of the shock wave front inside the substance;
• transition layers near interphase boundaries;
• processes of protein formation in cells;
• hierarchy of structures in living systems;
• processes of fission of heavy elements nuclei;
• thermonuclear fusion;
• as well as the behavior of groups of people.

2. Discretization in System State Space

Among new directions in distributed systems research connections between distributed systems theory on the one hand, and canonical problems in turbulence and statistical mechanics on the other could be suggested. In one class of problems, spatio-temporal dynamical analysis clarifies old and complex questions in the theory of shear flow turbulence. In another class of problems, structured, distributed control design exhibits dimensionality-dependence and phase transition phenomena similar to those in statistical mechanics.

Figure 1. Mesoscale system structure can occur as a result of external influence at the system with microscale structure. It exhibits the properties of both micro- and macroscale structure. On the one hand, the clusters at meso-scale can be considered as a set of different structures, as single particles at microscale but of larger size. On the other hand each cluster could be viewed as a rigid structure, as a macroscale object but of smaller size.

Figure 1. Mesoscale system structure can occur as a result of external influence at the system with microscale structure. It exhibits the properties of both micro- and macroscale structure. On the one hand, the clusters at meso-scale can be considered as a set of different structures, as single particles at microscale but of larger size. On the other hand each cluster could be viewed as a rigid structure, as a macroscale object but of smaller size.

Assume that (unidirectional) time t is introduced, and consider non-isolated systems consisting of elements. Evolution of each system over time is determined by the current states of both the system itself and other elements of the system. The evolution is also affected by external disturbing influences W (the absence of influence can be interpreted as zero impact). External influences W can be formally included in the general set of system states. The inclusion of W into of the system state can significantly complicate the descriptive model. External influences W naturally fall into two groups:

$$W=\left(\begin{array}{c}u\\ w\end{array}\right),$$

(1)

controlled ones u (or simply, control) and uncontrolled ones w.

It is usually assumed that at time t the system state $X\left(t\right)$ is finite-dimensional, and to describe the dynamics of the system a system of differential equations is used

$${\dot{x}}_i=g_i(X,W),X=\left\{x_i\right\},i\in \mathbb{M}=\{1,2,\dots ,m\}$$

(2)

with some functions $g_i(·)$ and external disturbances W. In models of complex systems of this type, consisting of a huge number of components, it is customary to assume a large dimension n of state space. But, on the one hand, the choice of the threshold for n significantly limits the “upper bound” of complexity, and, on the other hand, does not allow to take possible “flexibility” of the system into account during its changing process. We will assume that the system consists of a continuum of elements $X=\left\{x_{\gamma }\right\}$, parameterized by $\gamma \in [0,1]$, and the evolution of each of the elements is described by equation

$${\dot{x}}_{\gamma }=g_{\gamma }(X,W),X=\left\{x_{\gamma }\right\},\gamma \in \mathbb{M}=[0,1].$$

(3)

Such setups can occur, for instance, in stochastic games with many players [3] and mean-field games with almost infinitely many players [4]. Additionally, we assume that the external influence W for all its arbitrariness at each moment of time k has a structure $s_k$ of finite order. Finite structure of external influence after some transition process causes discretization of spatial elements (clustering) in the considered complex system

$${\mathcal{X}}_{s_k}=\{X_1,X_2,\dots ,X_{m_{s_k}}\}:X={\cup }_{i=1,2,\dots ,m_{s_k}}{X}_i,X_i\subset X.$$

(4)

Discretization occurs due to the self-organization of groups of elements and their synchronization. For clusters $i=1,\dots m_{s_k}$, a set of $m_{s_k}\in N$ variables ${\overline{x}}_i$ averaged over cluster $X_i$ is naturally introduced. Set of ${\overline{x}}_i$, $i=1,\dots ,m_{s_k}$ could be generalized as a set of some integrals over the clusters. Such approach is usually applied to simplify physical models (dimension reduction). The general integral characteristics of clusters of elements with similar properties are introduced and dynamic models in reduced state spaces are considered. Experiments show that such simplifications are justified and often give good results.

Typically, the process of clustering (self-organization) in a system is not “one-shot”, but is constantly reproduced due to changes in external influences and critical changes in internal states. But at the same time we will assume that a change in the structure of external influences does not occur permanently, but in some time instants $T_0,T_1,T_2,\dots .$ This leads to the necessity to consider the dynamic processes under condition of state space structure change over time.

So, when the structure of external influence changes, the discretization of spatial elements may change. Let us assume that this transient process takes a duration of time no longer than some $\delta \ge 0$ (see Figure 2). We will assume that $\delta$ is many times smaller than the intervals between successive changes in the structure of external influences:

$$\delta <<\zeta =\underset{k}{min}|T_{k+1}-T_k|$$

(5)

In addition to discretizing spatial variables we obtain time sampling, neglecting duration of transient process intervals. After such discretization in many practical applications the system of differential equations describes the dynamics of changes in the original complex system over the time interval from $T_k+\delta$ to $T_{k+1}$:

$$\dot{\overline{x_i}}={\overline{g}}_i(\overline{X},u,w,{\theta }_{s_k}),i=1,2,\dots ,m_{s_k},$$

(6)

where ${\overline{x}}_i$ is aggregated state $\left\{x_i\right\}$ or $\left\{x_{\gamma }\right\}$ from cluster $X_i$, $\overline{X}=col({\overline{x}}_1,{\overline{x}}_2,\dots ,{\overline{x}}_{m_{s_k}})$, ${\theta }_{s_k}$ is a finite set of current parameters at time interval $[T_k+\delta ,T_{k+1})$.

3. Control Problem and Discretization on Time

Consider a control problem of choosing the strategy of control minimizing the cost function comprised of local cost functions computed at different parts of the system:

$$L\left(\left\{u\right\}\right)=\sum _{i\in \mathbb{M}}l(x_i,u)\to \underset{u}{min}$$

or, in case the system consists of continuum elements (3):

$$L\left(\left\{u\right\}\right)={\int }_{\mathbb{M}}l(x_{\gamma },u)d\gamma \to \underset{u}{min}.$$

Previously it was assumed perturbation W has a “finite structure” at each moment, and the structure of $s_k$ changes at times $T_0,T_1,T_2,\dots$, causing clustering of the state space:

$${\mathcal{X}}_{s_k}=\{X_1,X_2,\dots ,X_{m_{s_k}}\}:X={\cup }_{i=1,2,\dots ,m_{s_k}}{X}_i,X_i\subset X$$

(7)

where ${\overline{x}}_i$ is aggregated state $\left\{x_i\right\}$ or $\left\{x_{\gamma }\right\}$ of cluster $X_i$, $\overline{X}=col({\overline{x}}_1,{\overline{x}}_2,\dots ,{\overline{x}}_{m_{s_k}})$, ${\theta }_{s_k}$ is finite set of current parameters. We assume that for any k the dimension of ${\theta }_{s_k}$ is bounded by d. For $t\in [T_k+\delta ,T_{k+1})$ due to integral (sum) additive property the loss function could be changed to:

$$L={\int }_{\mathbb{M}}l(x_{\gamma },u)d{x}_{\gamma }\approx \sum _i^{m_{s_k}}{\overline{l}}_k({\overline{x}}_k,u,{\theta }_{s_k})={\overline{L}}_k(\overline{X},u,{\theta }_{s_k}).$$

Assume that control strategy u are piece-wise constant and changing in the end of each time interval of length h,

$$u\left(t\right)=u_n,t\in [(n-1)h,nh).$$

The feedback u will be computed on the base of noised observations of the loss function $\tilde{L}$. After sampling in time and space, we obtain an observation model for the loss function:

$$y_n={\tilde{L}}_k({\overline{X}}_n,u_n,{\theta }_{s_k})+{\xi }_n,$$

(8)

where $t_n\in [T_K+\delta ,T_{k+1})$, ${\tilde{L}}_k(·)$ are functions from ${\overline{X}}_n=\overline{X}\left(t_n\right),u_n=u\left(t_n\right),{\theta }_{s_k},$ ${\xi }_n={\xi }_n^{\prime }+{\xi }_n\left(s_k\right)"$ is discrepancy (error) composed of some random noise ${\xi }_n^{\prime }$ independent of current system structure $s_k$, and systematic error ${\xi }_n\left(s_k\right)"$ which is, in general, is some function of the current system state.

4. Parameter Identification Problem

Estimation of system parameters values can be formulated as an optimization problem. The discrepancy between estimated and real parameter values could be expressed in terms of some loss function value and thus to solve the system identification problem one has to minimize given loss function. To identify the system structure it is required to choose such control formation strategy $\left\{u\right\}$, which minimizes some loss function.

Suppose that for a given parameter vector ${\theta }_{s_k}\in {\mathbb{R}}^d$ the optimal control strategy

$$\left\{u_n\right\}=\mathcal{U}\left({\theta }_{s_k}\right)$$

is known. After substituting it in (8), we get the problem of minimizing the function

$$f_{s_k}\left(\theta \right)={\tilde{L}}_k({\overline{X}}_n,\mathcal{U}\left(\theta \right),{\theta }_{s_k})$$

when observing its values against under the influence of noise ${\xi }_n$. Under the assumptions made, the minimum of the function $f_{s_k}\left(\theta \right)$ is reached when $\theta ={\theta }_{s_k}$. To solve the formulated problem, the method from [5] can be used.

The distributed optimization task can be formulated in terms of finding the method of constructing cluster (meso-) control, in which the same control action is applied to all elements of the cluster. In this case, the discretized loss function can be represented as a distributed functional

$$f_{s_k}\left(\theta \right)=\sum _{i=1}^{m_{s_k}}{\tilde{l}}_k({\overline{X}}_n^i,{\mathcal{U}}^i\left(\theta \right),{\theta }_{s_k}).$$

This problem could be solved by appliying the stochastic optmimization type method from [15].

The considered problem setting is a particular case of more general problem of minimizing a non-stationary differentiable function with respect to $\theta$. Let ${\mathcal{F}}_{n-1}$ be the $\sigma$-algebra of all probabilistic events which happened during time interval $n-1$ before start of time interval n. Hereinafter ${\mathbb{E}}_{{\mathcal{F}}_{n-1}}$ is a symbol of the conditional mathematical expectation with respect to the $\sigma$-algebra ${\mathcal{F}}_{n-1}$, $\mathbb{E}$ is a symbol of the mathematical expectation. The minimum point ${\theta }_{s_k}$ of function

$$F_n\left(\theta \right)={\mathbb{E}}_{{\mathcal{F}}_{n-1}}{f}_{s_k}\left(\theta \right)\to \underset{\theta }{min}$$

needs to be estimated.

More precisely, using the observations $y_1,y_2,\dots ,y_n$ and inputs ${\theta }_1,{\theta }_2,\dots ,{\theta }_n$, construct an estimate ${\widehat{\theta }}_n$ of an unknown vector ${\theta }_{s_k}$ minimizing the time-varying mean-risk functional.

4.1. Assumptions

Let us formulate Assumptions about disturbances and functions $f_{s_k}\left(\theta \right),F_n\left(\theta \right)$.

• For $n=1,2,\dots$, the successive differences ${\overline{\xi }}_n={\xi }_n^{+}-{\xi }_n^{-}$ of observation noise are bounded: $|{\overline{\xi }}_n|\le c_{\xi }<\infty$, or $\mathbb{E}{\overline{\xi }}_n^2\le c_{\xi }^2$ if a sequence $\left\{{\xi }_n\right\}$ is random, where ${\xi }_n^{+}$, ${\xi }_n^{-}$ are observation noises occurred during the same time interval n but at different time instants.
• Functions $F_n(·)$ have unique minimum points ${\theta }_{s_k}$ and $\forall \theta \langle \theta -{\theta }_{s_k},{\mathbb{E}}_{{\mathcal{F}}_{n-1}}\nabla f_{s_k}\left(\theta \right)\rangle \ge \mu {\parallel \theta -{\theta }_{s_k}\parallel }^2$ with a constant $\mu >0$. Here and further $\langle ·,·\rangle$ is a scalar product of two vectors.
• The gradient $\nabla f_{s_k}$ is uniformly bounded in the mean-squared sense at the minimum points ${\theta }_t:E{\parallel \nabla f_{s_k}\left({\theta }_t\right)\parallel }^2\le g^2$
• $\forall s_k\in S$ the gradient $\nabla f_{s_k}\left(\theta \right)$ satisfies the Lipschitz condition: $\forall {\theta }^{\prime },{\theta }^{\prime \prime }$

  $$\parallel \nabla f_{s_k}\left({\theta }^{\prime }\right)-\nabla f_{s_k}\left({\theta }^{\prime \prime }\right)\parallel \le M\parallel {\theta }^{\prime }-{\theta }^{\prime \prime }\parallel$$

  with a constant $M\ge \mu$.
• $\forall n\ge 1$ random vector ${\Delta }_n$ does not depend on ${\overline{w}}_n$, random vectors ${\overline{w}}_n,{\Delta }_n$ do not depend on ${\overline{w}}_1,\dots ,{\overline{w}}_{n-1}$; if $\left\{{\overline{v}}_n\right\}$ are random variables, then ${\overline{w}}_n,{\Delta }_n$ also do not depend on ${\overline{v}}_1,\dots ,{\overline{v}}_n$.

Using available observations, it is necessary to construct a sequence of estimates ${\widehat{\theta }}_n$ of the unknown vector ${\theta }_{s_k}$ that minimizes the function $f_{s_k}\left(\theta \right)$. To solve the problem, we will use an iterative algorithm with two measurements.

Let the trial simultaneous disturbance ${\Delta }_n,n=1,2,\dots$ be an observable (set or user-controlled) sequence of independent random vectors with known distribution functions ${\mathrm{P}}_n(·)$ — and specified vector functions ${\mathcal{K}}_n(·):{\mathbb{R}}^d\times {\mathbb{R}}^d\to {\mathbb{R}}^1$, satisfy the conditions

$$\int {\mathcal{K}}_n\left(x\right){\mathrm{P}}_n\left(dx\right)=0;\int {\mathcal{K}}_n\left(x\right)x^{\mathrm{T}}{\mathrm{P}}_n\left(dx\right)=\mathrm{I},$$

(9)

$$\underset{n}{sup}\int {\parallel {\mathcal{K}}_n\left(x\right)\parallel }^2{\mathrm{P}}_n\left(dx\right)<\infty ,n=1,2,\dots .$$

(10)

4.2. Algorithm

Let’s choose an arbitrary initial estimate vector ${\widehat{\theta }}_0\in {\mathbb{R}}^d$ and scalar parameters $\alpha$, $\beta$ for an iterative algorithm

$$\left\{\begin{array}{cc}\hfill & {\theta }_n^{+}={\widehat{\theta }}_{n-1}+\beta {\Delta }_n,{\theta }_n^{-}={\widehat{\theta }}_{n-1}-\beta {\Delta }_n\hfill \\ \hfill & y_n^{+}=f_n\left({\theta }_n^{+}\right)+{\xi }_n^{+},y_n^{-}=f_n\left({\theta }_n^{-}\right)+{\xi }_n^{-}\hfill \\ \hfill & {\widehat{\theta }}_n={\widehat{\theta }}_{n-1}-\frac{\alpha }{2\beta }{\mathcal{K}}_n\left({\Delta }_n\right)\left(y_n^{+}-y_n^{-}\right).\hfill \end{array}\right.$$

(11)

4.3. Main Result

Let us introduce the following notation:

$$\nu =\alpha \left(0.5-\mu -\frac{L\alpha C_{\tau }}{2}\right),\varphi =\alpha \gamma +\frac{L}{2}{\alpha }^2{C}_{\sigma },\gamma =0.5{\left(C_1\beta \right)}^2,\psi =\frac{\varphi }{\nu }$$

Here $C_{\tau }$, $C_{\sigma }$ could be chosen to satisfy inequation $\mathbb{E}\parallel \frac{1}{2\beta }(f_n(\theta +\beta {\Delta }_n)+{\xi }_n^{+}-f_n(\theta -\beta {\Delta }_n)-{\xi }_n^{-}){\parallel }^2\le C_{\sigma }+C_{\tau }{\parallel \theta -{\theta }_{s_k}\parallel }^2$ and $C_1\ge \int \parallel {\mathcal{K}}_n{\left(x\right)\parallel M\parallel x\parallel }^2{\mathrm{P}}_n\left(dx\right)$ and $\alpha$ chosen to satisfy following conditions:

$$\begin{array}{cc}\hfill & 0\le \alpha \le \frac{4\mu -2}{L{C}_{\tau }},\alpha \le \frac{2\mu -1-\sqrt{{(2\mu -1)}^2-2L{C}_{\tau }}}{L{C}_{\tau }}\hfill \\ \hfill & \alpha \ge \frac{2\mu -1+\sqrt{{(2\mu -1)}^2-2L{C}_{\tau }}}{L{C}_{\tau }}.\hfill \end{array}$$

(12)

 Theorem 1. 

Let Assumptions 1–5 and conditions for kernels $\mathcal{K}$ (9)-(10) and α (12) be satisfied. Set ${\widehat{\theta }}_0$, choose interval size parameter k

$$\mathbb{E}\{\parallel {\widehat{\theta }}_n-{\theta }_{s_k}{\parallel }^2\}\le \mathbb{E}\{\parallel {\widehat{\theta }}_0-{\theta }_{s_k}{\parallel }^2\}{(1-{\nu }_i)}^n+\psi (1-{(1-{\nu }_i)}^n).$$

(13)

Proof. To analyze the convergence of the algorithm (11) estimates, a method from [10] is used. Choose

$$V\left({\widehat{\theta }}_n\right)=\frac{1}{2}{\parallel {\widehat{\theta }}_n-{\theta }_{s_k}\parallel }^2$$

(14)

as Lyapunov function. To prove the theorem (13) is true it is sufficient to show the following six propositions are satisfied.

Proposition 1.  The iterative process $Y_n$, which defines the direction of the estimate change, is a Markov process, i.e. the distribution of the random vector $Y_n$ depends only on ${\widehat{\theta }}_n$ and n:

$$Y_n=\frac{1}{2\beta }\mathcal{K}\left({\Delta }_n\right)(y_n^{+}-y_n^{-})$$

(15)

For algorithm (11) we have

$$\mathbb{E}\left\{Y_n\right\}=\mathbb{E}\left\{{\mathcal{K}}_n\left({\Delta }_n\right)\frac{y_n^{+}-y_n^{-}}{2\beta }\right|{\mathcal{F}}_{n-1}\},$$

where the right hand side depends only on ${\widehat{\theta }}_n$ and n in the sense that ${\Delta }_n$ does not depend on any other random variables. Proposition is true.

Proposition 2.  $V\left({\widehat{\theta }}_n\right)\ge 0$, $infV\left({\widehat{\theta }}_n\right)=0$, $V\left({\widehat{\theta }}_n\right)$ has first-order derivative, and its gradient satisfies Lipshitz condition:

$$\parallel \nabla V\left(x\right)-\nabla V\left(\theta \right)\parallel \le L\parallel x-\theta \parallel \forall x,\theta \in {\mathbb{R}}^d.$$

The proposition is valid due to the choice of Lyapunov function (14).

Proposition 3.   Pseudo-gradient condition:

$$\langle \nabla V({\widehat{\theta }}_n),\mathbb{E}\left\{Y_n\right\}\rangle \ge {\delta }_{n}V\left({\widehat{\theta }}_n\right)-{\gamma }_n,{\delta }_n>0,{\gamma }_n\ge 0.$$

(16)

At first consider $\mathbb{E}\left\{Y_n\right\}$. Due to (11) pseudo-gradient (15) after using Assumption 5 becomes

$$\mathbb{E}\left\{Y_n\right\}=\mathbb{E}\left\{{\mathcal{K}}_n\left({\Delta }_n\right)\frac{1}{2\beta }\left(f_{s_k}\left({\theta }_n^{+}\right)-f_{s_k}\left({\theta }_n^{-}\right)\right)\right|{\mathcal{F}}_{n-1}\}$$

(17)

Consider the expression under mathematical expectation. Using Taylor series representation it could be written as:

$$\begin{array}{cc}\hfill & \frac{1}{2}\nabla F_n\left({\widehat{\theta }}_n\right)+\frac{1}{2\beta }\int {\mathcal{K}}_n\left(x\right)x^T{\int }_0^1\left({\nabla }_x{f}_{s_k}({\widehat{\theta }}_n+t\beta x)-\nabla f_{s_k}\left({\widehat{\theta }}_n\right)\right)dt{P}_n\left(dx\right)+\hfill \\ \hfill & +\frac{1}{2}\nabla F_n\left({\widehat{\theta }}_n\right)-\frac{1}{2\beta }\int {\mathcal{K}}_n\left(x\right)x^T{\int }_0^1(\nabla f_{s_k}({\widehat{\theta }}_n-t\beta x)-\nabla f_{s_k}\left({\widehat{\theta }}_n\right))dt{P}_n\left(dx\right)\hfill \end{array}$$

Estimate absolute value of the sum of integral elements in the obtained expression. After using (9), Assumption 4 we get

$$|\int (·)P_n\left(dx\right)|+|\int (·)P_n\left(dx\right)|\le \frac{2{\beta }^2}{2\beta }\int \parallel {\mathcal{K}}_n\left(x\right)\parallel \parallel x\parallel M\parallel x\parallel P_n\left(dx\right)\le C_1\beta .$$

(18)

Substitute the estimate, elements containing gradients of $F_n$ and (14) into (16), regard the relation $\int (·)P_n\left(dx\right)\ge -|\int (·)P_n\left(dx\right)|$ and Cauchy–Bunyakovsky–Schwarz inequality:

$$\langle {\widehat{\theta }}_n-{\theta }_{s_k},\mathbb{E}\left\{Y_n\right\}\rangle \ge \langle {\widehat{\theta }}_n-{\theta }_{s_k},\nabla F_n\left({\widehat{\theta }}_n\right)\rangle -\parallel {\widehat{\theta }}_n-{\theta }_{s_k}\parallel C_1\beta .$$

Apply Assumption 2 and estimate $\parallel {\widehat{\theta }}_n-{\theta }_{s_k}\parallel C_1\beta \le \frac{1}{2}(\parallel {\widehat{\theta }}_n-{\theta }_{s_k}{\parallel }^2+{\left(C_1\beta \right)}^2):$

$$\langle {\widehat{\theta }}_n-{\theta }_{s_k},\mathbb{E}\left\{Y_n\right\}\rangle \ge (2\mu -1)\frac{1}{2}{\parallel {\widehat{\theta }}_n-{\theta }_{s_k}\parallel }^2-\frac{1}{2}{C}_1^2{\beta }^2$$

Proposition is true for $\mu >1/2$.

Proposition 4.

$$\mathbb{E}\{\parallel Y_n{\parallel }^2\}\le {\sigma }_n^2+\tau V\left(x\right),{\sigma }_n\ge 0,{\tau }_n\ge 0.$$

Using (10), Cauchy–Bunyakovsky–Schwarz inequality, Assumptions 1 and 3 it could be shown that

$$\begin{array}{cc}\hfill \mathbb{E}\{\parallel Y_n{\parallel }^2\}& \le \frac{1}{2}\underset{x}{sup}{\mathcal{K}}_n{\left(x\right)}^2\mathbb{E}\{{\left({\xi }_n^{+}\right)}^2+{\left({\xi }_n^{-}\right)}^2|{\mathcal{F}}_{n-1}\}+\hfill \\ \hfill & +\int {\left(f_{s_k}\left({\widehat{\theta }}_n^{+}\right)-f_{s_k}\left({\widehat{\theta }}_n^{-}\right)\right)}^2\parallel {\mathcal{K}}_n{\left(x\right)\parallel }^2{\mathrm{P}}_n\left(dx\right)\le C_2{\beta }^2(\parallel {\widehat{\theta }}_{n-1}-{\theta }_{s_k}{\parallel }^2)+C_3{\beta }^4+C_4{\xi }_n^2\hfill \end{array}$$

Proposition holds true with $\tau =2C_2{\beta }^2$, ${\sigma }_n^2=C_3{\beta }^4+C_4{\xi }_n^2$.

Proposition 5.   $\mathbb{E}V\left({\widehat{\theta }}_0\right)<\infty .$

Proposition is valid due to arbitrariness of initial approximation ${\widehat{\theta }}_0$ choice and an assumption regarding final order $s_k$ of the external disturbance W affecting the system and thus the final order of the system state vector.

Proposition 6.   $0\le \nu \le 1$; ${\sum }_n\nu =\infty ,n\to \infty$

The first inequality could be met by choice of α and the second one is true since ν is constant. Proposition is true.

The fulfillment of the given propositions allow to prove the theorem using the result in [10]. □

 Remark 1. 

After the algorithm converges the parameter estimates ${\widehat{\theta }}_n$ continue to fluctuate around the true parameter value ${\theta }_{s_k}$ until the new change of system structure and value of ${\theta }_{s_k}$.

 Remark 2. 

The most common problem setup for unknown-but-bound disturbances in existing works is formulated as follows. It is required to minimize the value of the objective function $f\left(\theta \right)$ with some adversarial deterministic noise $\xi \left(\theta \right)$ such that $\left|\xi \right(\theta \left)\right|\le \xi$ and $\xi >0$):

$$\tilde{f}\left(\theta \right)=f\left(\theta \right)+\xi \left(\theta \right).$$

The considered setup implies that external noise depends on system state $s_k$ which is not a direct function of θ but rather is affected by values of ${\theta }_t$ during some time period $[t-d,t)$

$$\tilde{f}\left(\theta \right)=f\left(\theta \right)+\xi \left(s_k\right).$$

Value $s_k$ is formed on basis of previous values of θ.

5. Application for Orientation Improvement of Radio Telescope’s Elements

For astronomers, an urgent task was to obtain images of objects that can’t be recorded using optical telescopes, which are situated on earth or in space. This problem was largely solved using radio telescopes [11]. The main tasks of the telescope are: to collect radiation that falls on the mirror system with minimal losses, and also obtain the most accurate image of the object [6]. There are various ways to solve this problem, consider for example, [8]. One of them is improving the quality of the device [12], which collects radiation for obtaining images. Another one is combining radio telescopes into systems [9]. If the radiation is collected with significant errors, then the image will have disturbances [27].

An important and time-consuming part of image acquisition is the precise tuning of the radio telescope antenna (or systems consisting of such antennas) [13]. The quality of the image obtained on a radio telescope directly depends on the quality of the construction of a reflecting system of mirrors that focuses the radiation coming from outside. To improve the image quality, it is necessary to focus the radiation of the device in such a way that it works as accurately as possible, especially if it is located in space [14]. Traditional antenna tuning algorithms are sufficient for the task. However, they lose their effectiveness under uncontrolled unpredictable external influences. We consider the case when these are deformations of the radio telescope shields that arise due to environmental influences such as temperature changes, wind and other influences.

In practice, radiation is subject to various distorting influences, and as a result the quality of the observed image decreases, despite the presence of various stabilizers and filters [11]. One of the ways to solve such a problem is using randomized stochastic optimization algorithms [15]. A method for improving image quality by improving the tuning characteristics of the radio telescope mirror system is considered in [25].

In an ideal antenna system, the signal is reflected from different mirror plates and assembled at one point. Deformations of radio telescope structures, external temperature, wind, mechanical influences lead to deviation of the optical path lengths of the rays from the required ones. As a result, the focus point on the plane of the receiver shifts [12]. To improve reflection accuracy on the surface of radio telescope mirrors, the following methods are used: autocollimation [22], telescope calibration by spectral density radiation flux, synchronous calibration method [21], laser geodetic measurements [6], improvement of the kinematics of antenna elements [19], radio holography physical method [24]. The goal is to develop a stochastic optimization type algorithm to improve the quality of the settings of the radio telescope mirror system model which could be used in a system similar to Radioastron [20,26]. The main criterion for efficiency is the recording power of the desired signal and the time required for adjusting the parameters.

The antenna segments can be set to the optimal position to improve the quality of image recording. Consider an irradiator (radiation generator), a receiver and a mirror system of a radio telescope, consisting of identical plates that reflect the incoming signal ($i=1,\dots ,N,N=895$ in the RATAN-600 installation). Radiation is created in the irradiator, which, falling on the plates and focused in the receiver. Let’s divide the number axis into time intervals of duration $\delta$, starting from some moment $t_0$, where $T_k$ is the k-th time interval k is the index of the time interval. Assume that we know:

1) the position (orientation) of each i-th plate, which is specified by the vector of parameters ${(a_i,b_i,c_i)}^T$, where $a_i$ is the rotation angle of the i-th reflective element horizontally; $b_i$ is the vertical rotation angle, $c_i$ is the forward horizontal displacement of the i-th reflecting element. Let ${\theta }_k$ be a vector that contains all parameters of the mirror system in a given time interval k, ${\theta }_k={(a_1,a_2,...,a_N,b_1,b_2,...,b_N,c_1,c_2,...,c_N)}^T,$

2) radiation coming from each mirror ($z_i$),

3) the common signal coming from all mirrors to the receiver $Z\left({\theta }_k\right)={\sum }_{i=1,\dots ,N}\left(z_i\right)$,

4) characteristics of the signal in the generator.

The front of a signal is the sum of harmonics that comes to us with different phases from a certain direction. The perfect placement of the plates brings all radiation from the objects into focus. We obtain the signal as a sum of sines with phases. Different $z_i$ arrive at different times with different phases. Let us evaluate the difference between the signals from an ideal antenna and from a real one (with deformations). Signals reflected from ideal mirror segments will have the same phase ${\varphi }_i$ (${\varphi }_1={\varphi }_2=\cdots ={\varphi }_N$). The signals from segments with deformations will look like this:

$$z_1=sin(\omega t+{\varphi }_1),z_2=sin(\omega t+{\varphi }_2),\cdots z_i=sin(\omega t+{\varphi }_i)$$

(19)

The objective function of the problem $\left(F\right(\theta )$ is the signal power) is defined as follows:

$$F\left(\theta \right)=\overline{{lim}_k}{\int }_{t\in T_k}{\left|Z(\theta ,t)\right|}^{2}dt.$$

(20)

It is required to maximize the objective function:

$$F\left(\theta \right)\to \underset{\theta }{max}$$

(21)

We consider the problem of optimizing the position of mirrors in the limit over time (not over a specific time interval).

The reflective elements of the antenna are made exactly the same, so they provide equivalent observations in all directions.

At the same time, if one moves along an ideal reflective surface, its local characteristics change. Therefore, a real reflective surface composed of identical elements will repeat deviations from the ideal surface from element to element. These deviations will be greater, if the shape of the surface of the element differs from the shape of that portion of the ideal surface which this element should represent. If the size of the reflective elements increases, then the deviations naturally increase. These deviations are an error distributed over the reflector of a variable profile antenna, which, at large values, creates unacceptable distortions reducing the efficiency of the antenna. We will call $v_k$ the signal power measurement errors arising due to unknown and uncontrolled deformations in the reflective elements of the antenna caused by weather, wind, and temperature changes.

After conventional procedures for adjusting the inclination angles and positions of the radio telescope mirrors, we obtain an initial approximation ${\widehat{\theta }}_0$ to the optimal tuning values. Then there is still the possibility of “tuning” in a certain neighborhood T containing ${\widehat{\theta }}_0$ and the optimal value ${\theta }^{★}$ corresponding to the maximum power of the received signal. We use a stochastic optimization algorithm with two measurements per iteration, which allows us to reduce the negative impact of various disturbances on the power of the recorded signal [7]:

Select ${\theta }_0.$

$n-1\to n$

We sequentially generate vectors ${\Delta }_n$ with components from +/-1, which are chosen with equal probability.

We measure the power values for two positions of the antenna system:

$({\theta }_n+\beta {\Delta }_n)$ and $({\theta }_n-\beta {\Delta }_n).$

Measurements are obtained with noise $v_{2n}$ and $v_{2n-1}$ :

$$y_{2n}=P_{2n}({\theta }_{2n}+\beta {\Delta }_n)+{\xi }_{2n};$$

(22)

$$y_{2n-1}=P_{2n-1}({\theta }_{2n-1}-\beta {\Delta }_n)+{\xi }_{2n-1};$$

(23)

Next, we form the following estimate ${\theta }_{2n}$ according to the rule:

$${\theta }_{2n}={\mathcal{P}}_{\mathcal{T}}({\theta }_{2n-1}+\frac{\alpha K_n\left({\Delta }_n\right)}{2\beta }(y_{2n}-y_{2n-1})),$$

(24)

where $\alpha ,\beta$ are the parameters of the algorithm, ${\mathcal{P}}_{\mathcal{T}}$ is the projection onto the set $\mathcal{T}$.

To construct the kernels $K_0(·)$ and $K_1(·)$ on the interval $[-1/2,1/2]$ orthogonal Legendre polynomials could be used. In this case, for initial values $\ell =1,2$ (i.e. $2\le \gamma \le 3$) the type of kernels is as follows: $K_0\left(q\right)=12q,K_1\left(q\right)=1,\mid q\mid \le 1/2,$

for $\ell =3,4$ (i.e. $3<\gamma \le 5$):

$$\begin{gathered} K_0\left(q\right)=5q(15-84q^2), \\ {K}_1\left(q\right)=9/4-15q^2,\mid q\mid \le 1/2, \end{gathered}$$

and for $\left|q\right|>1/2$ both functions are equal to zero.

The test disturbance is formed in such a way that $\forall n\ge 1$ random vector ${\Delta }_n$ does not depend on ${\overline{v}}_1,\dots ,{\overline{v}}_k$ and $\mathrm{E}\{{(v_{2n}-v_{2n-1})}^2/2\}\le {\sigma }_2^2,(\mathrm{E}\left\{v_n^2\right\}\le {\sigma }_1^2);$

The main element that reflects the incoming signal is the segment of the mirror system. It is important to configure these shields so that the signal comes into focus. The system allows to customize the required dimensions of the reflective shield and its curvature. Stochastic optimisation algorithm (11) can be applied for tuning a system which consists of a large number of telescopes in space under conditions of interference and also when individual elements of the system are deformed. About 3000 parameters of the mirror surfaces are adjusted, during focusing of the mirror system.The faster setting up of the mirrors is very important for the accuracy of observations. For acceleration, it is proposed to use the Nesterov acceleration method in distributed form [16,17,18]. In Figure 3 on page 11, the dependence of the signal power on the number of iterations of the algorithm is given.