Numerica: A Modeling Language for Global Optimization

1 Introduction

Many science and engineering applications require the user to find solutions to systems of nonlinear constraints over real numbers or to optimize a nonlinear function subject to nonlinear constraints. This includes appli­ cations such the modeling of chemical engineering pro­ cesses and of electrical circuits, robot kinematics, chemi­ cal equilibrium problems, and design problems (e.g., nu­ clear reactor design). The field of global optimization is the study of methods to find all solutions to systems of nonlinear constraints and all global optima to optimiza­ tion problems. Nonlinear problems raise many issues from a computation standpoint. On the one hand , de­ ciding if a set of polynomial constraints has a solution is NP-hard. In fact, Canny [Canny, 1988] and Rene­ gar [Renegar, 1988] have shown that the problem is in PSPACE and it is not known whether the problem lies in NP. Nonlinear programming problems can be so hard that some methods are designed only to solve problems up to, say, 20 variables. On the other hand, computing over real numbers raises numerical problems because of the finite nature of computers.

NUMERICA [Van Hentenryck et al., 1997c] is a mod­ eling language for global optimization which makes it possible to solve nonlinear problems written in a form close to the statements traditionally found in textbooks and scientific papers. In addition, and contrary to most nonlinear programming tools, NUMERICA provides many guarantees on its results (modulo implementation er­ rors):

• Correctness: NUMERICA never produces any wrong solution;
• Completeness: Under reasonable assumptions, NUMERICA is guaranteed to isolate all solutions to nonlinear equation systems and all global optima to unconstrained and constrained optimization prob­ lems.
• Finiteness: NUMERICA is guaranteed to converge;
• Certainty: NUMERICA can prove the existence of solutions and the absence of solutions.

These functionalities should be contrasted with tradi­ tional numerical methods (e.g., quasi-Newton methods).

ear programming problems involving thousands of vari­ ables. However, there are many· applications where the additional functionalities of NUMERICA are needed, ei­ ther because of the nature of the application, or because the problem is too hard for local methods, or simply be­ cause the robustness of the approach simplifies the task. This is especially true for small-scale highly nonlinear problems as those found in chemical and electrical engi­ neering where traditional methods are likely to diverge, are unable to locate all solutions or to prove the absence of solutions (a requirement in these problems). Refer­ ence [Gehrke and Marquardt, 1996] in fact indicates that progress in chemical engineering increases the need for these functionalities.

The rest of this extended abstract illustrates the ad­ vantages of NUMERICA and contrasts it with tradi­ tional methods. More information about NUMERICA can be found in [Van Hentenryck et al., 1997c; 1997a; Michel and Van Hentenryck, 1997; Van Hentenryck et al. , 1997b].

2 What IS Possible and What IS Not?

Today's computers can manipulate and store only a fi­ nite amount of information. Since the solution of a non­ linear problem may be a real number that cannot be represented in finite space or displayed on a screen in finite time, the best we can hope for in general is a point close to a solution (preferably with some guarantee on its proximity to the solution) or an interval enclosing a solution.

Computer methods for solving nonlinear problems typically use floating-point numbers to approximate real numbers. Since there are only finitely many floating­ point numbers, these methods are bound to make nu­ merical errors. These errors, although probably small considered in isolation, may have fundamental implica­ tions on the results. Consider, for instance, Wilkinson's problem, which consists in finding all solutions to the equation

in the interval [-20.4, -9.4]. When p = 0, the equation obviously has 11 solutions. When p = 2- 23 , it has no solution. Wilkinson's problem clearly indicates that a small numerical error (e.g., assume that pis the output of some numerical computation) can have fundamental implications for the results of an application. These nu­ merical issues require users of numerical software to ex­ ercise great care when interpreting their results. With this in mind, consider the combustion problem, which consists in finding positive values for x; (1 :S i :S 10)

satisfying the equations

Using (0.5, . . . , 0.5) as starting point and the default set­ ting of the system, a well-known commercial system pro­ duces a point, say a. In the same conditions but with the defaults set to obtain the highest numerical precision, the same commercial system produces another point, say b, and prints a warning that the machine precision is not sufficient to achieve the desired accuracy. It is not obvious in this case how to interpret these results in a meaningful way.

It is also interesting to mention the common belief that proving the existence or uniqueness of solutions is outside the scope of computer algorithms. For instance, Dennis and Schnabel in their excellent text [Dennis and Schnabel, 1983] present the three functions

and state

It would be wonderful if we had a general­ purpose computer routine that would tell us: "The roots of h (x) are x = 0, 3, 4, and 5; the real roots of fz(x) are x = 1 and x := 0.888; h(x) has no real roots."

