EFFICIENT COOPERATIVE SOLVERS FOR NONLINEAR ...

4 downloads 92207 Views 245KB Size Report
the form of trigonometric functions, which are similar to problems appearing in robot control ... nally, in Section 5, we conclude the paper and discuss future work . ... In this paper, when max = 1 we call it static penalty method (SPM), i.e.,.
EFFICIENT COOPERATIVE SOLVERS FOR NONLINEAR CONTINUOUS CONSTRAINT PROBLEMS  YI SHANG, MARKUS P.J. FROMHERZ, AND LARA S. CRAWFORD

Palo Alto Research Center 3333 Coyote Hill Road Palo Alto, CA 94304

fyshang,fromherz,[email protected]

In the past decade, signi cant progress has been made in understanding problem complexity of discrete constraint problems. In contrast, little similar work has been done for constraint problems in the continuous domain. In this paper, we study the complexity of several popular methods for nonlinear constraint problems and present cooperative solvers with improved performance. The methods include a sequential quadratic programming (SQP) method, a penalty method with a xed penalty function, a penalty method with a sequence of penalty functions, and an augmented Lagrangian method. In our experiments, we apply these methods to solve a series of continuous constraint problems with increasing constraint-tovariable ratios and obtain novel results on complexity phase transition phenomena. Speci cally, for constraint satisfaction problems, the SQP method is the best on weakly constrained problems, whereas the augmented Lagrangian method is the best on highly constrained ones. Although the static penalty method performs poorly by itself, by combining it with the SQP method, we show a cooperative solver that is signi cantly better than any of the individual methods on problems with moderate to large constraint-to-variable ratios.

1 Introduction

Many important real-world applications in domains such as control, planning, and diagnosis of physical systems can be regarded as constraint satisfaction (CSPs) and constrained optimization problems (COPs) 1 . E ective constraint-based techniques should be aware of the complexity of real-world constraint problems by dynamically selecting and adapting solvers to the structure of the problem. Toward this end, this paper studies the complexity of solution methods on continuous constraint problems. In the past decade, signi cant progress has been made in understanding problem complexity of discrete constraint satisfaction and combinatorial optimization problems, such as satis ability (SAT), planning, graph coloring, and the traveling salesman problem. One particularly well-studied phenomenon is the complexity phase transition 2 . For example, in solving randomly gen THIS WORK HAS BEEN SUPPORTED IN PART BY DARPA UNDER CONTRACT F33615-01-C-1904. PRESENTED AT ICMS'02, BEIJING, CHINA, AUG. 2002.

Submitted to World Scienti c on April 8, 2002.

1

erated SAT problems, the average time to nd a solution or determine that none exists is short when the ratio of clauses to variables is small (i.e., the problem is weakly constrained) as well as when this ratio is large (where the problem is highly or overconstrained), while solving time is largest in the intermediate case. This phenomenon has been observed for a variety of complete and incomplete search algorithms, including the Davis-Putnam method 3 and GSAT 4 . Towards understanding the complexity of continuous constraint problems, in this paper, we empirically study the complexity phase transition phenomena of several popular methods in solving these problems with varying constraintto-variable ratios. The methods include a sequential quadratic programming (SQP) method, a penalty method with a xed penalty function, a penalty method with a sequence of penalty functions, and an augmented Lagrangian method. They represent major approaches for nonlinear constraint satisfaction and constrained optimization problems: solving a sequence of quadratic programming sub-problems, or minimizing one or a sequence of unconstrained function. The test cases consist of randomly generated nonlinear constraints in the form of trigonometric functions, which are similar to problems appearing in robot control applications. In our experiments, we study the complexity of the methods in solving a series of continuous constraint problems with increasing constraint-to-variable ratios. In addition, we propose cooperative solvers that combine multiple methods and show their performance improvement. The rest of this paper is structured as follows. In Section 2, we rst introduce individual methods, and then discuss cooperative solvers. In Section 3, we present the test-case generators for generating continuous constraint problems used in our study. In Section 4, we show experimental results on problems with di erent constraint-to-variable ratios (or constraint ratios). Finally, in Section 5, we conclude the paper and discuss future work. 2 Constraint Solvers

