APPLICATION OF RANDOM SEARCH TO A PROBLEM IN NONLINEAR PROGRAMMING

V. V. Zakharov, V. P. Ivanchenkov, V. A. Filinova

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

For solving nonlinear programming problems of small dimension and with a small number of constraints, it seems natural to employ methods of statistical search, which provide simple machine algorithms. In the present article it is proposed to apply, for this purpose, the random search method [1] within the framework of one scheme of penalty functions [2].

1. Problem statement

Let \(F(x)\) and \(f_j(x)\) \((j \in I,\ I=\{1,2,\ldots,m\})\) be convex functions of the vector argument \(x=(x_1,\ldots,x_n)\). The paper considers the problem of finding

\[ \min_X \{F(X)\mid f_j(X)\geq 0,\ j\in J\}. \tag{1} \]

It should be emphasized that the convergence of the considered variant of the penalty-function method was proved in [3] for convex functions \(F, f_j\). Nevertheless, in that same work examples were considered of solving nonconvex problems of “good behavior,” i.e., such that the curvature of the surfaces \(F(x)=\mathrm{const}\), \(f_j(x)=\mathrm{const}\) was not too large. The penalty-function method was chosen in the present work in combination with random search motion because it made it possible to hope for an approximate solution of the practical problem considered below, with respect to whose convexity there was no certainty.

2. Description of the method

The solution of the problem in [2], [3] is reduced to the stage-by-stage minimization of the functions