It is unlikely that there will ever be such a rou­ tine. In general, the questions of existence and uniqueness -does a problem have a solution and is it unique?- are beyond the capabilities one can expect of algorithms that solve non­ linear problems. In fact, we must readily ad­ mit that for any computer algorithm there ex­ ist nonlinear functions (infinitely continuously differentiable, if you wish) perverse enough to defeat the algorithm. Therefore, all a user can be guaranteed from any algorithm applied to a nonlinear problem is the answer, "An approx­ imate solution to the problem is ... " or "No approximate solution to the problem was found in the allocated time."

This statement is not correct in general and applies mainly to local methods. Such wonderful procedures in fact exist (within the limits imposed by the finite nature of computers) and one of them is used in NUMERICA.

Let us conclude this section by showing the behavior of NUMERICA on the above examples. On Wilkinson's

Pragma: precision= 1e-12; Constant : range idx = [1 .. 10]; Variable : x:array[1 .. 10] in [0 .. 1]; Body: solve system all x[2] + 2 * x[6] + x[9] + 2*x[10] = 1e-5; x[3] + x[S] = 3e-5 ; x[1]+x[3]+2*x[5]+2*x[8]+x[9]+x[10]=5e-5; x[4] + 2 * x[7] = 1e-5; 0.5140437e-7 * x[5] = x[1]-2; 0.1006932e-6 * x[6] = 2 * x[2]-2; 0.7816278e-15 * x[7] = x[4]-2; 0.1496236e-6 * x[S] = x[1]*x[3]; 0.6194411e-7 * x[9] = x[1]*x[2]; 0.2089296e-14 * x[10] = x[1]*x[2]-2;

problem, NUMERICA returns -20,-19, ... , -10 as solu­ tions and proves their existence when p = 0; it proves the absence of solutions when p = 2- 23 · On the com­ bustion problem whose statement is depicted in Figure 1, NuMERICA returns the unique positive solution and proves its existence in about 0.1 second. Solution b pro­ duced by the commercial system mentioned previously is close to being contained in this output box. On the functions

and for an initial range [-10 8 ·. 10 9 ], NUMERICA returns the four ranges enclosing the solutions and proves the existence of a solution in each of them for h; it returns two ranges and proves the existence of a solution in each of them for !z; it shows the absence of solutions for fg. The computation times for these examples are negligible. More precisely, the NUMERICA statement

Variable: x in [0 .. 1e8] ; Body: solve system all x-4 - 12 * x-3 + 47 * x-2 - 60 * x + 24 = 0;

produces the following output boxes:

Solution: 1 [SAFE] x = 0.8883057790717 + [0.4e-13 , 0.6e-13] Solution: 2 [SAFE] x = 1.0 + [-0.1e- 13 , +0.1e-13]

3 Local Versus Global Optimum

Traditional globally convergent methods when applied to a minimization problem converge to a local optimum

Input: int n : "Number of variables"; Constant: range idx = [1 .. n]; Variable: x : array[idx] in [-10 .. 10]; Function: y(i in idx) = 1 + 0 .25 * (x[i]-1); Body: minimize 10 * sin(pi*y(1))-2 + (y(n) - 1)-2 + Sum(i in [1 .. n-1]) (y(i) -1)-2 * (1 + 10 * sin(pi*y(i+1))-2);

from almost all starting points. They are unable however to isolate all local optima and the global optima. This limitation is well illustrated by the minimization of the function

and the maximization of the function

as well as by the minimization of the function f(xi, ... , Xn) defined as

These functions have many local minima. For instance, the last function has 10 10 local minima when n = 10 but only a single global minimum. It is unlikely that a local method will converge towards a global minimum with­ out external kno~ledge about these problems. 1 Also, a local method will never be able to prove that the global minimum has been found. In contrast, NUMERICA iso­ lates all global optima to these functions without diffi­ culty. Figure 2 is the NUMERICA statement for the last problem and it involves several of the features of the languages: input constant, minimization, function, and summation. In addition, it uses a trigonometric function sin and a predefined constant pi. NUMERICA seems to be essentially quadratic in the number of variables on this problem, as shown in Table 1.

4 Convergence

In addition to the above theoretical limitations, local methods also suffer from practical problems when im­ plemented on a computer. One of the main problems is of course convergence. An interesting example in this context is the transistor modeling example of Ebers and

1 Of course, there always exists a starting point that will con­ verge towards the global optimum.