Di erent methods have di erent complexity behaviors in solving continuous constraint problems. In this section, we introduce some popular methods and discuss alternatives of combining multiple methods in a cooperative solver. 2.1 Quadratic Penalty Method

One fundamental approach to constrained optimization is to replace the original problem by a penalty function or a sequence of penalty functions 5 . The approach is popular because of its simplicity and intuitive appeal. One of Submitted to World Scienti c on April 8, 2002.

2

the simplest penalty function is the quadratic penalty function. For the inequality-constrained problem (the type of problem we tested in the experiments) min f (x) subject to ci (x)  0; i = 1;    ; m (1) x where x = (x1 ; x2 ; : : : ; xn ), the quadraticmpenalty function Q(x; ) is 1 X(max(0; c (x)))2 Q(x; ) = f (x) + (2) i 2 i=1

where  > 0 is the penalty parameter. Based on the penalty function, the quadratic penalty method is (Section 17.1 of Numerical Optimization 5 ) Given 1 > 0, r 2 (0; 1), starting point x1 for k = 1;    ; kmax Find a minimizer xk of Q(x; k ), starting at xk if a solution satisfying the requirement is found, stop with xk k+1 = r  k , xk+1 = xk 0

0

0

end

In this paper, when kmax = 1 we call it static penalty method (SPM), i.e., solving a xed penalty function. When kmax > 1, we call it dynamic penalty method (DPM), i.e., solving a sequence of penalty functions with di erent penalty parameters. In our experiments, we use a quasi-Newton solver, from the Matlab Optimization Toolbox, to nd a minimizer of Q(x; k ) in each iteration. For nonlinear constrained problems, the penalty method may be trapped by local optima and may not be able to nd a feasible solution. When that happens, we restart it from another random point. fminunc

2.2 Augmented Lagrangian Method

The augmented Lagrangian method (ALM) improves the quadratic penalty function by using explicit Lagrange multiplier estimates to avoid illconditioning. For the inequality-constrained problem in Eq. (1), the augmented Lagrangian function LA (x; ; ) is (Section 17.4 of Numerical Optimization 5 ) m X LA (x; ; ) = f (x) + (ci (x); i ; ) (3) i=1

where  = (1 ;    ; m ) are Lagrange multipliers,  > 0 is the penalty parameter, and  1 2 ci (x) + i  0 (ci (x); i ; ) = i ci(x2 ) + 2 ci (x) ifotherwise (4) 2

i

Submitted to World Scienti c on April 8, 2002.

3

Thus, the augmented Lagrangian method is Given 1 > 0, r 2 (0; 1), starting point x1 and 1 for k = 1;    ; kmax Find a minimizer xk of LA (x; k ; k ), starting at xk if a solution satisfying the requirement is found, stop with xk Update Lagrange multipliers: ki +1 = max(ki + ci (xk )=k ; 0) k+1 = r  k , xk+1 = xk 0

0

0

0

end

In our experiments, we use to nd a minimizer of LA(x; ; ) in each iteration. For nonlinear constrained problems, the augmented Lagrangian method may not be able to nd a feasible solution. When that happens, we restart it from another random point. fminunc

2.3 Sequential Quadratic Programming Method

Sequential quadratic programming (SQP) methods are among the most e ective methods for nonlinearly constrained optimization. The essential idea of SQP is to model the constrained problem at each iteration as a quadratic programming subproblem and to use the minimizer of this subproblem to de ne a new iteration. One example is to apply Newton's method to the Karush-KuhnKucker optimality conditions for the constrained problem. This iteration is repeated until a constrained locally optimal solution is found. The SQP solver used in our experiments is from the Matlab Optimization Toolbox. An iteration consists of three main stages: updating of the Hessian matrix of the Lagrangian function, quadratic programming problem solution, and line search and merit function calculation. When fails to nd a feasible solution, we restart it from another random point in the experiments. fmincon

