IDENTIFICATION OF A STATIC OBJECT USING INTEGRAL PROJECTIONS

V. V. Zakharov

(Presented by the scientific seminar of the Laboratory of Computer Technology and Automation, Research Institute of Nuclear Physics)

For the construction of mathematical models (identification) of static objects, the method of regression analysis [1] and various modifications of the method of optimal experimental design [2] are currently widely used. In [3], by means of a technique called coordinatewise integral projection, a nonlinear object was partially identified (the coefficients at the second powers of the variables were obtained). This same method gives a very simple algorithm for constructing a linear model and is easily generalized to the case of a nonlinear dependence of the variables.

1. Basic idea of the method. Suppose it is required to find the coefficients of the linear form

\[ F(\vec{x})=\sum_{i=1}^{n} b_i x_i . \tag{1.1} \]

from the results of observations \(F_j=F(\vec{x}^{\,j})\), where \(\vec{x}^{\,j}=(x_1^j,\ldots,x_n^j)\),

\[ \vec{x}^{\,j}\in[\alpha,\beta],\quad \alpha=(\alpha_1,\ldots,\alpha_n),\quad \beta=(\beta_1,\ldots,\beta_n),\quad j\in J,\quad J=\{1,2,\ldots,N\}. \]

Denote

\[ I(\alpha,\beta)=\int_{\alpha}^{\beta} F(\vec{x})\,d\vec{x}. \]

Obviously,

\[ I(\alpha,\beta)=\prod_{i=1}^{n}(\beta_i-\alpha_i)\sum_{i=1}^{n} b_i\frac{\alpha_i+\beta_i}{2}. \tag{1.2} \]

Rewrite (1.2) in another form and call this form the integral projection onto the \(k\)-th coordinate axis:

\[ I(\alpha_k,\beta_k)=\prod_{i\ne k}(\beta_i-\alpha_i)(\beta_k-\alpha_k) \left(\sum_{i\ne k} b_i\frac{\alpha_i+\beta_i}{2} +\frac{\alpha_k+\beta_k}{2}\right). \]

Compute \(I(\alpha_k,\alpha_k+\delta_k)\), \(I(\beta_k-\delta_k,\beta_k)\), where \(\delta_k=(\beta_k-\alpha_k)/(t+1)\), \(t\ge 1\).

The difference of these quantities gives

\[ I(\beta_k-\delta_k,\beta_k)-I(\alpha_k,\alpha_k+\delta_k) = \delta_k^2 t b_k \prod_{i\ne k}(\beta_i-\alpha_i), \]

whence

\[ b_k= \frac{I(\beta_k-\delta_k,\beta_k)-I(\alpha_k,\alpha_k+\delta_k)} {\delta_k^2 t \prod_{i\ne k}(\beta_i-\alpha_i)} \quad k=1,2,\ldots,n. \tag{1.3} \]

Thus, in order to compute the coefficients \(b_k\), it is necessary to find an effective method for computing the integrals in (1.3). The Monte Carlo method [4] is very convenient in this respect.

Let \(I(a,\beta)\) be a statistical estimate of \(J(\alpha,\beta)\). Then we have

\[ I(a,\beta)=\prod_{i=1}^{n}(\beta_i-\alpha_i)\frac{1}{N}\sum_{j=1}^{n}F_j, \tag{1.4} \]

whereas

\[ \left|I(a,\beta)-I(a,\beta)\right|\le c_1\sqrt{\frac{\sigma_F^2}{N}}, \tag{1.5} \]

where \(\sigma_F^2\) is the variance of \(F(\vec{x})\), a function of the uniformly distributed random variables \(x_k\) \((k=1,\ldots,n)\), and the quantity \(C_1\) is determined by the desired probability that inequality (1.5) will hold.

Substituting (1.4) into (1.3) and denoting \(t\delta_k\) by \(\tau_k\), we obtain for \(b_k\)

\[ b_k \simeq b_k^0= \frac{1}{\tau_k} \left( \frac{1}{N_{2k}}\sum_{s\in J_{2k}}F_s - \frac{1}{N_{1k}}\sum_{r\in J_{1k}}F_r \right). \tag{1.6} \]