{ (1- x1xz)x3 [ exp(x5(glk- g3kx110- 3 - g5kXs1e- 3))- 1 ]- g5k + g4kxz = 0 (1 :S k :S 4) (1- x1xz)x3 [ exp(x6(glk- gzk- g3kx110- 3 + g4kXg10- 3 ))- 1]- g5kX1 + g4k = 0 (1 :S k :S 4) X1X3 - XzX4 = 0

Moll [Ebers and Moll, 1954]. The problem is to find a solution to the system of nonlinear equations depicted in Figure 3 where the variables x; must take their values in [0, 10] and the constants are given by

The article [Ratschek and Rokne, 1993] summarizes various attempts to find a solution to this problem using local methods and states

In 1974, Cutteridge [Cutteridge, 1971] com­ bined local damped Newton-Raphson steps with the conjugate gradient method and a second­ order gradient-descent method with eigenvalue determination where the two latter methods were applied to the least squares problem [ . . . ] Cutteridge emphasized that only the sophisti­ cated combination of the three methods had led to a positive result, i.e., it did not suffice to only use the first two approaches mentioned above [. .. ]

NuMERICA finds the unique solution to the transistor modeling problem in the box [0.10] 9 and proves its exis­ tence and the absence of other solutions in less than 40 minutes. The previous interval solution required more than 14 months on a network of workstations.

Another important practical problem is convergence to an undesired solution, i.e., a solution that fails to sat­ isfy some external constraints not included in the prob­ lem statement. Globally convergent algorithms are guar­ anteed to converge to some solution or some local mini­ mum from almost any starting point, but they may fail to produce a given solution. For instance, a traditional quasi-Newton method applied to the transistor model­ ing problem almost always converges to a solution in which some variables have negative values. Solution a produced by the commercial system on the combustion problem has some negative components. Morgan [Mor­ gan, 1987] also mentions that these undesired conver­ gences are typical of chemical equilibrium systems:

The other day an electrochemist friend came by my office with a problem. He was trying to work out part of a battery-plate manufacturing pro­ cess. He had set up a math model to deter­ mine the amounts of various metal compounds that would be present in the plating bath at var­ ious times. He had ended up with a system of 10 polynomial equations in 10 unknowns. His problem was that Newton's method kept con­ verging to nonphysical solutions. [ ... ] This in­ cident has been repeated in various guises many times.

5 Practicality

The functionalities of NUMERICA of course come at a price. The intractable nature of nonlinear program­ ming precludes any guarantee on the computation times of interval methods. Conventional wisdom claims that interval methods are too slow to be of practical use and that their guarantees and ease of use come at too high a price. The performance of NUMERICA indicates that, for a rich collection of nonlinear problems, the price to pay is reasonable. Moreover, even when the full func­ tionality of global methods is not needed, NUMERICA avoids the tedious work necessary to tune local meth­ ods and find suitable starting points. As a consequence, NUMERICA's ease of use and robustness frequently com­ pensate for a longer running time and may even reduce the actual time to obtain a solution. In this context, it may be useful to mention that NUMERICA takes essen­ tially linear time in the number of variables to isolate the zeros of the Broyden banded function, a traditional benchmark from numerical analysis, even when the ini­ tial range of the variable is as large as or larger than, say, [ -10 8 , 10 8 ]. See Figure 4 for a description of the

Input:

int n : "Number of variables: "; Constant: range idx = [1 .. n]; Set : J[i in idx] = { j in [max(1,i-5) .. min(n,i+1)] Variable: x : ' array [idx] in [ -10e8 .. 10e8] ; Body: solve system all [i in idx]: j <> i }; 0 = x[i] * (2 + 5 * x[i]-2) + 1- Sum(k in J[i]) x[k] * (1 +x[k]);

Broyden banded function in NUMERICA and Table 2 for experimental results.

In addition, NUMERICA compares well and frequently outperforms continuation methods on their benchmarks. This good performance comes from a novel combination of interval analysis methods (e.g., Hansen-Sengupta's operator) and constraint satisfaction techniques. The combination of these orthogonal techniques gives sur­ prisingly good results on many problems, although un­ derstanding its strengths and limitations more formally requires further research.

Of course, there are also classes of problems for which interval methods are not appropriate at this point be­ cause interval evaluations may lose too much precision. For instance, nonlinear least-squares problems are not amenable to effective solution with the interval methods of which we are aware. Interval methods converge, of course, on these applications but they do not compare well in efficiency with local methods.

6 Challenges and Opportunities

There are many possible ways to improve global meth­ ods for nonlinear programming and we mention some of them without trying to be exhaustive. A particularly in­ teresting research avenue (studied by F. Benhamou and D. Kapur for instance) is the combination of symbolic and numerical methods. New pruning techniques with a more global view of the problem is also of paramount importance to improve the pruning when far from a solu­ tion. Similarly, it would be interesting to study ways of collecting global information beyond intervals. Finally, constraint satisfaction techniques have been a driving force behind the development of NUMERICA but only a tiny fraction of the existing research is exploited in Nu­ MERICA. It is an exciting field and it is likely to evolve substantially in the coming years.

Acknowledgment NUMERICA was developed jointly with Laurent Michel and Yves Deville. NUMERICA is also based on previous work on Newton with Frederic Benhamou, Deepak Kapur, and David McAllester. Thanks to all of them for their invaluable contributions.

Part of this research was supported by the Office of Naval Research under grant ONR Grant N00014-94-l-1153and a NSF National Young Investigator Award.