fmincon

fmin-

con

2.4 Cooperative Solvers

To take advantage of the respective strengths of multiple methods, various forms of meta-heuristics have been proposed to combine them. One popular approach is to combine global and local search methods, in which local search may follow global search in sequence, or local search may be used as a re nement component in global search. The goal of the former strategy is to let the global solver move towards feasible regions or global optimal solutions, from which the local solver converges quickly to local minima 6 . The second strategy is for the global solver to use a local solver to improve its intermediate points to true local optima 7 . Submitted to World Scienti c on April 8, 2002.

4

Another approach is to combine di erent types of methods, such as penalty methods and SQP methods. This is the approach we take in this paper. In our experiments on constraint satisfaction problems, we found that the static penalty method can nd solutions with small penalty values very quickly, even for highly constrained problems. However, it has diÆculty in nding a solution that satis es all constraints simultaneously. On the other hand, the SQP method are very e ective when starting from an initial points close to the feasible region. Therefore, we propose a cooperative solver (called SPM+SQP) that combines the static penalty method and the SQP method to solve CSPs as follows: 1. Run the static penalty method from a starting point xs and obtain a solution x0 . If x0 is a feasible solution, then stop. 2. Run the SQP method from x0 and obtain a solution x00 . If x00 is a feasible solution, then stop. 3. Generate a new random starting point and go to step 1. In a similar way, we propose cooperative solvers DPM+SQP and ALM+SQP that combine the dynamic penalty method or the augmented Lagrangian method with the SQP method, respectively. The procedure of SPM+SQP is slightly di erent for COPs as follows: 1. Run the static penalty method from a starting point xs to minimize constraint violation and obtain a solution x0 . 2. Run the SQP method from x0 and obtain a solution x00 . If x00 is a feasible local minimum, then stop. 3. Generate a new random point and go to step 1. The di erence is that for COPs, SQP is always run to optimize the objective since the static penalty method only nds a feasible solution. Another variation of the procedure is to let the static penalty method minimize the penalty function in Eq. (2) that includes the original objective. 3 Test Problems

Since our performance analysis is experimental in nature, it is important to choose appropriate problem ensembles for testing algorithms, as is done in the study of discrete problems 2 . Regarding continuous constraint problems, there is an unful lled need for generic benchmark problem sets, in spite of the work in the mathematical programming and evolutionary computing communities. In the mathematical programming eld, test problem sets such as CUTE (Constrained and Unconstrained Testing Environment) 8 and COPS Submitted to World Scienti c on April 8, 2002.

5

(Constrained Optimization ProblemS) 1 provide very limited ability to generate problems with di erent characteristics by changing problem parameters. In the evolutionary computing eld, a exible test-case generator proposed recently 9 has the major limitation that it de nes a landscape that consists of a collection of functions that are only piece-wise continuous and are de ned on di erent subspaces of equal size. Because derivatives do not exist on the boundary of two subspaces, methods requiring continuous derivatives do not work well on these problems. To study constraint problems of di erent sizes with various numbers and types of constraints, we need a scalable test suite, in which many of the problem parameters can be changed easily. It should allow one to easily compare di erent optimization techniques developed in di erent elds. It should be useful in investigating the relationships between various characteristics (e.g., constraint ratio) and complexity of constraint problems, as well as the computational complexity of various solving techniques. In this section, we present a exible test-case generator for continuous constraint problems. It is based on trigonometric functions and generates constraint problems similar to those appearing in robot control applications. It also borrows ideas from randomly generated SAT problems and generates constraints using either a xed-clause-length model or a constant-density model 3 . In the xed-clause-length model, each constraint has a xed number of variables that are randomly chosen. In contrast, in the constant-density model, each variable may appear in a constraint with some equal probability. Our test-case generator G(n; m; ks ; kp ; ps ; pp ; t) is mainly parameterized to generate di erent forms of constraints, with optional objectives. In this paper, the objective f (x) of the test cases is either 0 or a quadratic function n X min f (x) = (xi ai )2 (5) x i=1

