MathOptimizer Professional: 
New Features and Applications

Presented at the 2006 Wolfram Technology Conference

Frank Kampas
WAM Systems Inc., Plymouth Meeting, PA, USA
fkampas@msn.com

János D. Pintér
Pintér Consulting Services Inc., Halifax, NS, Canada
jdpinter@hfx.eastlink.ca
www.pinterconsulting.com

Introduction

Overview

MathOptimizer Professional (MOP) is a third party product that has been in development and use since 2003 (Pintér and Kampas, 2003). MOP solves constrained nonlinear – global and local – optimization models, under very general structural requirements. The general model form handled by MOP is

min f(x)   subject to xϵD   D:={x: xl≤x≤xu  g(x)≤0}⊂(R^n).

Here xϵ(R^n) is the vector of decision variables; f:(R^n)->R is the continuous objective function; D is the set of feasible decisions. The (by assumption, non-empty) set D is defined by finite (component-wise) lower and upper bounds xl and xu, and by a collection (m-vector) of continuous constraint functions g:(R^n)->(R^m).

This model-class is very general and includes great many difficult instances. For a simple illustration, consider the model (in Mathematica-style notation)

min (Sin[x*y - y]^2 + 
            Cos[x - (y^2)]^2)  {x + Sin[y] ≤ 
            0.3,  (x^3) y - (y^2) ≤ 0}, {{x, -5, 
          2}, {y, -1, 8}}

The objective function of this model instance is shown below. (Similar figures can be generated for the constraints, including combined function plots.)

Plot3D[Sin[x*y - y]^2 + Cos[x - (y^2)]^2, {x, -5, 
            2}, {y, -1, 8}, ColorFunction -> Hue, 
          PlotPoints -> 200,  Mesh  ->  False,  
          ImageSize  ->  500];

Obviously, local search will not work in such – or even in much simler – models, and hence one needs to apply suitable global scope search techniques.

MathOptimizer Professional combines the model development power of Mathematica (Wolfram, 2003; Wolfram Research, 2006) with the Lipschitz Global Optimizer (LGO) solver suite (Pintér, 1996, 2002, 2005). MOP automatically converts the optimization model formulated in Mathematica, making use of a (supported) C or Fortran compiler. The generated model code is then linked to the LGO external solver engine. The numerical solution produced by LGO is seamlessly reported back to the calling Mathematica document. (The user does not need to write or manipulate the interim C or Fortran code and the generated LGO input/output files.) The functionality outlined above is summarized by the next figure.

MathOptimizer Professional offers sophisticated application development tools and an external solver engine based functionality comparable to other efficient – compiler-based or optimization modeling language-related – implementations. For further key technical information and professional details related to MOP, consult the references, and visit  www.pinterconsulting.com and www.wolfram.com/products/applications/mathoptpro/.

Topics Covered

In this presentation, we illustrate several new features of MathOptimizer Professional. These features include the following:

• Increased solution speed

• Increased solvable model sizes: in default usage mode, MOP can now handle models with up to 2000 variables and 2000 general constraints, in addition to the explicit variable bound constraints

• Support for an extended range of C and Fortran compilers, for Windows and Linux platforms

• New option: support for repeating a MathOptimizer Professional program run, without the need to recompile the generated C or Fortran code

• New option: ShowMOPProgress

• New option: Auxiliary Functions

• New technique: Logical Combination of Constraints

We also review several recent challenging and large-scale MOP applications, including further application perspectives. For more details regarding these, consult (Castillo, Kampas and Pintér, 2006; Kampas and Pintér, 2004a, b, 2006; Mason et al., 2006; Pintér and Kampas, 2005; Pintér, 2006).

MathOptimizer Professional: Highlighted Features

Supported Compilers

