Vol. 277

1977

DIRECT MODELING OF A STOCHASTIC PROGRAMMING PROBLEM

V. V. Zakharov, V. G. Smirnova

(Presented by the scientific seminar of the Laboratory of Computing Technology and Automation, Research Institute of Nuclear Physics, Siberian Branch, Academy of Sciences)

In [1] a problem of optimal control of a static object was considered in the presence of constraints on the quality of control. It proved possible to formulate and solve this problem as a stochastic programming problem. Under certain assumptions, the latter reduces to a problem of linear deterministic programming. In the present work, for a more general case, a special method of solution is proposed by modeling the motion toward a constrained minimum under conditions of noise. At this stage of solving the problem, the method of integral smoothing described in [2] is used. The methodology of integral identification also makes it possible to construct a linear model of the object being optimized and to correct it over time according to the results of its normal operation. In principle, however, there will be nothing new here in comparison with the method of regression analysis usually applied for the same purposes. Therefore the question of identifying a specific object is not considered here. It is nevertheless worth noting that the use of essentially one and the same program for identifying the object and for finding the optimal conditions of its functioning is undoubtedly attractive when using computers with a small memory capacity.

1. Problem formulation

A feature of some automatic-control objects is the limited possibility of recording the controlled variables. This applies in particular to chemical production processes.

Let the object operate in a static mode, and let the purpose of control be to attain the minimum of some function \(F(\vec{x})\). At the same time, the object must not violate certain production regimes, specified by inequalities \(f_j(\vec{x}) \ge 0,\ j \in J,\ J=\{1,2,\ldots,m\}\). The problem of optimal control of such an object is formulated as follows: find

\[ \min_x \{F(\vec{x}) : f_j(\vec{x}) \ge 0,\ j \in J\}, \tag{1} \]

where

\(x=(x_1,\ldots,x)\) is the vector of variables by means of which it is possible to control the object (controlled variables).

A feature of optimal control of chemical-technological processes is the partial controllability of the variables \(x_i\). Thus, in the automation of cellulose production at the Barnaul chemical combine, the deviation of the actual value of \(x_i\) from the nominally established one reached 30–40%. Consequently, it is meaningful to speak of a probabilistic form of the formulation of problem (1):

\[ \operatorname{E}\left[\min_x\{F(\vec{x}),\mathrm{P}[f_j(\vec{x})\geq 0]\geq a_j,\; j\in J\}\right], \tag{2} \]

where

\(0<a_j<1\); \(\operatorname{E}\), \(\mathrm{P}\) are the symbols of mathematical expectation and probability, respectively.

The formulation of the problem in such a general form causes great difficulties for its numerical solution. Therefore one usually passes to the modified formulation:

\[ \min_x\{\operatorname{E}[F(\vec{x})]/\mathrm{P}[f_j(\vec{x})>0]\geq a_j,\; j\in J\}, \tag{3} \]

although the latter is, obviously, approximate in character. A survey of methods of stochastic programming is given in [3].

In the case where the form of the functions \(F\), \(f_j\) is linear, and the variables \(x_i\) are normally distributed, it has been possible to reduce (3) to a deterministic linear programming problem in a region narrower than the original one [1]:

\[ \min_{\operatorname{E}(\vec{x})=\vec{\bar{x}}} \left\{ \sum_{i=1}^{n} c_i \bar{x}_i \;\middle|\; \sum_{i=1}^{n} a_{ij}\bar{x}_i \geq \Phi^{-1}(a_j) \sqrt{\sum_{i=1}^{n} a_{ij}^{2}s_{x_i}^{2}} \right\}, \tag{4} \]

where

\(\Phi\) is the integral distribution function of the Gaussian distribution,

\(s_{x_i}^{2}\) is the variance of the random variable \(x_i\).

If the control variables \(x_i\) have arbitrary distributions (in the specific case of control of cellulose production, the distributions of \(x_i\) were specified by histograms), then the only way to solve (2), (3) remains the method of direct modeling.

2. Initial assumptions of the modeling