i

where ai is either a constant or a random number. When f (x) is 0, the problem becomes a constraint satisfaction problem. The constraints of the test cases are generated according to a single base function with random coeÆcients. For the experiments in this paper, the base function of the test cases is a nonlinear function, in the form of sums of products of sine and cosine functions: NS NP X Y  1 c sin 2x + 2Æ  (6) i

PNS QNPi

i=1

j =1 cij i=1 j =1

ij

rij

ij

where rij 2 f1; : : : ; ng is the index of a randomly chosen variable; 1  cij  1 is a random number specifying the magnitude of the sine function; 0  Æij  1 Submitted to World Scienti c on April 8, 2002.

6

is a random number specifying the phase; N S is the number of sum terms; N Pi is the number of product terms in each sum term; and 1    1 speci es the threshold for satisfying the constraint. Based on constraints de ned in Eq. (6), the parameters of the generator G(n; m; ks ; kp ; ps ; pp ; t) are the number of variables (n), the number of constraints (m), the maximum number of sum terms in a constraint (ks ; ks  1), the maximum number of product terms in any sum term of a constraint (kp ; kp  1), the probability determining the number of sum terms (up to ks ) in a constraint (ps ), the probability determining the number of product terms (up to kp ) in a sum term (pp ), and the average tightness of the constraints (t; 0  t  1). The tightness t of a (continuous) constraint is de ned as the ratio of the size of the space made infeasible by the constraint to the size of the total problem space. For example, t = 0 means that the constraint leaves the entire problem space feasible, and t = 1=2 means that it makes half of the problem space infeasible. The parameter  in Eq. (6) is determined based on the tightness t. The  value that achieves the average tightness speci ed by t is found experimentally through random sampling. All the parameters of the generator can be changed independently. Constraints are generated with randomly chosen coeÆcients. By changing the tightness parameter t we can control the diÆculty of the constraints, which in turn a ects the complexity of the whole problem. 4 Experimental Results

In the experiments, the test cases are generated by G(n; m; 3; 3; 1; 1; t). A simple bound 1  x  1 is imposed. To guarantee a randomly generated problem has feasible solutions, we let the origin be a feasible point of all its constraints. When we generate a constraint for which the origin is not a feasible point, the constraint is simply discarded. Since the constraint ratio has been identi ed as a key indicator of complexity for discrete CSPs, in this paper we compare the performance of solvers on continuous problems with varying constraint ratios. For each constraint ratio, 100 random instances are generated and the performance of the solvers is compared based on the median numbers of function evaluations. We use medians instead of means because the means of the number of function evaluations are heavily in uenced by a small number of extremely large values. Our interest is in what the majority of instances from the distribution are like, not in the unusual or extreme cases. As the median is more robust in the presence of such outliers, it appears to be a more informative statistic for our current purpose. Submitted to World Scienti c on April 8, 2002.

7

4.1 Results for CSPs