\[ A_m(X)=\sum_{j=1}^{m'} f_j(X), \tag{2} \]

which represent a penalty for approaching the boundary of the domain of definition \(R\) of the function \(F(X)\).

In the present work, at the \((k+1)\)-st stage of the search for \(\min_{X\in R} A_m(X)\) (the domain \(R=\{f_j(X)\geq 0,\ j\in J\}\)), function (2) was supplemented by the constraint:

\[ f_{m+1}(X)=F(X^{0k})-F(X)\geq 0, \]

where \(X^{0k}\) is the point of the search trajectory at the \(k\)-th stage at which the function \(A_{m+1}(X)=A(X)\) attains its least value.

The modification of the penalty method used consisted also in the fact that at the \(k\)-th stage the problem of finding \(\min A(X)\) with high accuracy was not posed. Moreover, exit from the \(k\)-th minimization stage occurred when the resolving power of the random-search method (the method step) had decreased to some fixed value \(\tau\).

3. Algorithm

Let us consider the algorithmic implementation of the method at a certain \(k\)-th stage. The main task of the stage is to obtain a point \(X^{0,k+1}\) such that \(F(X^{0,k+1})<F(X^{0k})\). The descent to this point is carried out over the penalty surface \(A(X)\) by the random-search method with step \(\tau\). It is not difficult to see that \(f_{m+1}(X^{0k})=0\) and \(A(X^{0k})=+\infty\). Therefore it makes sense to take the inequality for the additional \((m+1)\)-st constraint in the form

\[ f_{m+1}(X)=F(X^{0k})-F(X)+\varphi\geq 0, \]

where \(\varphi>0\).

Now one can already seek \(\min A(X)\) by the random-search method.

In accordance with the known random-search algorithm with re-crossing [1], we take the point \(X^k\) as the initial point of the search trajectory at the \(k\)-th stage: \(X^k=X_0^k\). In a cubic neighborhood of the point \(X_0^k\) (the cube edge is equal to \(\tau\)), we perform a uniform scattering of points until at least one of them is found at which the value of \(A(X)\) is less than \(A(X_0^k)\). We denote this successful point by \(X_1^k\), again take it as the center of a cube with edge \(\tau\), and again carry out uniformly distributed trials in this cube, and so on.

For the parameter \(\tau\) we introduce adaptation to the search conditions: the value \(\tau\) will be decreased by a factor of \(q\) if at least one of the \((m+1)\) constraints is violated; the value \(\tau\) will be increased by a factor of \(q\) in the opposite case. The initial value was assumed equal to \(\bar{\tau}\). If \(\tau<\tau\) was satisfied, then the \(k\)-th stage was considered completed and \(X_L^k=X^{k+1}\) was set; moreover, obviously,

\[ A(X^k=X_0^k)>A(X_1^k)>\cdots>A(X_{L-1}^k)>A(X_L^k=X^{k+1}). \]

It may be assumed (and practical computations confirm this) that at the final stages the search will proceed in a narrow ravine. As a result, in the cube with center \(X_i^k\) it is necessary to carry out an extremely large number of trials before \(A(X_i^k)>A(X_{i+1}^k)\) is satisfied. To counter such a situation, the following modification of the computational algorithm was introduced. If the number of trials in the neighborhood of \(X_i^k\) exceeded \(N_z\) (cycling), then from the point \(X_0^k\) a secondary descent is carried out according to the described algorithm until it is not satisfied

\(\tau<\bar{\tau}\) (i.e., the natural termination of the stage) or cycling will again occur at the point \(X_{l_2}^k\). Let \(A(X_{l_2}^k)<A(X_{l_1}^k)\). In this case, as \(X_L^k\) the point chosen is

\[ X_L^k=\frac{X_{l_2}^k-X_R^k}{2}, \]

where the point \(X_R^k\) is the point of intersection of the boundary of the domain \(R\) with the straight line \(X_{l_2}^k X_{l_1}^k\) in the direction of decrease of \(A(X)\) along this straight line. At this point the \(k\)-th stage also ends.

The entire procedure for solving problem (1) continues until the prescribed number of evaluations of the function \(F(X)\), equal to \(N\), has been exhausted.

4. Flowchart and Its Description

The structure of the algorithm can be clarified in more detail from the flowchart presented. Let us explain it.

Block 1 carries out the formation of the program commands, taking into account the number of variables \(n\) and the number of constraints \(m\).

Block 2. The indicator \(d\) determines the number of descents from the initial point \(X_L^k\); if \(d=1\), then from the point \(X_L^k\) the descent has already been performed twice. The working step \(\tau\) at the beginning of the \(k\)-th search stage is assigned the value equal to \(\bar{\tau}\). In block 2 the initial values of \(F(X)\) and \(A(X)\) are computed.

Block 3. The point \(X^k\) is the search center at the \(k\)-th stage. The initial value \(X^k\) is equal to \(X_0^k\).

Block 4. The counter \((\tau)\) is needed for counting the number of consecutive points at which \(A(X^\xi)>A(X_l^k)\), where \(X^\xi\) is a random point uniformly distributed in the cube with center \(X^k\).

Blocks 5–6 perform the computation of \(f_j(X^\xi)\). If the condition \(f_j(X^\xi)\ge 0\) is not satisfied, then the step \(\tau\) is reduced by a factor of \(q\) (block 17).

Block 7. The counter \((F)\) controls the number of evaluations of \(F(X)\); if \(\operatorname{ct}(F)=N\), then the found \(F_{\min}\) and \(X_{\min}\) are printed. The search terminates (block 18).

Blocks 8–9. \(M\) is the current minimum value of \(F(X)\), and \(X^M\) is the point corresponding to \(M\).

Block 10. The working step \(\tau\) is increased by a factor of \(q\) if \(m\) constraints have been satisfied for the point \(X\).

Block 11 checks: \(A(X^\xi)<A(X_l^k)\). If yes, then \(\to\) block 26; \(A(X_{l+1}^k):=A(X^\xi)\), \(X_{l+}^k:=X^\xi\), otherwise \(\operatorname{ct}(\tau):=\operatorname{ct}(\tau)+1\).

Blocks 12–16. If \(\operatorname{ct}(\tau)=N_\xi\), then the step \(\tau\) is reduced by a factor of \(q\); otherwise (in block 14) it is checked whether \(d=0\).

If yes, then \(X_{l_1}^k:=X^k\) and \(d=1\).

With step \(\tau=\tau/2\), a descent is once again carried out until the point \(X_{l_2}^k\) is obtained.

Blocks 19–22. If \(d\ne 0\), then \(X_{l_2}^k:=X_l^k\).

The check \(A(X_{l_1}^k)<A(X_{l_2}^k)\) is performed in order to choose the direction of motion along the straight line \(X_{l_1}^k X_{l_2}^k\).

Algorithm flowchart

1. Formation of commands,

\[ \operatorname{cnt}(F)=0 \]

1. \[ d:=0,\quad \tau:=\bar{\tau} \]

Computation of \(F(X_L^k)\) and \(f_j(X_L^k)\).

1. \[ A(X_l^k):=A(X_L^k),\quad X^k:=X_L^k \]

2. \[ \operatorname{cnt}(\tau):=0 \]

3. \[ X^\xi:=X_l^k+\tau(2\xi-1) \]

4. Computation of \(f_j(X^\xi)\).

\[ X^\xi\in R? \]

no / yes

1. Computation of \(F(X^\xi)\),

\[ \operatorname{cnt}(F):=\operatorname{cnt}(F)+1,\quad \operatorname{cnt}(F)=N? \]

yes / no

1. \[ F(X^\xi)<M? \]

2. \[ M:=F(X^\xi),\quad X^M:=X^\xi \]

3. \[ \tau:=\tau\cdot q \]

4. \[ A(X^\xi)<A(X_l^k)? \]

   yes / no

5. \[ \operatorname{cnt}(\tau):=\operatorname{cnt}(\tau)+1 \]

   \[ \operatorname{cnt}(\tau)=N_\vartheta? \]

   yes / no

6. \[ \tau:=\tau/q,\quad \tau<\underline{\tau}? \]

   yes / no

7. \[ d=0? \]

   yes / no

8. \[ d:=1,\quad X_{l_1}^k:=X_l^k \]

9. \[ \tau:=\bar{\tau}/2 \]

10. \[ \tau:=\tau/q,\quad \tau<\underline{\tau}? \]

    yes / no

11. Print \(F_{\min}\) and \(X_{\min}\). Stop.

12. \[ X_{l_2}^k:=X_l^k \]

13. \[ A(X_{l_1})<A(X_{l_2})? \]

    yes / no

14. \[ X_{l_1}^k\text{—beginning of motion} \]

15. \[ X_{l_2}^k\text{—beginning of motion} \]

16. Motion along the straight line \(X_{l_1}^k X_{l_2}^k\) with step \(\tau/Q\).

17. \[ X_{l_1}^k,\ X_{l_2}^k\in R? \]

    yes / no

18. \[ X_L^k=\frac{X_{l_1}^k+X_R^k}{2} \]

19. \[ X_{l+1}^k:=X^\xi,\quad A(X_{l+1}^k):=A(X^\xi) \]

Algorithm flowchart

Blocks 23–26 carry out motion along the straight line \(X_{l_1}^k X_{l_2}^k\) with step \(\tau/Q\) up to the boundary of the region \(R\), after which the initial point of the next stage of the search for \(\min A(X)\) is determined:

\[ X^{k+1}=X_L^k=\frac{X_{l_1}^k+X_R^k}{2}. \]

5. Description of the Program

The program was compiled for the BESM-4 computer.

a) Parameters used in the method.

Table

The programmer may change these parameters at his discretion.

b) Memory allocation and calling the method program.