Here \(r\in J_{1k}\) if \(x_k^r\in(\alpha_k,\alpha_k+\delta_k)\); \(s\in J_{2k}\) if \(x_k^s\in(\beta_k-\delta_k,\beta_k)\). It is clear that the set \(J_{1k}\) contains \(N_{1k}\) elements, and the set \(J_{2k}\)—\(N_{2k}\) elements, with \(J_{1k}\cup J_{2k}=J\) and \(N_{1k}+N_{2k}=N\).

Of special note is the fact that, in order to compute all the coefficients of the linear form (1.1) by formula (1.6), one and the same \(N\) observations \(F_j\) are used, only projected onto different coordinate axes.

2. Refinement of the coefficients. The variance \(\sigma_F^2\), and hence the error of formula (1.6), may in particular cases be rather large. The problem arises of reducing the value of \(\sigma_F^2\) and thereby the error in computing \(b_k\), without increasing the number of observations. Such a possibility is provided by the method of extracting the principal part of the function integrated by Monte Carlo, described in [4].

As the principal part of \(F(x)\), consider the function

\[ Y^0(\vec{x})=\sum_i b_i^0 x_i. \tag{2.1} \]

The expression

\[ I^0(a,\beta)= \prod_i(\beta_i-\alpha_i)\frac{1}{N}\sum_{j=1}^{N} \left[F_j-Y^0(\vec{x}^{\,j})\right] + \int_{\alpha}^{\beta}Y^0(\vec{x})\,d\vec{x}. \tag{2.2} \]

gives a more accurate value of \(I(\alpha,\beta)\) than (1.4). Substituting this approximate value into (1.3) instead of the corresponding exact value, we obtain a simple iterative formula for computing the coefficients \(b_k\):

\[ b_k \approx b_k^m=b_k^0+b_k^{m-1}-\frac{1}{\tau_k} \left( \frac{1}{N_{2k}}\sum_s Y_s^{m-1} - \frac{1}{N_{1k}}\sum_r Y_r^{m-1} \right),\quad m=1,2,\ldots \tag{2.3} \]

3. Influence of noise. Estimation of the variances of the coefficients. It can be shown that if the observed values \(F_j\) are distorted by superposed additive noise with variance \(\sigma^2\), then the coefficients \(b_k\) will be computed with variances

\[ \sigma^2(b_k) \leq c_2 \frac{4}{\tau^2}\frac{\sigma^2}{N}, \tag{3.1} \]

where \(c_2\) is a constant characterizing the probability that inequality (3.1) holds.

These variances have the same order as in linear orthogonal planning [2]. Consequently, the use of formula (2.3) for computing the gradient vector of an arbitrary function \(F(x)\) in the presence of noise is promising.

4. Generalization of the form of the approximating polynomial. Note that

\[ \int_{-\delta}^{\delta}\prod_{i\ne k} x_i^{m_i}\cdot x_k^r\,d\vec{x}=0, \tag{4.1} \]

where \(m_i\) is a real quantity, and \(r\) is odd. (In (4.1) the integration is multidimensional.)

Suppose that the hypothesis \(F(\vec{x})\approx Q(\vec{x})\) has been proposed, where

\[ Q(\vec{x})=\sum_i a_{ii}x_i^2+\sum_{i\ne k}a_{ik}x_{ik}+\sum_i a_{i0}+a_{00} \]

and it is required to establish the form of \(Q(\vec{x})\) (i.e., to find \(\{a_{ik}\}, i,k=\overline{0,n}\)) from the results of observations \(\{x^j,F_j\},\ j\in J\).

Taking (4.1) into account, we obtain

\[ \int_{-\delta}^{\delta}Q(\vec{x})\,d\vec{x} = (2\delta)^n \left( \frac{\delta^2}{3}\sum_{i=1}^n a_{ii}+a_{00} \right), \tag{4.2} \]