This set of experiments is based on CSPs with n = 25 and 50, t = 1=4 and 1=8, and varying m. The initial and restart points are generated with the following pre-sampling scheme: n random points are sampled within the bounds of the constrained problems, i.e., [ 1; 1]n, and the best of them is chosen as the initial or restart point. The best sample point is the one with the smallest constraint violation, which may be close to or even inside a feasible region. The pre-sampling scheme improves the solvers' performance over just picking a single random point within the bounds 6 . Four individual methods (SQP, static penalty, dynamic penalty, and augmented Lagrangian method) and three cooperative solvers (SPM+SQP, DPM+SQP, and ALM+SQP) are used to solve the problems. In the static penalty method,  = 1=2. In both the dynamic penalty method and the augmented Lagrangian method, kmax = 5, 1 = 10, and r = 0:1. 1 = 0 in the augmented Lagrangian method. All methods run using only constraint and/or objective values, without user supplied gradients. Figure 1 shows the numbers of function evaluations and restarts used by the four individual methods in solving 25 and 50 variable instances of tightness t = 1=4 or 1/8. Their complexity behaviors are quite di erent. In terms of number of restarts, they all have a phase transition from at, or 1, to exponential. The SQP and the static penalty method have similar transition points, whereas the dynamic penalty and augmented Lagrangian method have similar transition points, which are much larger than the ones of the SQP and the static penalty method. In the exponential phase, the number of restarts used by the static penalty method increases at a much slower rate than that of SQP, but much larger than those of the dynamic penalty and the augmented Lagrangian method. In terms of number of function evaluations, they also have phase transitions from polynomial to exponential, corresponding to their phase transitions in number of restarts. For problems of smaller constraint ratios, SQP uses the fewest number of function evaluations, while the other three methods are similar. On the other hand, for problems of larger constraint ratios, SQP is the worst, the static penalty method is in the middle, and the augmented Lagrangian and the dynamic penalty method are much better. Similar results were obtained on 50 variable instances. For these methods, when the tightness parameter doubles, i.e., from 1/8 to 1/4, the problems become harder and the number of function evaluations needed is nearly doubled. Figure 1 also shows the performance improvement of the cooperative solver SPM+SQP over the individual methods on highly constrained probSubmitted to World Scienti c on April 8, 2002.

8

n

= 25; t = 1=8 3

5

10

10

SQP Static Penalty (SPM) Dynamic Penalty (DPM) Augmented Lagrangian (ALM) SPM+SQP

4

10

2

Median # of restarts

Median # of function evaluations

SQP Static Penalty (SPM) Dynamic Penalty (DPM) Augmented Lagrangian (ALM) SPM+SQP

3

10

10

1

10

2

10

0

1

2

3

4

5

6

7 8 9 Constraint ratio

10

11

12

13

n 5

10

14

3

4

5

6

7 8 9 Constraint ratio

10

11

12

13

14

4

10 SQP Static Penalty (SPM) Dynamic Penalty (DPM) Augmented Lagrangian (ALM) SPM+SQP

SQP Static Penalty (SPM) Dynamic Penalty (DPM) Augmented Lagrangian (ALM) SPM+SQP

3

10 4

Median # of restarts

10

3

10

2

10

1

10

2

10

0

1

2

3

4 Constraint ratio

5

6

10

7

n 5

1

3

4 Constraint ratio

5

6

7

3

10 SQP Static Penalty (SPM) Dynamic Penalty (DPM) Augmented Lagrangian (ALM) SPM+SQP

SQP Static Penalty (SPM) Dynamic Penalty (DPM) Augmented Lagrangian (ALM) SPM+SQP

4

2

Median # of restarts

10

3

10

2

10

2

= 50; t = 1=8

10

Median # of function evaluations

2

= 25; t = 1=4

10

Median # of function evaluations

1

10

1

10

0

1

2

3

4

5 6 Constraint ratio

7

8

9

10

10

1

2

3

4

5 6 Constraint ratio

7

8

9

10

Figure 1. Complexity of four individual methods and one cooperative solver SPM+SQP in solving CSPs of 25 and 50 variables. The constraint tightness parameter t is 1/8 or 1/4. SPM+SQP is the best on highly constrained problems.

Submitted to World Scienti c on April 8, 2002.

9

lems. Compared with the best individual method, the augmented Lagrangian method, SPM+SQP has about the same number of restarts, but uses much fewer function evaluations. Similarly we have tried two other cooperative solvers, DPM+SQP and ALM+SQP. They improve over the dynamic penalty method and the augmented Lagrangian method, respectively, on highly constrained problems by reducing the number of restarts needed and the number of function evaluations. 4.2 Results for COPs