The method program is located in cells 0030 through 0330. Cells 0400 through 1277 are occupied by the working field for intermediate and final results. The program uses the constants IS-2.

The program is called by the commands:

where

\(n \le 63\) is the number of variables,

\(m\) is the number of constraints on the optimized parameters,

\(F\) is the cell at the beginning of the count of \(f_j(X)\) and \(F(X)\),

\(n\), \(m\), \(F\) are specified in octal form.

The programs for computing \(f_j(X)\) and \(F(X)\) may be placed anywhere in memory, beginning with cell 1300. Cell \(F-1\) is intended for return to the method program and must be free.

The components of the vector \(X\) occupy cells \(0601\)—\(0600+n\). Here the programmer also specifies the initial point \(X^0\).

It is recommended first to compute \(f_j(X)\) and enter these values in cells \(2701\)—\(2700+m\), and then to check \(f_j(X)\ge 0,\ j\in J\). If at least one of the constraints is not satisfied, it is necessary to transfer control

in cell 0076; otherwise compute \(F(X)\). The value \(F(X)\) must be written to cell 0600. Upon completion of the calculation of \(F(X)\), transfer control to cell \(F-1\).

The current minimal value of the function and the point corresponding to it are stored in cells \(0400\)–\(0400+n\). The working step \(t\) is stored in cell 0060.