Let us consider one method of directly solving problem (2), (3), consisting in modeling the search process. In modeling the solution of a stochastic programming problem (for brevity we shall henceforth say “stochastic program”), the form in which the densities of the random variables \(x_i\) are specified is, to a considerable extent, immaterial (this may be either an analytic expression or histograms). The distribution density of the quantities \(f_j(\vec{x})\) (we shall denote it by \(\varphi_j(f)\)) need not be specified a priori, but may be obtained in the course of the modeling. The functions \(F(x)\), \(f_j(x)\) may have a nonlinear form.

In principle, with simulation one can obtain an arbitrarily accurate solution of the problem, albeit at the cost of an enormous amount of computation. The obvious fact, however, is that practical needs are hardly such that we should seek a solution with high accuracy. Indeed, we are considering the case where even the distribution densities \(\varphi_i(t)\) are specified so roughly that there is no point in seeking their analytic form.

There is one more circumstance which, in our case, contributes to a substantial reduction in the number of trials when implementing the search model. It consists in the fact that we may consider problem (2), i.e., a more general problem than (3), and consequently may count on a more economical solution. At the same time, in the search process there is no need to seek \(\mathrm E[\min F(\vec{x})]\), but one may confine oneself to seeking \(\min \mathrm E[F(\vec{x})]\) at each iteration. It may be hoped that the error arising thereby will be compensated for in subsequent iterations. The result of the solution at each iteration has no economic meaning for us, since, in the terminology of [3], in formulations (2), (3) we are dealing with a one-stage stochastic programming problem. Therefore, maintaining the adequacy of formulation (2) and finding \(\mathrm E\left[\min\limits_x F(\vec{x})\right]\) at a separate iteration makes sense only at the final stages of the search. Taking this circumstance into account, the scheme for simulating the problem in formulations (2), (3) was constructed.

Let us introduce some definitions and make remarks explaining the prerequisites and components of the simulation algorithm that will be given below.

A. Let the distributions of the components of the vector \(\vec{x}\) be given by histograms. We shall simulate the random variable \(x_i\) in accordance with the given histogram. We shall denote the center of trials at the \(k\)-th iteration of the simulation process by \(\vec{x}^{\,k}\), and realizations of the vectors by \(\vec{x}^{\,kl}\) \((l=1,2,\ldots,\nu_k)\). We shall call the point \(\vec{x}^{\,k}\) stochastically feasible if the inequality \(\mu_{kj}/\nu_k \ge \alpha_j\) is satisfied, where \(\nu_k\) is the total number of points \(\vec{x}^{\,kl}\), and \(\mu_{kj}\) is the number of points that did not violate the inequalities \(f_j(\vec{x}) \ge 0\).

B. Denote by \(Y_d\) the set of feasible points of the original deterministic program. We shall write this fact briefly in the form

\[ Y_d=\{\vec{x}\mid f_j(\vec{x})\ge 0,\ j\in J\}. \]

Using this notation, the set of stochastically feasible points can be denoted by

\[ Y_{st}=\{\vec{x}\mid \mathrm P[f_j(\vec{x})\ge 0]\ge \alpha_j,\ j\in J\}. \]

Introduce the set

\[ Y_d^{*}=\{\vec{x}\mid f_j(\vec{x})-\delta\ge 0,\ j\in J\}. \]

In [4] an analogous set is introduced in order to eliminate cycling of the process. In our case it is also necessary in order to reduce the number of checks for stochastic feasibility of points

\(\vec{x}^{\,k1}\), which will accelerate convergence of the process far from the boundary of the region \(Y_{st}\). The quantity \(\delta\) may be taken in advance as a multiple of the mean (or largest) variance of the variables \(x_i\), for example, \(\delta=3\sigma_{\max}\). One may also set

\[ \delta_1=\max_{j\in J} f_j(\vec{x}^{\,0}),\quad \delta_2=\min_{j\in J} f_j(\vec{x}^{\,0}),\quad \delta=(\delta_1+\delta_2)/2 . \]