In this section, we show results on COPs with a simple quadratic objective function in Eq. (5), where each ai is a random number in [-0.5, 0.5]. The constraints are generated in the same way as in the last section. Figure 2 compares the complexity of SQP and SPM+SQP on COPs with increasing constraint ratios. In terms of function evaluations, their cost curves have a up-down-up shape. One explanation is that the initial up phase happens when feasible regions are large and easy to nd. The solvers spend a long time on optimizing the objective. As the constraint ratio increases, the boundaries of feasible regions become more complicated, requiring more effort in nding a feasible local minimum. In the down phase, feasible regions get fewer and smaller, and their boundaries may become simpler. It becomes easier to nd a feasible local minimum. No restart is necessary so far. In the second up phase, feasible regions become rare and hard to nd. The solvers used increasingly more restarts. The cost curves of the two solvers are almost identical on problems with small constraint ratios. The cooperative solver is slightly better on problems with large constraint ratios. They are similar in number of restarts. 5 Conclusions

In this paper, we study the complexity of several representative methods for continuous constraint problems and obtain various novel observations. By combining two methods, we propose some cooperative solvers that improve over existing methods on harder problems by reducing the number of restarts. We show that problem parameters such as constraint ratio and tightness are important indicators of the diÆculty of the constraint problems. Based on these parameters, adaptive solvers that select appropriate methods dynamically and automatically may achieve the best-possible performance for the problem at hand. Submitted to World Scienti c on April 8, 2002.

10

5

10

1

10

4

10

Median # of restarts

Median # of function evaluations

SQP, t=1/4 SQP, t=1/8 SPM+SQP, t=1/4 SPM+SQP, t=1/8

3

10

SQP, t=1/4 SQP, t=1/8 SPM+SQP, t=1/4 SPM+SQP, t=1/8

2

10

0

0

5

10 Constraint ratio

15

10

0

5

Constraint ratio

10

15

Figure 2. Comparison of SQP and SPM+SQP on 25-variable COPs, where the tightness parameter t = 1=8 or 1=4.

References

1. A. S. Bondarenko, D. M. Bortz, and J. J. More. COPS: Large-scale nonlinearly constrained optimization problems. Technical Memorandum ANL/MCS-TM-237, Argonne National Laboratory, Argonne, Illinois, 1998. 2. Tad Hogg, Bernardo A. Huberman, and Colin Williams. Phase transitions and the search problem. Arti cial Intelligence, 81:1{15, 1996. 3. Bart Selman, David Mitchell, and Hector J. Levesque. Generating hard satis ability problems. Arti cial Intelligence, 81:17{29, 1996. 4. M. Yokoo. Why Adding More Constraints Makes a Problem Easier for Hill-Climbing Algorithms: Analyzing Landscapes of CSPs. In Proc. of CP'97, number 1330 in LNCS, pages 357{370. Springer-Verlag, 1997. 5. J. Nocedal and S. J. Wright. Numerical Optimization. Springer, 1999. 6. Y. Shang, M. P. J. Fromherz, T. Hogg, and W. B. Jackson. Complexity of continuous, 3-SAT-like constraint satisfaction problems. In Proc. IJCAI01 Workshop on Stochastic Search Algorithms, pages 49{54, Seattle, WA, 2001. 7. W. E. Hart. Adaptive global optimization with local search. PhD thesis, University of California, San Diego, 1994. 8. CUTE. Constrained and unconstrained testing environment, 2001. http://www.cse.clrc.ac.uk/Activity/CUTE. 9. Z. Michalewicz, K. Deb, M. Schmidt, and T. Stidsen. Test-case generator for constrained parameter optimization techniques. IEEE Transactions on Evolutionary Computation, 4(3):197{215, September 2000.

Submitted to World Scienti c on April 8, 2002.

11