6. Example of applying the method to a production problem

The problem consisted in determining the parameters of the charged-particle transport channel of an accelerator that are optimal from the point of view of the minimum cost of the structure as a whole. The channel consisted of 13 sections: 7 quadrupole lenses and 6 particle-drift sections. Each section is uniquely determined by a set of parameters in [4]

\[ k_i,\ l_i,\ \beta_i,\ \alpha_i,\ \gamma_i \quad (i=\overline{1,13}). \]

The quantities \(\beta_i,\ \alpha_i,\ \gamma_i\) are computed by the following formulas.
For the horizontal plane of lens focusing:

\[ \beta_i=\gamma_{0i}\sin^2\Theta_i -2\alpha_{0i}\frac{1}{\sqrt{k_i}}\sin\Theta_i\cos\Theta_i +\beta_{0i}\cos^2\Theta_i, \]

\[ \alpha_i=\gamma_{0i}\frac{1}{\sqrt{k_i}}\cos\Theta_i\sin\Theta_i +\alpha_{0i}(\cos^2\Theta_i-\sin^2\Theta_i) +\beta_{0i}\sqrt{k_i}\sin\Theta_i\cos\Theta_i, \]

\[ \gamma_i=\gamma_{0i}\cos^2\Theta_i +2\alpha_{0i}\sqrt{k_i}\cos\Theta_i\sin\Theta_i +\beta_{0i}k_i\cos^2\Theta_i,\quad (i=1,\ldots,7). \]

For the particle-beam drift sections:

\[ \left. \begin{aligned} \beta_i&=\beta_{0i}-2\alpha_{0i}l_i+\gamma_{0i}l_i^2,\\ \alpha_i&=\alpha_{0i}-\gamma_{0i}l_i,\\ \gamma_i&=\gamma_{0i} \end{aligned} \right\} \quad k=0.\quad (i/2=1,\ldots,6). \]

For the vertical plane of lens focusing:

\[ \beta_i=\gamma_{0i}\frac{1}{k_i}\operatorname{ch}^2\Theta_i -2\alpha_{0i}\frac{1}{\sqrt{k_i}}\operatorname{ch}\Theta_i\operatorname{ch}\Theta_i +\beta_{0i}\operatorname{ch}^2\Theta_i, \]

\[ \alpha_i=-\gamma_{0i}\frac{1}{\sqrt{k_i}}\operatorname{ch}\Theta_i\operatorname{ch}\Theta_i +\alpha_{0i}(\operatorname{ch}^2\Theta_i+\operatorname{ch}^2\Theta_i) -\beta_{0i}\sqrt{k_i}\operatorname{ch}\Theta_i\operatorname{ch}\Theta_i, \]

\[ \gamma_i=\gamma_{0i}\operatorname{ch}^2\Theta_i -2\alpha_{0i}\sqrt{k_i}\operatorname{ch}\Theta_i\operatorname{ch}\Theta_i +\beta_{0i}k_i\operatorname{ch}^2\Theta_i,\quad (i=1,\ldots,7). \]

Here, in the formulas, everywhere \(\Theta_i=l_i\sqrt{k_i}\).
In this notation the objective function \(F(X)\) has the form:

\[ F(X)=4\sum_{i=1}^{7}l_i\sqrt{C\cdot \beta^*_{\mathrm{hor}.i}\cdot \beta^*_{\mathrm{vert}.i}},\quad C=\mathrm{const}. \]