As outlined above, MathOptimizer Professional provides a Mathematica interface to the LGO solver suite. The Windows version of LGO requires an input file (referred to as LGO.in): this text file describes fixed characteristics of the problem such as model/variable/constraint names, number of variables, number and type of constraints, and solver options settings. In addition to this information, LGO iteratively communicates with a dynamic link library (dll) which evaluates the objective function and the constraints for algorithmically chosen settings of the variables. Frequently, global optimization requires a large number of such model function evaluations. Windows dlls can be generated by a number of FORTRAN and C compilers and return results to LGO with very little overhead. We currently have versions of MOP for the Intel, Lahey and Salford FORTRAN compilers and for Microsoft Visual C 6.0, Visual Studio 2003 and Visual Studio 2005 C compilers, and the Borland C compiler.  The free SDKs (software development kits) for Visual Studio 2003 and 2005 are also supported. There is also a version of MOP for 64 bit Linux, using the GFORTRAN compiler.  Under Linux, a FORTRAN program is compiled and linked to an LGO object file.

MOP Operations: Illustrative Examples

Key MathOptimizer Professional Steps

When performing an optimization, MOP (Windows version) carries out the following steps:

1.  Put the problem into a standard form.

2.  Generate and export the lgo.in file.

3.  Generate C or FORTRAN source code giving the objective function and constraint violations (using Mark Sofroniou's package Format.m).

4.  Call the external compiler to generate the dll.

5.  Call LGO, which calls the dll thousands to millions of times.

6.  Read LGO.sum, which gives the result of the optimization, and convert the variable values to Mathematica rules.

A MathLink executable linking to the C function syscall is used to perform the system calls (since the Mathematica Run function flashes a DOS window on the screen, which could become annoying).

This process is illustrated by a calculation which finds the smallest sphere into which spheres of radii Range[10]^(-1/3) can be packed.  The Salford Fortran95 compiler (available as a free download for personal use) is used in this example. Note that this formulation of the problem would not work with a C compiler since the C max function takes only 2 arguments, whereas the FORTRAN version can take any finite integer number.

Activate and Check MOP Functionality

Needs["MathOptimizerProS95F`callLGO`"];

Options[callLGO]

As a functionality test, let us solve the model defined in the introductory section. The numerical solution found (below) is obviously a global optimum; there could be also other global solutions.

callLGO[
      Sin[x*y - y]^2 + 
        Cos[x - (y^2)]^2,  {x + Sin[y] ≤ 
          0.3,  (x^3) y - (y^2) ≤ 0}, {{x, -5, 
          2}, {y, -1, 8}}]

A Non-uniform Sphere Packing Problem in Three Dimensions

We are given a set of three-dimensional spheres. We wish to find (numerically estimate) the radius of the smallest sphere that includes all the given spheres, in a mutually non-overlapping configuration. This model extends our earlier studies related to the planar equivalent of this problem (Kampas and Pintér, 2004), (Pintér and Kampas, 2005). The model development is briefly described below.

outerdist[r] is the maximum distance from the origin to a sphere of radius r, centered at {x[r],y[r],z[r]}

outerdist[
        r_] := r +(Sqrt[x[r]^2 + y[r]^2 + 
              z[r]^2])

pairdist[{r1,r2}] is the minimum distance between two spheres with radii r1 and r2, centered at { x[r1], y[r1], z[r1] } and { x[r2], y[r2], z[r2]}

pairdist[{r1_, 
          r2_}] :=   (Sqrt[(((x[r1] - x[r2]))^2) + (((y[r1] - y[r2]))^2) + (((z[r1] - z[r2]))^2)]) - r1 - r2

packspheres[r, max, opts] finds the optimized packing of spheres with radii ri inside a sphere of radius r.  argument max is an upper bound on the radius of the enclosing sphere. The results are the radius of the enclosing sphere, rules giving the sphere center positions and the maximum constraint violation.

packspheres[r_List, max_, 
        opts___] := 
Module[
{num  =  
            Length[r], maxouter, pairs,  cons, vars, bds, lb, ub, soln, 
          rcirc, pd}, 
vars = 
          Join[x  /@  r, y  /@  r,  
            z  /@  
              r]; 
bds  =  (max - #) &  /@  
            r; 
ub  =  
          Join[bds, bds, 
            bds]; 
lb = -ub; 
vwb = 
          Transpose[{vars, lb, ub}]; 
maxouter  =  
          Max[outerdist  /@  r]; 
pairs  =  
          Select[Flatten[ Outer[List, r, r], 
              1], #[[1]] < #[[2]] &]; 
cons =  
          Join[Thread[
              pairdist  /@  
                  pairs  ≥    0], {maxouter  ≤ 
                max}]; 
callLGO  @@  {maxouter, cons, 
            vwb, opts}
]

showspheres[soln] generates a 3D graphic of the packspheres solution, with the enclosing sphere shown as a wire frame.

<< Graphics`Shapes`

showspheres[
        soln_] := 
Module[
{r  = 
            Cases[soln[[2, All, 1]], 
              x[n_] ->  n], t1, t2, p1, 
          p2}, 
t1  =  
          Graphics3D @ 
            Table[TranslateShape[
                Sphere[r[[i]], 10, 
                  10], {x[r[[i]]], 
                    y[r[[i]]], 
                    z[r[[i]]]} /. 
                  soln[[2]]] , {i, 
                Length[r]}]; 
t2  =  
          WireFrame[
            Graphics3D[
              Sphere[soln[[1]], 20, 
                20]]]; 
p1 =  
          Show[t1, 
            DisplayFunction  ->  
              Identity]; 
p2  =  
          Show[t2,  
            DisplayFunction  ->  Identity]; 
Show[
          p1, p2,  DisplayFunction  ->  $DisplayFunction,  
          Boxed  ->  False]]

AbsoluteTiming[
      s4 = packspheres[Range[10]^(-1/3), 2.25, 
          MaxEvals  -> (10^4)]]

showspheres[s4];

Option ShowMOPProgress (New)

To illustrate MOP features, including the new options, we will solve instances of the 3D non-uniform sphere packing problem formulated above. The next statement solves a 10-sphere instance. As indicated below, the compilation of the model code is the most time-consuming operation listed. The external LGO solver runtime will be retrieved in the next section.

AbsoluteTiming[
      packspheres[Range[10]^(-1/3), 2.25, MaxEvals  -> (10^4), 
        ShowMOPProgress  ->  True]]

Options ShowLGOInputs and ShowSummary

The LGO input and output files can be seen by setting these options to True. For our example, see the automatically generated Fortran code of the model function file, the LGO.in file, and the LGO.sum summary result file below. Observe that this fairly dificult model has 30 decision variables, and 45 non-convex constraints, and that it is solved in seconds. The results also show that all constraints are met with more than (10^(-8)) precision, see the "Constraint function values at optimum estimate" entries in the summary result file below.

packspheres[Range[3]^(-1/3), 2.25, MaxEvals  -> (10^4), 
      ShowLGOInputs  ->  True,  
      ShowSummary  ->  True]

Larger Models Handled (New)

The LGO solver engine in MOP now handles models with up to 2000 variables and constraints.

vars = Array[x, {2000}];

The next line solves the simple test problem min (∑_(i = 1)^2000) ((Subscript[x, i])^2)    -1≤(Subscript[x, i])≤3. We illustrate this with a local search and will extract the solution after the next line is executed, to avoid seeing all 2000 variable results.

AbsoluteTiming[soln  =  
          callLGO[Total[vars . vars], {}, Thread[{vars, -1, 3}],  
            Method -> LS];]

soln[[1]]

soln[[2, 100]]

soln[[2, 1000]]

Option Recompile (New)

The C or FORTRAN dll need not be recompiled if only the optimization options are changed.

AbsoluteTiming[packspheres[Range[10]^(-1/3), 2.25, 
          MaxEvals  -> 
            (10^4)][[1]]]

For this calculation, nearly half of the elapsed time went into managing the compilation of the source code for the dll.

AbsoluteTiming[packspheres[Range[10]^(-1/3), 2.25, 
          MaxEvals  -> (10^4), 
          Recompile  ->  
            False][[1]]]

Option LGORandomSeed (New)

Using multiple random seeds is often an effective way of getting a better solution. If multiple random seeds are used, the model is only compiled for the first run. The total time for doing 10 random seeds is approximately 10 times larger than the elapsed time for the evaluation without recompilations. Notice the best result found (an embedding circle with radius 1.88297).

AbsoluteTiming[packspheres[Range[10]^(-1/3), 2.25, 
          MaxEvals  -> (10^4), 
          LGORandomSeed  ->  2*Range[10]][[All, 
        1]]]

New Built-in Stopping Mechanism

The maximum number of model function evaluations in the global search phase is now enforced more strictly. Observe the improved numerical values resulting from more global search effort.

AbsoluteTiming[
      First @ (s5 = 
            packspheres[Range[10]^(-1/3), 2.25, 
              MaxEvals  -> (10^5), 
              Recompile  ->  False])]

showspheres[s5];

AbsoluteTiming[
      First @ (s6 = 
            packspheres[Range[10]^(-1/3), 2.25, 
              MaxEvals  -> (10^6), 
              Recompile  ->  False])]

showspheres[s6];

AbsoluteTiming[
      First @ (s7 = 
            packspheres[Range[10]^(-1/3), 2.25, 
              MaxEvals  -> (10^7), 
              Recompile  ->  False])]

showspheres[s7];

Reporting Intermediate Results

Now a more useful result is obtained and reported, when the set TimeLimit is reached.

AbsoluteTiming[packspheres[Range[10]^(-1/3), 2.25, 
          MaxEvals  -> (10^7), Recompile  ->  False,  
          TimeLimit  ->  
            30][[1]]]

Definition of Auxiliary Functions (New)

Auxiliary functions can be defined to simplify model building (currently, in FORTRAN compiler based versions only). Functions of more than one variable can be defined, although that is not illustrated with the example shown below. This example is taken from Trefethen's hundred-dollar, hundred-digit numerical challenge (Trefethen, 2002); for solutions, see for example the topical entry of MathWorld (Weisstein, 2006) with references.

callLGO[
      Exp[Sin[50*x]] + Sin[60*Exp[y]] + Sin[70*Sin[x]] + Sin[Sin[80*y]] - 
        Sin[10*(x + y)] + (((x^2) + (y^2))/4), {}, {{x, -3, 
          3}, {y, -3, 3}},  MaxEvals  ->  (10^5)]

callLGO[
      e[s[50*x]] + s[60*e[y]] + s[70*s[x]] + s[s[80*y]] - 
        s[10*(x + y)] + (((x^2) + (y^2))/4), {}, {{x, -3, 
          3}, {y, -3, 3}}, 
      AuxFuncs  ->  {{e, {x}, Exp[x]}, {s, {x}, Sin[x]}}, 
      MaxEvals  ->  (10^5)]

By looking at the FORTRAN source code, we see how the auxiliary functions are coded (generated).

callLGO[
      e[s[50*x]] + s[60*e[y]] + s[70*s[x]] + s[s[80*y]] - 
        s[10*(x + y)] + (((x^2) + (y^2))/4), {}, {{x, -3, 
          3}, {y, -3, 3}}, 
      AuxFuncs  ->  {{e, {x}, Exp[x]}, {s, {x}, Sin[x]}}, 
      MaxEvals  ->  (10^5),  
      ShowLGOInputs  ->  True]

Logical Combinations of Constraints (New)

This isn't really based on a new feature of MOP, but a mental breakthrough on our part.  In short, And is equivalent to taking the maximum of the constraint violations; Or is equivalent to taking the minimum of the constraint violations.

cnvrt[cons_List] :=  
      Module[{r1 =  {GreaterEqual[x_, y_] -> y - x, 
              Greater[x_, y_] -> y - x, 
              Equal[x_, y_] -> Abs[x - y], 
              LessEqual[x_, y_] -> x - y, Less[x_, y_]  -> x - y}, 
          r2  =  {And  ->  Max,  
              Or  ->  Min}}, (# ≤  
                0) &  /@  ((cons  /.  r1) /.  
              r2)]

For example, minimize (x^(2 ))subject to x>= 3 or x <= -2

cnvrt[{Or[x ≥ 3, x ≤ -2]}]

callLGO[
      (x^2), {Min[3 - x, 2 + x] ≤  0}, {{x, -5, 
          5}}]

Find the shortest distance to either of two circles with radius = 1, centered at {1,1} and {-1,-1} to a point at {-.5, +.6}

cnvrt[{(((x - 1))^2) + (((y - 1))^2) ≤  
            1  ||  (((x + 1))^2) + (((y + 1))^2) ≤  
            1}]

sln = 
      callLGO[(((x +  .5))^2) + (((y -  .6))^2), {Min[-1 + (((-1 + x))^2) + (((-1 + y))^2), -1 + (((1 + x))^2) + (((1 + y))^2)] ≤ 0}, {{x, -2, 2}, {y, -2, 2}}]

The point is on the circle centered at { 1, 1 }, since its distance to that point is 1.

{(((x - 1))^2) + (((y - 1))^2) ,  (((x + 1))^2) + (((y + 1))^2) } /.  sln[[2]]

The distance to the other circle is larger.

callLGO[(((x +  .5))^2) + (((y -  .6))^2), {(((x + 1))^2) + (((y + 1))^2) ≤  1}, {{x, -2, 2}, {y, -2, 2}}]

The solution is shown in the following plot.

MathOptimizer Professional: Some Recent Applications

The Zimmermann (RecMath) Circle Packing Contest

See http://www.recmath.org/contest/CirclePacking/index.php

Pack N non-overlapping circles with radii from 1 to N into as small a circle as possible, with N having values from 5 to 50.

The scoring process used was very nonlinear:
Q=(best submitted radius)/(your submitted radius)
Score=(0.7*Q)+(0.1*Q^16)+(0.1*Q^128)+(0.1*Q^1024)

One of us (FK) came in 47th out of 155 competitors,using MOP with a score of 35.6776 of a possible 46, having submitted 62 solutions. Several other MOP users also competed. The best score by the team of Addis, Locatelli and Schoen was 45.1379 with 765 submissions. As a side note, AL&S are our good colleagues, and personal communication between Schoen and JP revealed that they spent several months(!) on solving these problems. By contrast, we simply used MOP in "blind" illustrative runs, without further model structure studies.

The following graph shows Q as a function of number of circles for FK's solutions.Typically the MOP result was within several percent of the best solution found.

After the contest, methods were revealed and we found out that a good method for improving solutions is to repeatedly swap circles adjacent in size.  This was implemented in MOP, using local search and the option Recompile -> False, since the only change for each search was the initial values of the unknowns (Castillo, Kampas, and Pintér, 2006).

As an example, consider the following plot of the packing of 10 circles of radii Range[10]^(-1/2).  The circumscribing circle has a radius of 1.94151.

By repeatedly swapping circles of adjacent size, the radius of the enclosing circle is reduced to 1.9266.

The following graph shows results for repeated swapping of adjacent spheres for packing of 10 spheres of radii Range[10]^(-1/3).  The technique appears to work even better for spheres than circles.

Optimization of Gas Production at Shell IEP

The following brief notes are based on our joint modeling and numerical study with Shell's researchers as described in (Mason, Bagirov, Emelle, Kampas, Pintér and van Berkel, 2006).

We analyze and then optimize the blending of gas produced from 5 fields ((Subscript[F, 1]).. (Subscript[F, 5])) to supply different quality gas at six processing plans ((Subscript[P, 1]).. (Subscript[P, 6])) through a converging-diverging gas gathering/distribution network.

The model considered the reservoir pressures, the gas inflow performance relationship, the gas vertical flow performance, the well pressure-rate performance, the pressure drop along horizontal pipes and the network pressure-rate equations.  The goal of the calculation was the maximization of LPG content delivered to plant 6 while meeting carbon dioxide and nitrogen concentration limits. Using the LGO solver (essentially the same solver that is built into MOP) has resulted in significant production improvement when compared to results with Shell's previously applied optimizaiton methodology.

Global Optimization: Application Perspectives

Global optimization is needed and brings tangible benefits to an increasing range of applications. These include e.g. the following areas (in many cases based on actual LGO or MOP software client/user needs):

• agricultural modeling 
• air traffic control 
• aircraft trajectory optimization 
• artificial intelligence 
• chemical and phase equilibrium problems
• chemical and process engineering 
• computational chemistry 
• computer experiments 
• constraint satisfaction 
• object cutting and packing 
• crystallography 
• environmental sciences 
• experimental design 
• facility location 
• fuzzy control 
• finite element modeling combined with (global) optimization
• inverse problems in hydrology and geophysics 
• layout design 
• optimized lighting of a given area or volume 
• machine learning 
• material design 
• mechanical engineering design 
• model fitting do data (calibration) 
• molecular structure modeling 
• network flows 
• neural networks (training) 
• nonlinear equations and inequalities
• nonlinear regression 
• process synthesis 
• protein folding studies 
• quadratic assignment problems 
• quantum mechanics 
• radiotherapy 
• robotics 
• structural system identification 
• VLSI circuit design 
• wastewater engineering

Please consult the illustrative references listed in the next section, with further pointers therein. For instance, the edited volume (Pintér, 2006) presents a collection of nearly two dozen engineering and scientific case studies.

References

Castillo, I., Kampas, F.J., and Pintér, J.D. (2006) Solving circle packing problems by global optimization: numerical results and industrial applications. (Submitted for publication.)

Kampas, F.J. and Pintér, J.D. (2004) Generalized circle packings: model formulations and numerical results. Proceedings of the 2004 International Mathematica Symposium, Banff, AB.

Kampas, F.J. and Pintér, J.D. (2006) Configuration analysis and design by using optimization tools in Mathematica. The Mathematica Journal 10 (1), 128-154.

Mason, T.L., Bagirov, A., Emelle, C., Kampas, F.J., Pintér, J.D., and van Berkel, J. (2006) Application of Global Optimization Methods for Production Optimization. Presented at the 21st European Conference on Operations Research (Reykjavik, Iceland, July 2006). Paper to appear in: Journal of Industrial and Management Optimization.

Pintér, J.D. (1996) Global Optimization in Action. Kluwer Academic Publishers, Dordrecht.

Pintér, J.D. (2002) Global optimization: software, test problems, and applications. In: Pardalos P.M. and Romeijn, H.E. (Editors), Handbook of Global Optimization, Vol. 2, pp. 515-569. Kluwer Academic Publishers, Dordrecht.

Pintér, J.D. (2005) LGO – An Integrated Model Development and Solver Environment for Continuous Global Optimization. User Guide. Pintér Consulting Services, Inc., Halifax, NS, Canada. (Published with updates since 1996.)

Pintér, J.D. and Kampas, F.J. (2003) MathOptimizer Professional – An Advanced Modeling and Optimization System for Mathematica Users with an External Solver Link. User Guide. Pintér Consulting Services, Inc., Halifax, NS, Canada.

Pintér, J.D. and Kampas, F.J. (2005) Global optimization in Mathematica with MathOptimizer Professional. Mathematica in Education and Research 10 (2), 1-18.

Pintér, J.D. and Kampas, F.J. (2006) MathOptimizer Professional: key features and illustrative applications. In: Liberti, L., and Maculan, N. (Editors), Global Optimization: From Theory to Implementation, pp. 263-279. Springer Science + Business Media, New York.

Pintér, J.D. (Editor) (2006) Global Optimization – Scientific and Engineering Case Studies. Springer Science + Business Media, New York, 2006.

Trefethen, N. (2002) A Hundred-Dollar, Hundred-Digit Challenge.SIAM News 35, No. 1, Jan./Feb. 2002 issue. http://www.siam.org/siamnews/01-02/challenge.pdf.

Weisstein, E.W. (2006) Hundred-Dollar, Hundred-Digit Challenge Problems. From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Hundred-DollarHundred-DigitChallengeProblems.html.

Wolfram, S. (2003) The Mathematica Book. (5th Edition.) Wolfram Media, Inc., Champaign, IL.

Wolfram Research (2006) Mathematica. http://www.wolfram.com/.