B. Similarly to how this is done in [4], we shall consider a direction \(s^k\) at the point \(\vec{x}^{\,k}\) to be stochastically feasible if \(\vec{x}^{\,k}+\lambda \vec{s}^{\,k}\in Y_{st}\) and \(\lambda>0\). If, along this direction, the value of the function \(F(x)\) decreases, then it is stochastically suitable.

G. Let \(\Theta\) denote the method for modeling problem (2), (3) as a whole. Then the method \(\Theta\) will consist of a set of modeling rules \(\Theta_q\), \(Q\) in number, which we shall denote as follows: \(\Theta=(\Theta_1,\ldots,\Theta_Q)\). One may consider several methods \(\Theta\) intended for one or another kind of problem, for the density distribution of the random variables \(x_i\), and for variants of the formulations (2), (3). In the present work we shall consider two of them, which have a universal character.

We shall now consider the scheme for modeling problem (2) in general terms, and then detail it. The basis of the method \(\Theta\) is the rule \(\Theta_1\)—the choice of a direction of decrease of the objective function \(F(\vec{x})\) in the region \(Y_{st}\), i.e., a stochastically suitable direction. Motion along this direction must, sooner or later, make the point \(\vec{x}^{\,k+1}\) stochastically infeasible. For this case the action of rule \(\Theta_2\) is provided for—the “ricochet” from the boundary \(Y_{st}\) into the interior of the region.

The last rule \(\Theta_3\) in the proposed algorithm is the rule for terminating the computations. Since it is known in advance that, under modeling, only an approximate solution of the problem can be obtained, rule \(\Theta_3\) recommends terminating the computation at the moment when a “cycling” of the solution process has occurred, i.e., when for a comparatively long interval of time there is no decrease in the objective function.

