Multidimensional Constrained Global Optimization in Domains with Computable Boundaries

Abstract

Multidimensional constrained global optimization problem with objective function under Lipschitz condition and constraints generating a feasible domain with computable boundaries is considered. For solving this problem the dimensionality reduction approach on the base of the nested optimization scheme is used. This scheme reduces initial multidimensional problem to a family of one-dimensional subproblems and allows applying univariate methods for the execution of multidimensional optimization. Sequential and parallel modi cations of well-known information-statistical methods of Lipschitz optimization are proposed for solving the univariate subproblems arising inside the nested scheme in the case of domains with computable boundaries. A comparison with classical penalty function method being traditional means of taking into account the constraints is carried out. The results of experiments demonstrate a signi cant advantage of the methods proposed over the penalty function method.

Keywords

• global optimum
• multidimensional problems cursive optimization
• computable boundaries

-

The global optimization problem [1{3] is considered in the following form: y 2 Q

RN ; D = y 2 RN : ai yi bi; 1 i

N ; Q = y 2 D : gi(y) 0; 1 i m ; is de ned by constant ( 2 ) and functional ( 3 ) constraints. The objective function is supposed to satisfy in the domain Q the Lipschitz condition ( 1 ) ( 2 ) ( 3 ) with the Lipschitz constant L > 0, where the denotation k k signi es the Euclidean norm in RN .

The traditional way of solving the constrained problems of mathematical programming consists in solving the unconstrained problem where

F (y) ! min;

y 2 D;

F (y) = f (y) + CG(y): An auxiliary function G(y) (penalty function) satis es the condition G(y) = 0; y 2 Q;

G(y) > 0; y 2= Q; and the constant C > 0 (see [1, 2, 4]).

If the constant C is su ciently large and functions f (y) and G(y) are, for instance, continuous, solutions of the problems ( 1 ) and ( 5 ) coincide. For some function classes there is a nite penalty constant providing the coincidence of solutions (Eremin-Zangwill exact penalty functions [4, 5]).

As an example of penalty function one can take the function

G(y) = max 0; g1(y); : : : ; gm(y) :

If all the functions gj (y); 1 j m are continuous in D, the function ( 8 ) is continuous as well.

On the one hand, the problem ( 5 ) is simpler than the initial problem ( 1 ) because of simplifying the feasible domain D. On the other hand, a choice of the penalty constant is rather di cult. If it is insu ciently large the global minimizer of the problem ( 5 ) can fall out of feasible domain. If the penalty constant is too large, it worsens the properties of the function F (y) in comparison with the initial function f (y) (the function F (y) can have a ravine surface, the Lipschitz constant of F (y) can increase considerably, etc.).

Another factor which in uences the complexity of the optimization significantly is the dimension N of the problem. To overcome this complexity, the approaches connected with reducing the multidimensional problem to one or several univariate subproblems are often applied. We consider one of approaches to the dimensionality reduction based on the nested optimization scheme which replaces solving the multidimensional problem ( 1 ) by solving a family of univariate subproblems connected recursively. In the framework of this approach in combination with di erent univariate global search methods [6{12] many sequential and parallel multidimensional algorithms have been proposed [7, 8, 13{16] and applied to practical problems (see, for example, papers [17{19]). The other interesting approaches to parallelizing global search algorithms can be found in publications [21{24].

The scheme of nested optimization consists in the following [3, 7, 20]. Let us introduce the notations ui = (y1; : : : ; yi); i = (yi+1; : : : ; yN ); ( 9 ) ( 5 ) ( 6 ) ( 7 ) ( 8 ) allowing to write down the vector y as a pair y = (ui; i) for 1 i N 1. Assume that y = 0 if i = 0 and y = uN for i = N .

Let us de ne also two series of sets. The rst series contains sections of the domain Q:

S1 = Q;

Si+1(ui) =

i 2 RN i : (ui; i) 2 Q ; The second collection of sets consists of projections

Pi+1(ui) = yi+1 2 R : 9(yi+1; i+1) 2 Si+1(ui) ; i i

N N 1: 1; of the sections Si+1(ui) onto the axis yi+1.

Then, according to [3] and [8] the basic relation of the nested optimization scheme min f (y) = min min y2Q y12P1 y22P2(u1) : : :

min yN 2PN (uN 1) f (y) takes place.

Now let us introduce one more family of functions generated by the objective function f (y): f N (y)

f (y); f i(ui) = min f i+1(ui; yi+1) : yi+1 2 Pi+1(ui) ; i

N 1; de ning in the projections

Qi =

i ui 2 R : 9(ui; i) 2 Q; ; 1

N; of the domain Q onto the coordinate axes y1; : : : ; yi.

As it follows from ( 11 ), in order to solve the problem ( 1 ) { ( 3 ) it is su cient to solve a one-dimensional problem f 1(y1) ! min; y1 2 P1

R1: f i(ui 1; yi) ! min; yi 2 Pi(ui 1); where the xed vector ui 1 2 Qi 1.

According to ( 14 ), each estimation of the function f 1(y1) at some xed point y1 2 P1 consists in solving a one-dimensional problem f 2(y1; y2) ! min; y2 2 P2(y1)

R1: This problem is a one-dimensional minimization with respect to y2 since y1 is xed, etc., up to solving the univariate problem with xed vector uN 1. f N (uN 1; yN ) = f (uN 1; yN ) ! min; yN 2 PN (uN 1);

On the whole, the solving of the problem ( 1 ) { ( 3 ) is reduced to solving a family of \nested" one-dimensional subproblems 1 0 1 i ( 10 ) ( 11 ) ( 12 ) ( 13 ) ( 14 ) ( 15 ) ( 16 ) ( 17 ) ( 18 ) ( 19 )

The approach of nested optimization can be applied to the problem ( 5 ) as well. In this case the domains of one-dimensional search in ( 19 ) are simple, namely, they are closed intervals Pi = [ai; bi]; 1 i N .

However, in general case subject to continuity of the penalty function G(y) the projection Pi(ui 1) is a system of non-intersecting closed intervals Pi(ui 1) = qi [ [ais; bis]: s=1 where the number of the intervals qi and their bounds ais; bis; 1 on the vector ui 1, i.e., s qi, depend qi = qi(ui 1); ais = ais(ui 1); bis = bis(ui 1):

If the domain Q is such that it is possible to obtain the explicit (analytical) expressions for the values qi; ais; bis as functions of ui 1 2 Qi 1 for all 1 i N , then the feasible set Q is called the domain with the computable boundaries. 2

Sequential and parallel algorithms of global search over a system of closed intervals Let us introduce a uni ed form for one-dimensional problems ( 19 ) arising inside the nested scheme '(x) ! min; x 2 X = q [ [ s; s]; s=1 ( 20 ) ( 21 ) ( 22 ) where q 1; s s; 1 s q; s < s+1; 1 s q 1.

One of the most e cient methods of univariate multiextremal optimization is the information-statistical algorithm of global search (AGS) proposed by Strongin [7] for the case ( 22 ) with q = 1. We will consider a sequential modi cation of AGS for the more general situation ( 22 ) when q > 1, i.e., the domain X of one-dimensional search consists of several intervals and call this modi cation as AGS-G. This method belongs to the class of characteristical algorithms [21] and can be described like these methods in the following manner.

Let the term \trial" denote an estimation of the objective function value at a point of the domain X, k be a trial number, xk be the coordinate and zk = '(xk) be the result of k-th trial.

The initial stage of AGS-G consists in carrying out 2q initial trials at points x1 = 1; x2 = 1; x3 = 2; x4 = 2; : : : ; x2q 1 = q; x2q = q, i.e., x2j 1 = j ; x2j = j ; 1 j q, with corresponding trial results zj = '(xj ); 1 j 2q.

The choice of a point xk+1; k > 2q, for implementation of a new (k + 1)-th trial consists in the following steps.

Step 1 Renumber by subscripts the points x1; : : : ; xk of previous trials in increasing order, i.e., 1 = x1 < x1 < < xk = q; ( 23 ) and juxtapose to them the values zj = '(xj ); 1 j k, obtained earlier, i.e., zj = '(xl), if xj = xl.

Step 2 Assign to each interval (xj 1; xj ); 2 j k, generated by the points from ( 25 ) a feasibility indicator j in the following way: if xj 1 = j and xj = i+1, or xj 1 = xj , then j = 0, otherwise, j = 1.

Step 3 For feasible intervals (xj 1; xj ) with indicator j = 1 calculate the divided di erences and the value Step 4 Determine an adaptive estimation j = jzj xj zj 1j ; xj 1 k = maxf j : 2 j

k; j = 1g: Mk = (r k; 1; k > 0; k = 0; xk+1 = xt + xt 1 2 zt

zt 1 ; 2Mk 2(zj + zj 1); (27) ( 24 ) ( 25 ) (26) (28) (29) where r > 1 { a parameter of the method.

Step 5 Juxtapose a numerical value (characteristic)

R(j) = Mk(xj xj 1) +

(zj Mk(xj zj 1)2 xj 1) to each feasible interval (xj 1; xj ); 2 j k; j = 1.

Step 6 Among feasible intervals select the interval (xt 1; xt) which the maximal characteristic R(t) corresponds to.

Step 7 Choose in the interval (xt 1; xt) the point

as the coordinate of (k + 1)-th trial and calculate the value zk+1 = '(xk+1). Step 8 Increase the iteration number k by 1 (k = k + 1) and go to Step 1.

The computational scheme described above generates an in nite sequence of trials. It can be truncated by introducing a termination criterion. For characteristical algorithms as such criterion the inequality xt xt 1 "; is used, as a rule, where " > 0 is a prede ned search accuracy, i.e., the search is considered to have been completed when the length of the interval with maximal characteristic is less than accuracy ".

As AGS-G is a simple modi cation of AGS it is easily to show that their convergence conditions (see [3]) are the same. In particular, convergence to global minima is provided by ful llment of the inequality Mk > 2L, where L is the Lipschitz constant of the function ( 22 ).

In order to describe a parallel version of AGS-G let us assume that at our disposal there are p > 1 independent processors, and modify the steps 6 { 8 of the algorithmic scheme presented above. Initial stage of executing trials in end points of intervals ( 22 ) can be realized either sequentially or in parallel { no matter.

Step 6 Arrange all characteristics (27) in the decreasing order

R(t1)

R(t2) : : : (30) and take intervals (xtj 1 ; xtj ); 1 j p, which the rst p characteristics in the series (30) correspond to.

Step 7 Within intervals with maximal characteristics determined at the previous step calculate points In order to compare the e ciency of the penalty function method and the method of explicit boundaries in the framework of nested optimization scheme a class of known two-dimensional multiextremal test functions has been considered and on the base of this functions taken as objective ones the problems of constrained optimization have been constructed with constraints forming a feasible domain with computable boundaries, which provide the complicated structure (non-convex and disconnected) of this domain.

For the experiment test functions two-dimensional [10, 25] as objective ones in the problem ( 1 ) have been choosen. The domain ( 2 ) is the square 0 y1; y2 1 and 5 non-linear constraints g1(y1; y2) = 0:5 exp( y1)

y2 g2(y1; y2) = g3(y1; y2) = 4(y1 4(y2 0:9)2 + y2 0:6)2 + y1 0:25 0:8 0:7 0; 0; 0; g4(y1; y2) = g5(y1; y2) = 10jy2 (y1 0:5y1 0:2)2 0:1j + j sin(7 y1)j (y2 0:8)2 + 0:1 0 0; form the feasible domain ( 3 ).

Hundred test problems of this type have been minimized by means of the nested optimization scheme with the algorithm AGS-G applied inside for univariate optimization ( 19 ) (method MAGS-G). Moreover, the same problems have been solved by the penalty function method (MAGS-P) with di erent penalty constants C. For illustration of methods behavior Figure 1 shows level curves of a test function, the feasible domain and trial distribution in the region ( 2 ) for MAGS-G and MAGS-P (trials are marked with crosses). Figure 1 demonstrates that all the trials of MAGS-G are placed in the feasible domain only while MAGS-P carried out a signi cant part of trials of the the feasible region.

For conclusive comparison of the optimization algorithms considered above the method of operating characteristics [3, 16] was used. In the framework of this method several test problems are taken and each of them is solved by an optimization algorithm with a given set of its parameters. After the experiment one can determine a number P of problems solved successfully and average number K of trials spent by the algorithm. Repeating such experiment for di erent sets of parameters we obtain a set of pairs (K; P ) called operating characteristic of the algorithm. It is convenient to represent operating characteristics of the algorithms compared on the plane with axes K and P . Such representation enables to compare e ciencies of the algorithms in a visual manner. If for the same K the operating characteristic of a method is located above the characteristic of another one, the rst method is better as it has solved successfully more problems than its rival spent the same computational resources.

Figure 2 contains the operating characteristics of MAGS-G and MAGS-P with penalty constants C = 1, C = 10, C = 100, obtained after minimization of 100 two-dimensional random functions [10, 25] with r = 3 and di erent values of accuracy " from (29).

The experiment con rms known situation when for small penalty constant MAGS-P are not able to nd the global minimum of certain constrained problems. On the other hand, if the penalty constant is overestimated MAGS-P spends too many trials for searching the global minimum. Moreover, the MAGSP e ciency worsens because of placing trials out of the feasible domain while MAGS-G executes trials at feasible points only.

For evaluation of e ciency of parallelizing the nested optimization scheme let us take again the test class (33) with the same 5 constraints and apply to optimization the nested scheme where for solving problem ( 16 ) the sequential method AGS-G is used, but the parallel algorithm PAGS-G with di erent numbers of processors performs solving the univariate problems for the coordinate y2. This parallel multidimensional method will be called PMAGS-G. Let K( 1 ) be the average number of trials executed by MAGS-G and K(p) be the average number of trials spent by PMAGS-G with p > 0 processors for optimizing 100 test functions. As a criterion of e ciency the speed-up in trials is taken. This criterion is de ned [21] as s(p) = K( 1 )=(pK(p)) and characterizes the speed-up under assumption that the time spent for realization of algorithmic scheme and data transmission is much less than the time of objective function evaluations. The graph of speed-up in trials is presented in Figure 3.

Experiments were performed using two-socket machine with Intel R XeonTM E5-2680 v3 processors (Haswell, 2.5GHz, 12 Cores; OS: Red Hat Enterprise Linux Server 6.6; compiler: GCC 4.8.2). Library Intel R Threading Building Blocks 4.3 Update 3 was used for implementing parallel version of the method. 4

Conclusion A class of multidimensional multiextremal problems with special type of constraints has been considered. Constraints of this type can generate non-convex and even disconnected feasible domains. For solving the problems under consideration sequential and parallel global optimization methods on the base of nested optimization scheme have been adapted. Decision rules of these methods provide estimating the objective function within the feasible domain only and do not require any parameters for taking into account the constraints as opposed to classical penalty function method. Computational experiments demonstrate advantages of the algorithms proposed in comparison with the penalty function method for constraints of the considered type.

As a further way of the development of researches it would be interesting to investigate di erent parallelization schemes inside the nested optimization approach in combination with other methods of univariate optimization. Acknowledgments This research was supported by the Russian Science Foundation, project No. 15-11-30022 \Global optimization, supercomputing computations, and applications".

References

1. Horst , R. , Pardalos , P. , Thoai , N.: Introduction to Global Optimization. Kluwer Academic Publishers, Dordrecht ( 1995 )

2. Pintr , J.: Global Optimization in Action (Continuous and Lipschitz Optimization: Algorithms, Implementations and Applications). Kluwer Academic Publishers, Dordrecht ( 1996 )

3. Strongin , R. , Sergeyev, Ya.: Global Optimization with Non-convex Constraints: Sequential and Parallel Algorithms . Kluwer Academic Publishers, Dordrecht ( 2000 )

4. Eremin , I. , Astafyev , N.: Introduction to the theory of linear and convex programming . Nauka , Moscow ( 1976 )

5. Di Pillo , G. , Grippo , L. : Exact penalty functions in constrained optimization . SIAM Journal on Control and Optimization 27 ( 6 ), 1333 { 1360 ( 1989 )

6. Piyavskij , S.: An Algorithm for Finding the Absolute Extremum of a Function . USSR Comput. Math. Math. Phys . 12 ( 4 ), 57 { 67 ( 1972 )

7. Strongin , R.: Numerical Methods in Multiextremal Problems (Information - Statistical Algorithms) . Nauka, Moscow ( 1978 )

8. Grishagin , V. , Strongin , R.: Optimization of multiextremal functions subject to monotonically unimodal constraints . Engineering Cybernetics 5 , 117{ 122 ( 1984 )

9. Sergeyev , Ya., Grishagin , V.: A parallel algorithm for nding the global minimum of univariate functions . Journal of Optimization Theory and Applications 80 ( 3 ), 513 { 536 ( 1994 )

10. Sergeyev , Ya., Grishagin , V. : Sequential and parallel global optimization algorithms . Optimization Methods and Software 3 , 111 { 124 ( 1994 )

11. Sergeyev , Ya.: An information global optimization algorithm with local tuning . SIAM Journal on Optimization 5 ( 4 ), 858 { 870 ( 1995 )

12. Gergel , V. , Strongin , R.: Parallel computing for globally optimal decision making . Lecture Notes in Computer Science (including subseries Lecture Notes in Arti cial Intelligence and Lecture Notes in Bioinformatics) 2763 , 76 { 88 ( 2003 )

13. Shi , L. , lafsson , S.: Nested Partitions Method for Global Optimization . Operations Research 48 ( 3 ), 390 { 407 ( 2000 )

14. Sergeyev , Ya., Grishagin , V. :, Parallel Asynchronous Global Search and the Nested Optimization Scheme . Journal of Computational Anaysis and Applications 3 ( 2 ), 123 - 145 ( 2001 )

15. Barkalov , K. , Gergel , V:. Multilevel scheme of dimensionality reduction for parallel global search algorithms . In: Proceedings of the 1-st International Conference on Engineering and Applied Sciences Optimization, Kos Island , 2111 { 2124 ( 2014 )

16. Gergel , V , Grishagin , V. , Isra lov , R.: Local Tuning in Nested Scheme of Global Optimization. Procedia Computer Science 51 , 865 { 874 ( 2015 )

17. Gergel , V.: A software system for multiextremal optimization . European Journal of Operational Research 65 ( 3 ), 305 { 313 , ( 1993 )

18. Dam , E. , Husslage , B. , Hertog , D. : One-dimensional Nested Maximin Designs . Journal of Global Optimization 46 ( 2 ), 287 { 306 ( 2010 )

19. Gergel , V. , Kuzmin , M. , Solovyov , N. , Grishagin , V. : Recognition of Surface Defects of Cold-Rolling Sheets Based on Method of Localities. International Review of Automatic Control (IREACO) . Theory and Applications 8 ( 1 ), 52 { 55 ( 2015 )

20. Carr , C. , Howe , C. : Quantitative Decision Procedures in Management and Economic: Deterministic Theory and Applications . McGrawHill , New York, ( 1964 )

21. He , J. , Verstak , A. , Watson , L. , Sosonkina , M. : Design and implementation of a massively parallel version of DIRECT . Comput. Optim. Appl . 40 ( 2 ), 217 { 245 ( 2008 )

22. Evtushenko , Yu., Malkova , V. , Stanevichyus , A. : Parallel global optimization of functions of several variables . Computational Mathematics and Mathematical Physics 49 ( 2 ), 246 { 260 ( 2009 )

23. Ferreiro , A.M. , Garcia , J.A. , Lopez-Salas , J.G. , Vazquez , C. An e cient implementation of parallel simulated annealing algorithm in GPUs . Journal of Global Optimization 57 ( 3 ), 863 { 890 ( 2013 )

24. Barkalov , K. , Polovinkin , A. , Meyerov , I. , Sidorov , S. , Zolotykh , N. SVM Regression Parameters Optimization Using Parallel Global Search Algorithm. Lecture Notes in Computer Science , 7979 , 154{ 166 ( 2013 )

25. Grishagin , V. , Sergeyev , Ya., Strongin , R.: Parallel Characteristical Algorithms for Solving Problems of Global Optimization . Journal of Global Optimization 10 ( 2 ), 185 { 206 ( 1997 )