Here \(i\) is the number of the channel section with a lens, and \(\beta^*_{\mathrm{hor}.i},\ \beta^*_{\mathrm{vert}.i}\) are the maximum values on section \(i\) in the vertical or horizontal plane of lens focusing. Let us relabel the vector of variables as

\[ X=(x_1=x_3,\ x_2=k_4,\ x_3=l_4,\ x_4=l_6,\ x_5=l_8,\ x_6=l_{10}). \]

The variables \(x\) were subject to the constraints:

\[ 0.09 \le x_1 \le 0.64;\quad 0.09 \le x_2 \le 1;\quad 4 \le x_3 \le 10; \]

\[ 4 \le x_4 \le 10;\quad 4 \le x_5 \le 9;\quad 4 \le x_6 \le 9; \]

\[ \sum_{i=3}^{6} X_i = 31.92. \]

The remaining parameters were not changed.
For the existing tract the values of the variables were equal to:

\[ x_1=x_2=0.438;\quad x_3=8;\quad x_4=7.96; \]

\[ x_5=7.96;\quad x_6=8;\quad F(X)=0.08230. \]

After minimizing \(F(X)\) by the method described above, we obtained:

\[ F(X)=0.00135; \]

\[ x_1=0.184;\quad x_2=0.375;\quad x_3=5.084; \]

\[ x_4=8.460;\quad x_5=8.895;\quad x_6=9.480. \]

The machine time spent searching for \(F(X^*)=0.00135\) from the given initial point was approximately 13 min. A suitable result by the random-search method without any modifications had previously been obtained in 4 hours of BESM-4 operation.

The algorithm of the method in the form in which it is presented in the flowchart is rather sensitive to ravines; therefore, when to the existing constraints there were added the constraints \(\beta_7(X)\le 0.0833;\ \gamma_7(X)\le 33.3\), which had a pronounced ravine form, the method quickly cycled.

Thus, in order to increase the effectiveness of the method, its resolving power in a ravine situation should be increased.

Appendix

Program of the Method

0100        56     —     0303   0056              1   1   32   0001   0130   7777
   1     4  55     —     7734   0072              2       00     —      —     0062
   2     1  14   0064   0072   0072              3       72     —     0072     —
   3     4  55     —     7732   0073              4       54   0107   0040   0053
   4     4  55     —     7731   0074              5       07   0053   0040   0053
   5        14   0114   0074   0074              6       13   0041   0053   0053
   6        00   0047     —     0400              7       23   0053     —     0040
   7        00     —       —     0067
0110        00     —       —     0061       0140       05   7762   0040   0054
   1        72     —     0072     —                1       02   0054   7761   0054
   2     5  00   0600     —     0500              2       05   0054   0060   0054
   3     1  32   0002   0112   7777              3   3   01   0054   0700   0600
   4        00     —       —     0063              4   1   32   0002   0134   7777
   5        00   0043     —     0060              5       72     —     0073     —
   6        72     —     0073     —                6   3   16   0147     —     7777
   7     3  16   0120     —     7777              7       01   7761   0067   0067
0120        56   0600   0300   0500       0150       15   0067   0046     —
   1     2  04   7761   2700   0075              1       36     —     0272     —
   2        01   0075   0065   0065              2       02   0600   0400     —
                                                   3       76     —     0157   0070

REFERENCES

1. L. A. Rastrigin. Statistical Search Methods. Moscow, “Nauka,” 1968.

2. C. W. Carroll. The CRST for Optimizing Nonlinear Restrained systems. Oprns. Res., v. 9, No. 2, 1961.

3. A. V. Fiacco, G. P. McCormick. Computational Algorithm for SUMT for Nonlinear programming. Manag. Sc., v. 10, No. 4, 1964.

4. V. P. Ivanchenkov. Questions of the application of correlation-extremal systems in problems of electron optics. Dissertation. Tomsk, 1969.