In the proposed algorithm for modeling a stochastic program, at different stages different modifications of \(\Theta_1\) are used, basically coinciding with the algorithm described in [2]. The difference is that far from the boundary there is no need to test the point \(\vec{x}^{\,k}\) for stochastic feasibility: with high probability it is stochastically feasible if \(\vec{x}^{\,k}\in Y_d^{*}\). Therefore, far from the boundaries of the region \(Y_d\), a stochastically suitable direction is chosen by formulas (2, 3) from [2] on the basis of uniformly distributed points in a hypercube with center at \(\vec{x}^{\,k}\) and edge \(r\) (rule \(\Theta_1'\)). Near the boundaries, where the probability that the next point \(\vec{x}^{\,k+1}\) will become stochastically infeasible increases, the points scattered around \(\vec{x}^{\,k}\) are modeled in accordance with the distribution law of their coordinates \(x_i\), and, in accordance with formulation (2), the stochastic feasibility of the point \(\vec{x}^{\,k}\) is checked. At the same time, before computing the direction by (2, 3), the \(F_j\) are normalized.

with respect to the probability \(P_j\), in accordance with the formula of item 4 from [2] (rule \(\Theta_1''\)).

In the case where it turns out that \(\vec{x}^{\,k}\notin Y_{st}\), on the segment \(\vec{x}^{\,k-1}\vec{x}^{\,k}\) one finds the point nearest to \(\vec{x}^{\,k}\) that is stochastically admissible. This action is carried out by rule \(\Theta_2'\)—the “small” ricochet. The point found is denoted by \(\vec{x}^{\mathrm{opt},l}\) and is stored. After this, rule \(\Theta_2''\)—the “large” ricochet—operates, which gives, as the initial point for the next search stage, the center of gravity of all points generated at the preceding stage from \(\vec{x}^{0}\) to \(\vec{x}^{\mathrm{opt},l}\) (this point is formed gradually, as search points are generated).

Let us introduce the set \(J_k\), consisting of the indices of the violated (in the stochastic sense) constraints \(f_i(\vec{x})\). It may be assumed that repeated (\(L\)-fold) identification of the set \(J_k\) with one and the same set \(J_*\) means that the sought optimal point of the stochastic program is located near the boundaries whose indices constitute the current set \(J_k\). It is precisely in this case that rule \(\Theta_3\) operates, according to which the point \(\vec{x}^{\mathrm{opt},L}\) is regarded as a solution of the stochastic program if, in addition, the inequality

\[ \left|F(\vec{x}^{\mathrm{opt},L})-F(\vec{x}^{\mathrm{opt},L-1})\right|<\varepsilon, \]

is satisfied, where \(\varepsilon\) is a small quantity specified a priori.

3. Modeling algorithm

On the basis of what has been said above, the algorithm for modeling the stochastic program (2) can be described as follows (the parameters \(\delta, r, \varepsilon, q_1, q_2\) are specified a priori):

p. 1. We have \(\vec{x}^{0}\in Y_d^{*}\), \(J_*:=0\), \(l:=0\).

p. 2. Apply rule \(\Theta_1'\) to the point \(\vec{x}^{0}\); obtain the point \(\vec{x}^{1}\) (symbolically \(\Theta(\vec{x}^{0})=\vec{x}^{1}\).)

In general,

\[ \Xi(\vec{x}^{\,k-1})=\vec{x}^{\,k}. \]

p. 3. If \(\vec{x}\in Y_d^{*}\), then \(rq_1:=r\) and \(\to\) p. 2; otherwise \(\to\) p. 4.

p. 4. \(\delta/q_2:=\delta,\quad r/q_1:=r.\)

p. 5. \(\Theta''(\vec{x}^{\,k})=\vec{x}^{\,k+1}\). If \(\vec{x}^{\,k+1}\in Y_{st}\), then \(\to\) 2; otherwise \(\to\) p. 6.

p. 6. \(J_*:=J_k\cup J_*;\quad \Theta_2'(\vec{x}^{\,k+1})=\vec{x}^{\mathrm{opt},l}\) (small ricochet).

p. 7. If \(J_*=J_k\), then \(l:=l+1\) and \(\to\) p. 8; otherwise \(\to\) p. 9.

p. 8. If \(l=L\), then \(\to\) p. 10; otherwise \(\to\) p. 9.

p. 9. \(\Theta_2''(\vec{x}^{0},\vec{x}^{1},\ldots,\vec{x}_k,\vec{x}^{\mathrm{opt},l})=\vec{x}^{0}\) (large ricochet), \(\to\) p. 2.

p. 10. If rule \(\Theta_3\) is satisfied, then \(\to\) p. 11; otherwise \(\to\) p. 9.

p. 11. End.

Here, in the cycle over \(l\), one essentially solves \(L\) problems under the expectation sign in (2). From the solutions obtained, the best is selected, and ta-

thereby obtaining an approximate solution of problem (2). This modeling scheme should be especially economical if \(F(\vec{x})\) is nonlinear.

4. Computational experiment

To verify the effectiveness of the scheme presented, the problem in [1] was solved. The following table gives an idea of the efficiency of the algorithm. The first column gives the results for the stochastic case from [1] (more precisely, the solution by the simplex method). The subsequent columns give the optimal value \(F^*\) and the probabilities of satisfying the constraints on \(f_j\), obtained as a result of modeling after five ricochets from the boundary \(Y_{st}\). Here the notation is

\[ \begin{gathered} \mathrm{P}[1129{,}26 \leqslant f_1]=P_{11}, \quad \mathrm{P}[f_1 \leqslant 1256{,}26]=P_{12},\\ \mathrm{P}[501{,}15 \leqslant f_2]=P_{21}, \quad \mathrm{P}[f_2 \leqslant 623{,}15]=P_{22},\\ \mathrm{P}[245{,}51 \leqslant f_3]=P_{31}, \quad \mathrm{P}[f_3 \leqslant 298{,}58]=P_{32}. \end{gathered} \]

The solution was obtained for \(\varepsilon=0{,}01;\ q_1=q_2=2,\ \delta=15\).

Table 1

It is clear from the table that the proposed modeling method is quite workable. At the same time, rule \(\Theta_2\), according to which the ricochet from the boundary of the region \(Y_{st}\) is performed, cannot be regarded as entirely successful. This is indicated by the spread of \(P_{ij}\) in the table.

It should be specially noted that the experiment was carried out on a linear objective function, for which the equality
\(\mathrm{E}\{\sum_i a_i x_i\}=\sum_i a_i \bar{x}_i\) holds. Therefore the advantages due to solving problem (1) in the more general formulation (2), rather than (3), did not show themselves. In the case of a nonlinear form of the objective function \(F(\vec{x})\), these advantages may prove to be significant.

5. On the domain of application of the method

First of all, let us note that, according to the method set forth, stochastic linear programming problems with a stochastic matrix of coefficients \(\|a_{ij}\|\) of the system of linear constraints can in fact be solved without changing the program that has been compiled. It is only necessary to compile a program that models each of the coefficients according to the corresponding distribution law.

It is possible, however, to indicate a more general direction for applying the method, connected with a stagewise implementation of formulation (2) and with stagewise modeling of the solution in this formulation. It is important to emphasize that in this case optimization of \(F(\vec{x})\) of a nonlinear stochastic form is admissible in principle. This is all the more important because methods for solving problems of nonlinear stochastic programming do not exist [3], and modeling remains the only way of obtaining practical results.

Let us list possible applications of the modeling method for finding a stochastically feasible direction, as set forth in the present paper.

A) In problems where the model of the object is specified by means of statistical modeling, optimization must take place under conditions of interference owing to the limited number of statistical trials [5], [6].

An example of a problem of this kind is the problem of placing counters of charged particles, optimal in the sense of maximizing the number of particles reliably distinguishable from background radiation (the physical formulation of the problem is given in [7]).

B) The optimal-control problem considered in [8] is, in its general form, posed as a problem of nonlinear mathematical programming. In this problem, in addition to the requirement of maximizing the radiation intensity, constraints are added on the particle energy and the time of their dumping. Taking into account the instability of the position of the extremum near the boundary of the nonzero radiation region and the presence of interference, we shall evidently have a stochastic programming problem.

C) There is a broad class of problems in which the analytic form of the function to be optimized is specified by an operator of such a complex kind (integro-differential, multidimensional integral, minimax operator) that the computation of a single value of the function with high accuracy takes an excessively large amount of machine time. In order to reduce the cost of computing functions of such a complex kind, one has to reduce the accuracy of their computation, which obviously leads to an increase in the level of interference in the process of searching for the extremum. Operators of this kind arise in large numbers in problems of optimizing structural designs [9].

In all these cases, application of the approach developed in the present paper, combining modeling of the problem in its initial formulation with stagewise integral smoothing in the search process, makes it possible to hope for obtaining a practically applicable solution.

REFERENCES

1. V. V. Zakharov, V. P. Tarasenko. Optimal control of an object in a static regime as a problem of stochastic programming. Proceedings of the 1st All-Union Symposium on Static Problems in Engineering Cybernetics (Nonlinear and Optimal Systems), Moscow, “Nauka,” 1971.

2. V. V. Zakharov. Identification of a static object by means of integral projections. Izvestiya TPI. Present collection.

3. J. Hadley. Nonlinear and Dynamic Programming. Moscow, “Mir,” 1967.

4. G. Zoutendijk. Methods of Feasible Directions. Moscow, IL, 1963.

5. N. P. Buslenko. Mathematical Modeling of Production Processes. Moscow, “Nauka,” 1964.

1. G. I. Marchuk, S. M. Ermakov. The Monte Carlo Method and Methods of Computational Mathematics. Moscow, Atomizdat, 1967.

2. V. V. Kuznetsov et al. Measurement of the asymmetry of photoproduction of \(\pi^+\)-mesons on hydrogen by linearly polarized \(\gamma\)-quanta. Yadernaya fizika, 1971, No. 13.

3. V. V. Zakharov et al. Algorithms for controlling a synchrotron in the “computer–adviser” mode. Izvestiya TPI. The present collection.

4. A. S. Volmir. Flexible Plates and Shells. Moscow, Tekhteorizdat, 1956.