whence we find \(a_{00}\), provided that \(a_{ii}\) \((i=\overline{1,n})\) have been found by the projection method described in [3].

Further,

\[ \int_{-\delta}^{\delta}Q(\vec{x})x_i\,d\vec{x} = (2\delta)^n a_{i0}\cdot \frac{1}{3}\delta^2. \tag{4.3} \]

\[ \int_{-\delta}^{\delta}Q(x)x_i x_k\,d\vec{x} = (2\delta)^n a_{ik}\cdot \frac{2\delta^4}{9},\quad i,k=\overline{1,n}. \]

Replacing formulas (4.2), (4.3) by their discrete analogue according to (1.4), we obtain estimates of \(a_{ik}\) on the basis of uniformly scattered observations \(F_j\):

\[ \frac{1}{N}\sum_{j=1}^{N}F_j \approx \frac{\delta^2}{3}\sum_{i=0}^{n}a_{ii}+a_{00}; \quad \frac{1}{N}\sum_{j=1}^{N}x_i^jF_j \approx a_{i0}\cdot\frac{\delta^2}{3}, \]

\[ \frac{1}{N}\sum_{j=1}^{N} x_i^j x_k^j F_j \approx a_{ik}\frac{2\delta^4}{9}; \quad i,\ k=\overline{1,n}. \tag{4,4} \]

Remark 1. It is not difficult to see that, by choosing in an appropriate manner the composition of the factors under the integral and their powers, one can identify polynomials of arbitrary degree.

Remark 2. The organization of iterative refinement of the coefficients \(a_{ik}\) \((i,k=\overline{0,n})\), using the procedure of extracting the principal part by formulas (2,2), presents no difficulty.

5. A more general form of the distribution of the variables.
Let the controlled variables \(x_k\) have an arbitrary (and not necessarily uniform) joint distribution density \(p(\vec{x})\), truncated by the hyperinterval \([\alpha,\beta]\).

Obviously,

\[ I(\alpha,\beta)=M\left[\frac{F(\vec{x})}{p(\vec{x})}\right] = \int_{\alpha}^{\beta} \frac{F(\vec{x})}{p(\vec{x})}\,p(\vec{x})\,d\vec{x}. \]

Therefore, everywhere in relations (1,6), (4,4), replacing the values of the function \(F_j\) by the ratio \(F_j/p(x^j)\), we obtain the possibility of using iterations of type (2,3) to find the coefficients \(a_{ik}\), without changing the formulas.

6. Field of application of the method.
a) In problems of active experiment, the proposed method makes it possible, on the basis of the same observations, to test successively linear, quadratic, cubic, etc., hypotheses concerning the mathematical description of the object. Thus, it turns out that a uniform scattering of tests is a universal and successively compositional experimental design.

b) The generalization in Sec. 5 makes it possible to apply the proposed method of integral identification to the determination of operating regimes of objects from the results of a passive experiment. In this case, knowledge of the distributions of the components of the vector \(x\) makes it possible to reduce the variance of the coefficients sought.

c) In some problems, refinement of the model coefficients from newly obtained observations of the process (model correction) is of essential importance. Such are problems of object control or problems of classification and pattern recognition (for example, in [5]). From formula (2,3) it is clear that the inclusion of new observations in the iteration process is possible at any time and does not cause any complication of the algorithm.

REFERENCES

1. Yu. V. Linnik. The Method of Least Squares and the Foundations of the Mathematical-Statistical Theory of Processing Observations. Moscow, Fizmatgiz, 1962.

2. V. V. Nalimov, N. A. Chernova. Statistical Methods of Planning Extreme Experiments. Moscow, “Nauka,” 1965.

3. V. V. Zakharov. The method of integral smoothing in multi-extremal and stochastic problems. Izv. AN SSSR, “Technical Cybernetics,” No. 4, 1970.

4. N. P. Buslenko et al. The Method of Statistical Tests. Moscow, Fizmatgiz, 1962.

5. E. S. Bozhanov, G. K. Krug. Classification of multifactor objects from experimental data. In: Experimental Design, Moscow, “Nauka,” 1966.