Getting Started with MathOptimizer Professional

Software development and documentation by János D. Pintér and Frank J. Kampas

Summary

The MathOptimizer Professional software package combines Mathematica’s optimization model development capabilities with the external Lipschitz Global Optimizer (LGO) solver suite.

MathOptimizer Professional solves (numerically) nonlinear optimization models, concisely stated as follows:

Minimize f(x) under the general constraints g(x)≤0, and the box constraints xl≤x≤xu.

Here x is the real n-vector of decision variables (a point of the n-dimensional vector space (R^n));

f(x): is the objective function (f:(R^n)->R);

g(x) is the m-vector of constraint functions (g:(R^n)->(R^m));

[xl, xu] defines an embedding n-interval (box) for the possible (initial) settings of x.

The suite of algorithms integrated by LGO is capable of handling the general nonlinear optimization model-class defined above under very general conditions. Essentially, only the continuity or Lipschitz-continuity of the optimization model functions f and g is required over the – by assumption, finite – box [xl, xu].

LGO – and hence, MathOptimizer Professional – can be used in global and local optimization solver modes; it operates making use of direct model function evaluations, without requiring higher-order (gradient or Hessian) information.

MathOptimizer Professional automatically transforms optimization models formulated in Mathematica into C or Fortran code, then the coded model is handed over to LGO for solution. Following a seamless optimization model compilation / linking / execution procedure, the optimization results are directly returned to the calling Mathematica (notebook) document. This approach can lead to significant program execution speedup which becomes more noticeable with increasing model sizes. MathOptimizer Professional can be used to handle models with thousands of decision variables and general constraints (as of 2015); algorithmic advances and more powerful computers will enable the handling of even larger models.

In the Mathematica document titled Getting Started with MathOptimizer Professional we review the key features of MathOptimizer Professional using a set of illustrative examples. We also provide a selection of references to some of our topical publications.

Please note that – due to the continuing development of MathOptimizer Professional and to the use of various compilers tested with the product – there could be small differences between the results cited here and the results obtained by your actual MathOptimizer Professional product copy.

Further information and examples will be provided to licensed users upon request.

Mathematica Platform and Date

$Version

DateList[]

Installing and Activating (Launching) MathOptimizer Professional

Installation

First, copy the MathOptimizer Professional installation files into the following directory:

Your-Mathematica-installation-directory\AddOns\Applications\MathOptimizerPro**

The symbols ** shown above will correspond to the compiler used in conjunction with MathOptimizer Professional: this aspect will be discussed below under Activation.

An example of a MathOptimizer Professional installation directory is C:\Mathematica\10\AddOns\Applications\MathOptimizerProLF.

Next, copy the Format.m file provided into Your-Mathematica-installation-directory\AddOns\Applications\ directory. This file is used by all MathOptimizer Professional product versions.

Activation

MathOptimizer Professional operations require access to a C or Fortran compiler, in order to generate intermediate files used by the LGO solver engine.

We assume that most of our users will have no problem to install a compiler of their choice, including the required environment variable settings. We will assist licensed users of MathOptimizer Professional with this task if needed/ We can also make available MathOptimizer Professional deliveries for other compilers than the ones listed below, upon request from licensed users.

Please note that direct C or Fortran code development is not required, since all such code is generated automatically by MathOptimizer Professional.

The currently supported (tested) compiler options include the following:

GFortran (GF)

GNU Compiler Collection (GCC)

Intel Fortran (IF)

Lahey Fortran (LF)

Microsoft Visual C/C++ (VC)

Salford Fortran (SF).

To reflect these options, there are small (secondary) setup differences specific to the compiler used in the background. None of these internal differences have any effect on the MathOptimizerPro user's work, or regarding MathOptimizer Professional functionality. However, the compiler used can affect (to some extent) the precision of the results as well as the program execution speed.

Please activate MathOptimizer Professional using exactly one of the statements below, depending on the compiler that you are using.

GCC

Needs["MathOptimizerProGCC`callLGO`"];

GFortran

Needs["MathOptimizerProGF`callLGO`"];

Intel Visual Fortran

Needs["MathOptimizerProIF`callLGO`"];

Lahey Fortran

Needs["MathOptimizerProLF`callLGO`"];

Microsoft Visual C (version 6 and higher)

Needs["MathOptimizerProVC`callLGO`"];

Salford Fortran

Needs["MathOptimizerProSF`callLGO`"];

Minimize, but do not close the system call “DOS box” if such a window appears (this is the case for earlier MathOptimizer Professional versions). The “DOS box” window will display compiler operations during the MathOptimizer Professional session. In current MathOptimizer Professional versions this “DOS box” is not used: hence, users receive information (only) from the MathOptimizer Professional output. More details can be obtained from the optionally generated LGO output files.

MathOptimizer Professional Definition and Options, with an Introductory Example

The basic definition and usage of the callLGO function can be retrieved by the following statement:

?callLGO

The callLGO solver options and their default settings can be queried by the following statement:

Options[callLGO]

To illustrate the basic usage of callLGO, we start with a simple optimization model example.

We want to minimize x^2+(y^2) under the (general) constraints {x+y-1==0,(x^2)+3*y≤2}, and the box constraints {{x,-2,3},{y,-3,2}}; we will also assume an initial solution guess: x=1 and y=-1. Note that providing such a guess is not required by LGO: in its absence the midpoint of the n-interval [xl, xu] is used.

callLGO[(x^2)+(y^2),{x+y-1==0,(x^2)+3*y≤2},{{x,-2,1,3},{y,-3,-1,2}}, SplitEqualityConstraints->True]

Notice that since the above model is convex, it suffices to use LGO only in its local search mode.

callLGO[(x^2)+(y^2),{x+y-1==0,(x^2)+3*y≤2},{{x,-2,1,3},{y,-3,-1,2}},LGOMethod-> "LocalSearch",SplitEqualityConstraints->True]

If you are interested, then now you can take a look at the automatically generated model function file LGOFCT.C or LGOFCT.FOR (depending on the C or Fortran compiler used) and the LGO input parameter file LGO.IN. The auto-generated output files are called LGO.OUT and LGO.SUM. All these files can be found in the MathOptimizer Professional installation directory. This information can be helpful in model development and debugging sessions.

The input files and the summary result (LGO.SUM) generated can be displayed directly in the calling notebook, before Mathematica's reply to callLGO, see the next statement. In some Mathematica implementations, the added delimiter symbols make it more difficult to read these files in the Mathematica environment, but their inspection by using a simple text editor – such as Notepad – is always possible.

result=callLGO[(x^2)+(y^2),{x+y-1==0,(x^2)+3*y≤2},{{x,-2,1,3},{y,-3,-1,2}},LGOMethod-> "LocalSearch", SplitEqualityConstraints->True,ShowLGOInputs->True]

A Simple Model Development Example: Box Design

First, we define the lists of decision variables, (optional) nominal values, lower and upper bounds:

variables={d, h, w}; 
varnom=Table[50., {3}]; 
varlb=Table[1., {3}]; 
varub=Table[300., {3}];

The next statement generates the list of variables, bounds, and the nominal value as required by the MathOptimizer Professional input format:

varswithbounds=Transpose[{variables,varlb,varnom,varub}];

Next, we define the list of model constraints (that may be empty in some models):

constraints={-2*(h*d+h*w+d*w)+888≤0,-h*d*w+1512≤0,h-0.718*w≤0,-h+0.518*w≤0,d*w-252≤0};

and the model objective function:

objective=0.1*(d*w+d*h)+0.2*w*h;

Observe that the model is non-convex, due to its constraints set. Therefore we shall apply the LGO solver in its (default) global search mode.
The next statement invokes LGO to solve the model.

callLGO[objective,constraints,varswithbounds]

This rather easy model can be solved also by applying only local search.

?LGOMethod

callLGO[objective,constraints,varswithbounds,LGOMethod->"LocalSearch"]

For additional information, next we display also the summary LGO result file LGO.SUM, as well as a more detailed runtime output file LGO.OUT. Please inspect these files (located in your MathOptimizer Professional installation directory) which summarize the operations of the underlying LGO solver suite.

callLGO[objective, constraints, varswithbounds, PrintLevel->2]

An Illustrative Application Model

The next small-scale model formulation is based on an engineering problem.

(This model has been sent to us for tests by a client, without background information or explanation. The solution obtained as shown below was confirmed to be as globally optimal by the client.)

obj=-0.5((Sqrt[qH])+(Sqrt[qL])-tH-tL);
cons = {tL-qL*βL≥0,tH-qH*βH≥0,tL-qL*βL≥ tH-qH*βL,tH-qH*βH≥ tL-qL*βH} /. { βH -> 2, βL -> 20};
varswithbounds = {{qL,0,1},{qH,0,1},{tH,0,1},{tL,0,1}};

The next sequence of results shows that even a small-scale global search effort leads to close-to-optimal solutions: see the next input line and its results.

callLGO[obj,cons,varswithbounds,DifficultyFactor->0.125]
callLGO[obj,cons,varswithbounds,DifficultyFactor->0.25]
callLGO[obj,cons,varswithbounds,DifficultyFactor->0.5]
callLGO[obj,cons,varswithbounds]
callLGO[obj,cons,varswithbounds,DifficultyFactor->2]
callLGO[obj,cons,varswithbounds,DifficultyFactor->5]
callLGO[obj,cons,varswithbounds,DifficultyFactor->10]

The next result shows that local search on its own leads to a suboptimal solution (which may not be acceptable also due to the fact that the client preferred positive solution components).

callLGO[obj,cons,varswithbounds,LGOMethod->"LocalSearch"]

A Difficult One-Dimensional Test Problem

Models that involve highly multimodal functions can be very difficult, even in low dimensions as illustrated below by a simple box-constrained example.

objective = 0.15* x+Sin[5 x] Cos[3 x]-Sin[3 x] Cos[2. x]+Sin[x] Cos[7 x];
constraints = {};
bounds = {{x,-3,10}};
result=callLGO[objective,constraints,bounds]
Plot[objective,{x,-3,10}]

In this case we can visually ‘verify’ (by inspecting the figure shown above) that the global solution has been found indeed. In higher dimensions (n>1), such verification is typically not possible.

Optimal objective function value:

result[[1]]

Solution (components):

result[[2]]

Maximal constraint violation by the solution found:

result[[3]]//InputForm

If we add a general constraint(s) then the problem could become (even) more difficult.

constraints1 = {(x^2)-x≤1};
result1=callLGO[objective,constraints1,bounds]

constraints2 = {(x^2)-x-1≤0,4 Sin[3 x]+ 2 Cos[2 x]+5.6≤ 0};
result2=callLGO[objective,constraints2,bounds]

Extending a One-Dimensional Test Problem to Two Dimensions: Curse of Dimensionality

Next we directly extend the model defined in the previous section to (just) two dimensions, to illustrate the phenomenon called ‘the curse of dimensionality’.

objective2 = 0.15* x+Sin[5 x] Cos[3 x]-Sin[3 x] Cos[2 x]+Sin[x] Cos[7 x]+0.15* y+Sin[5 y] Cos[3 y]-Sin[3 y] Cos[2 y]+Sin[y] Cos[7 y];
constraints2 = {};
bounds2 = {{x,-3,10},{y,-3,10}};
result=callLGO[objective2,constraints2,bounds2, DifficultyFactor -> 10]
Plot3D[objective2,{x,-3,10},{y,-3,10}]

Notice that this model has hundreds of local optima compared to the one-dimensional box-constrained model instance that has less than 20 local solutions.

Due to the special form of the model extension chosen, we can also directly see that the solution received is the global optimum indeed.

Optimal objective function value:

result[[1]]

Solution (components):

result[[2]]

Maximal constraint violation by the solution found:

result[[3]]//InputForm

A Small, but Not-so-simple Constrained Optimization Problem, with Added Model Visualization Features

It is straightforward to define an optimization model and then solve it in a single statement. Let us recall the basic calling statement structure as shown below:

callLGO[
objective function, 
{constraint list}, 
{list of {variable, lower bound, nominal value(optional), upper bound} lists}]

This usage is illustrated by the following example.

result=callLGO[(x^2)+(y^2)+10 Sin[x^3+x y],{(x^2)+y^3≤ 1,x*y≤10},{{x,-2,1,3},{y,-3,2,2}}]

Optimal objective function value:

result[[1]]

Solution components:

result[[2]]

Maximal constraint violation by the solution found:

result[[3]]//InputForm

In such a small 2-variable problem we are able to visualize the optimization problem.

Here is the 3d-plot of the objective function:

Plot3D[x^2+y^2+10 Sin[x^3+x y],{x,-2,3},{y,-3,2}]

Here is the contour plot of the same function:

cp=ContourPlot[(x^2)+(y^2)+10 Sin[(x^3)+x y],{x,-2,3},{y,-3,2}]

And here is the plot of the feasible region:

rp=RegionPlot[(x^2)+(y^3)≤ 1 && x*y≤10 ,{x,-2,3},{y,-3,2}]

sp=Graphics[{PointSize[0.05],Point[{-0.9641861921302124,0.4128044822848619}]}];

Show[{cp,rp,sp}]

These plots jointly indicate that the solution found is the (constrained numerical) global optimum.

A Scalable Global Optimization Model Class with Randomized Solutions

In the next example first we generate a random point that will become the solution of an induced global optimization problem. (Hence, every optimization run will lead to a different model and corresponding solution.)

xr=-5+10 Random[]
objective=(x-xr)^2+10.*Sin[5*(x-xr)]^2
bounds={x,-15,5}
callLGO[objective,{},{bounds}]
Plot[objective,{x,-5,5}]

Let us point out that changing the model parameters can lead to models of various difficulty. It is entirely possible to generate very tough challenges this way. Therefore it can happen that MathOptimizer Professional – or other global solvers – will be unable to find the true global solution. However, it can be duly expected that a global solver should find reasonably high quality feasible numerical solutions.

Finding Pythagorean Triplets by Global Optimization

The integrality of the solution is enforced by a set of ‘tricky’ constraints.

callLGO[Abs[(x^2)+(y^2)-(z^2)],{Sin[π x]==0,Sin[π y]==0,Sin[π z]==0},{{x,1,10},{y,1,10},{z,1,10}}]

Further Small-scale, but Non-trivial Illustrative Models

Even small models may be not trivial, as the next example shows (the objective and constraint functions are all non-convex).

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

You can modify this example by changing the model statement, then solve the resulting model, as shown by the example below.

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

Trefethen's Problem 4

The following summary description is cited from http://mathworld.wolfram.com/Hundred-DollarHundred-DigitChallengeProblems.html, written by Eric W. Weisstein.

“The hundred-dollar, hundred-digits challenge problems are a set of ten problems in numerical analysis published in the January/February 2002 issue of SIAM News (http://www.siam.org/siamnews/01-02/challenge.pdf). Nick Trefethen [L. N. Trefethen, FRS, Professor of Numerical Analysis, University of Oxford], the proposer of the problems, offered a $100 prize to the person or group obtaining the largest number of correct digits (up to a maximum of 10) for these problems by May 20, 2002. Trefethen underestimated the ingenuity of problem solvers, and 20 independent teams obtained 10 correct digits for all 10 problems. After an anonymous donor stepped in to help defray the larger than anticipated payoffs, prize checks went out to all winners by December 2002.”

Here we discuss only one of the HDHD problems which is a small-scale yet rather difficult numerical global optimization challenge with many local optima as shown by the surface and contour plots the objective function defined below. (We solved this problem using our own optimization software versions in 2002 before the solution was announced, and have been using it as one of our standard test problems ever since.)

objective=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);

Plot3D[objective,{x,-3,3},{y,-3,3}]

ContourPlot[objective,{x,-3,3},{y,-3,3}]

We try first local search: as it can be expected, (in general) a local optimizaiton approach will be inadequate for such a highly multi-extremal problem.

callLGO[objective,{},{{x,-3,0,3},{y,-3,0,3}},LGOMethod->"LocalSearch"]

In default solver mode, we get only an approximate solution: notice however that this approximate solution is significantly better than the one obtained by local search.

callLGO[objective,{},{{x,-3,3},{y,-3,3}}]

The model is solved successfully, when the default solver effort is increased. Notice that the number of search steps (direct model function evaluations, without requiring higher-order information) remains rather small in all cases shown below.

callLGO[objective,{},{{x,-3,3},{y,-3,3}},DifficultyFactor->5]

callLGO[objective,{},{{x,-3,3},{y,-3,3}},DifficultyFactor->10]

callLGO[objective,{},{{x,-3,3},{y,-3,3}},DifficultyFactor->20]

The solutions in the last two cases are almost identical. A more detailed direct printout of the results received shows that the last result is somewhat better:

-3.306865164664285   (obtained using DifficultyFactor->10)

-3.306868647475238   (obtained using DifficultyFactor->20)

Let us note that the second numerical solution displayed above agrees with the announced solution -3.306868647 to 10-digit precision.

For further details and background information, cf. e.g. http://mathworld.wolfram.com/news/2002-05-25/challenge/ by Eric W. Weisstein.

A Process Industry Model

The next model is based on process industry design problem which is often used in tests. This model is rather poorly scaled and several solution components lie on the boundary of the feasible region (a situation to avoid if at all possible – but without making the model irrealistic, of course).

Due to these reasons, the model shown below can give difficulty to nonlinear solvers, including MathOptimizer Professional.

We illustrate the use of several modeling options that are aimed at obtaining better quality solutions (if such solutions can be expected).

Notice that several variables could be eliminated, but we don't discuss this option here (also to follow the original model).

obj=-9x5-15  x8+6x1+16x2+10(x6+x7);
cons={x1+x2-x3-x4==0,0.03x1+0.01x2-x9(x3+x4)==0,x3+x6-x5==0,x4+x7-x8==0,
x9*x3+0.02x6-0.025x5≤0,x9*x4+0.02x7-0.015x8≤0};
varswithbounds={{x1,0,300},{x2,0,300},{x3,0,100},{x4,0,200},{x5,0,100},{x6,0,300},{x7,0,100},{x8,0,200},{x9,0,.03}};
callLGO[obj,cons,varswithbounds]

Constraint scaling (dividing each constraint by the largest multiplier when applicable).

obj=-9x5-15  x8+6x1+16x2+10(x6+x7);
cons1={x1+x2-x3-x4==0,(0.03x1+0.01x2-x9(x3+x4))/0.03==0,x3+x6-x5==0,x4+x7-x8==0,
(x9*x3+0.02x6-0.025x5)/0.025≤0,(x9*x4+0.02x7-0.015x8)/0.02≤0};
varswithbounds={{x1,0,300},{x2,0,300},{x3,0,100},{x4,0,200},{x5,0,100},{x6,0,300},{x7,0,100},{x8,0,200},{x9,0,.03}};
(* varswithbounds={{x1,0,10},{x2,0,200},{x3,0,10},{x4,0,200},{x5,0,10},{x6,0,10},{x7,50,100},{x8,180,200},{x9,0,.03}}; *)
callLGO[obj,cons1,varswithbounds]

The SplitEqualityConstraints option replaces the equality constraints by a pair of inequalites: this simple ‘trick’ can be helpful – together with constraint scaling – as it is illustrated below.

callLGO[obj,cons1,varswithbounds,SplitEqualityConstraints->True]

One can always use ‘brute force’ by increasing the value of DifficultyFactor. Here it does not lead to improved results since -400 is the global optimum indeed. The model has many suboptimal solutions.

callLGO[obj,cons1,varswithbounds,SplitEqualityConstraints->True,DifficultyFactor->10]

A Scalable Convex Quadratic Test Function Class, Solved in Local Solver Mode

Test setup 1
Fixed variable bounds, well-scaled objective, one added constraint. Optimal solution x[i]=1/n, for all i; n=1,2,3,... Objective function at optimal solution = 1/n.

vars[n_]:=Array[x,n]; 
vwb[n_]:=Thread[{vars[n],-n RandomReal[],n RandomReal[]}]; 
obfunc[n_]:=(∑_(i=1)^n)((x[i])^2);
cns[n_]:=(∑_(i=1)^n)x[i]==1;
setup[n_]:={obfunc[n],{cns[n]},vwb[n]};

AbsoluteTiming[s10=callLGO[Sequence @@ setup[10], LGOMethod -> "LocalSearch"]]

AbsoluteTiming[s100=callLGO[Sequence @@ setup[100], LGOMethod -> "LocalSearch"]]

Test setup 2
Random variable bounds, increasingly poorly scaled objective, no added constraints.  Optimal solution x[i]=0, for all i; n=1,2,3,...Objective function at optimal solution = 0.

vars[n_]:=Array[x,n]
vwb[n_]:=Thread[{vars[n],-1,RandomReal[],1}]
obfunc[n_]:=(∑_(i=1)^n)i ((x[i])^2)
setup[n_]:={obfunc[n],{},vwb[n]}

10-variable model instance.

AbsoluteTiming[s10=callLGO[Sequence @@ setup[10], LGOMethod -> "LocalSearch"]]

100-variable model instance.

AbsoluteTiming[s100=callLGO[Sequence @@ setup[100],LGOMethod -> "LocalSearch" ]]

1000-variable model instance.

AbsoluteTiming[s1000=callLGO[Sequence @@ setup[1000], LGOMethod -> "LocalSearch"];]

Objective Function

s1000[[1]]

2000-variable model instance

AbsoluteTiming[s2000=callLGO[Sequence @@ setup[2000], LGOMethod -> "LocalSearch"];]

s2000[[1]]

A Scalable Non-Convex Quadratic Test Function Class, Solved in Local and Global Solver Modes

Fixed variable bounds, well-scaled objective, no added constraints. Optimal solution x[i]=0 for even i, x[i]=-1 for odd i, for all i=1,...,n ; for n=1,2,3,...

vars[n_]:=Array[x,n]
vwb[n_]:=Thread[{vars[n],-1,0.789123456}]
obfunc[n_]:=(∑_(i=1)^n)(-1)^i*((x[i])^2)
setup[n_]:={obfunc[n],{},vwb[n]}

AbsoluteTiming[s10=callLGO[Sequence @@ setup[10], LGOMethod -> "LocalSearch"]]

Next, all three global solver modes are used: note that their solution overrules the initial local solver result, thus (on purpose) all solver models may give numerically different results. This can be rectified by increasing the global search effort thereby "approximating the convergence theory".

AbsoluteTiming[s10=callLGO[Sequence @@ setup[10], LGOMethod -> "BranchAndBound"]]

AbsoluteTiming[s10=callLGO[Sequence @@ setup[10], LGOMethod -> "BranchAndBound", DifficultyFactor -> 10]]

AbsoluteTiming[s10=callLGO[Sequence @@ setup[10], LGOMethod -> "RandomSearch"]]

AbsoluteTiming[s10=callLGO[Sequence @@ setup[10], LGOMethod -> "RandomSearch", DifficultyFactor -> 10]]

AbsoluteTiming[s10=callLGO[Sequence @@ setup[10], LGOMethod -> "MultiStart"]]

AbsoluteTiming[s100=callLGO[Sequence @@ setup[100], LGOMethod -> "LocalSearch"]]

AbsoluteTiming[s1000=callLGO[Sequence @@ setup[1000], LGOMethod -> "LocalSearch"]][[1]]

s1000[[1]]

A Scalable Optimization Problem Class Defined by Using Hilbert Matrices

Here we define and solve a quadratic optimization model with an embedded scalable Hilbert marix.

The theoretical optimum value is zero, at the zero vector, in all model instances.

As the model size grows, this model becomes more and more badly scaled: this can lead to numerical inaccuracies.

vars[n_]:= Table[x[i],{i,n}]
h[n_]:=Table[(1/(i+j-1)),{i,n},{j,n}]
Clear[varswithbounds];
varswithbounds[n_]:=Table[{x[i],-3,0,5},{i,n}]
obj[n_]:=vars[n].h[n].vars[n]

callLGO[obj[5],{},varswithbounds[5]]

callLGO[obj[20],{},varswithbounds[20]]

callLGO[obj[30],{},varswithbounds[30]]

Packing Non-Uniform Size Circles into an Optimized Circle

Given a set of (arbitrary size) circles, our objective is to find a ‘container’ circle with minimal radius that can hold a pairwise non-overlapping arrangement of the given set of circles.

This is a well-known and notoriously difficult instance of the broad class of object packing problems (which includes also far more difficult problems).

Model setup for the generic n-circle problem-instance.

Clear[obj,cons,vars]

circumscribecircles[n_Integer,func_,opts___]:=
Module[
{cons,cons1,cons2,vars,lb,ub,soln,rcirc,pd,obj},
Do[r[i]=func[i],{i,n}];
cons1=Flatten [Table[r[i]+r[j] ≤  (Sqrt[(((x[i]-x[j]))^2)+(((y[i]-y[j]))^2)]),{i,n-1},{j,i+1,n}]];
cons2 = Table[obj -r[i] ≥(Sqrt[ ((x[i])^2)+((y[i])^2)]),{i,n}];
vars=Join[Flatten[Table[{x[i],y[i]},{i,n}]],{obj}];
cons = Join[cons1,cons2];
lb=Join[Table[-2.5,{2*n}],{1}];
ub=Join[Table[2.5,{2*n}],{2.5}];
vwb=Transpose[{vars,lb,ub}];
callLGO[obj,cons,vwb,opts]
]

plt[soln_,n_Integer,func_, rmax_]:=
Module[
{radii,centers,t},
radii = func /@ Range[n];
radii = Join[{soln[[1]]},radii];
centers = Table[{x[i],y[i]},{i,n}]/. soln[[2]];
centers = Join[{{0,0}},centers];
Graphics[#, ImageSize -> Small,PlotRange -> {{-rmax,rmax},{-rmax,rmax}}]&  @ Table[Circle[centers[[i]],radii[[i]]],{i,n+1}]
]

Numerical solution of a 3-circle instance in which the radii are defined as r[i]=(i^(-1/2)). The generated model is also presented.

AbsoluteTiming @ (sln3=circumscribecircles[3,((#)^(-1/2))&,ShowLGOInputs -> True])

plt[sln3,3, ((#)^((-1)/2))&, 2.0]

Numerical solution of a 10-circle instance in which the radii are defined as r[i]=(i^(-1/2)). Notice the (possible, but not necessarily guaranteed) improvement of results with increasing search effort.

AbsoluteTiming @ (sln10=circumscribecircles[10,((#)^(-1/2))&])

plt[sln10,10, ((#)^((-1)/2))&, 2.0]

AbsoluteTiming @ (sln10a=circumscribecircles[10,((#)^(-1/2))&, DifficultyFactor -> 100])

plt[sln10a,10, ((#)^((-1)/2))&, 2.0]

AbsoluteTiming @ (sln10b=circumscribecircles[10,((#)^(-1/2))&, DifficultyFactor -> 300])

plt[sln10b,10, ((#)^((-1)/2))&, 2.0]

Numerical solution of a 20-circle instance in which the radii are defined as r[i]=(i^(-1/2)).

Here we use only default option settings, and leave further experimentation to the Reader.

AbsoluteTiming @ (sln20=circumscribecircles[20,((#)^(-1/2))&])

plt[sln20,20, ((#)^((-1)/2))&, 2.2]

Illustrative References

There are hundreds of textbooks, thousands of research articles, and an ever increasing number of web sites devoted to nonlinear (global and local) optimization. Here we only provide a selection our (more recent) topical research publications: further information is available upon request.

Pintér, J.D. (1996) Global Optimization in Action (Continuous and Lipschitz Optimization: Algorithms, Implementations and Applications). Kluwer Academic Publishers, Dordrecht. Now distributed by Springer Science +  Business Media, New York.

Pintér, J.D. (1999...2015) LGO – A Model Development and Solver System for Global-Local Nonlinear Optimization. User’s Guide. Published by Pintér Consulting Services, Inc., Halifax, NS, Canada. www.pinterconsulting.com.

Pintér, J.D. (2002) Global optimization: software, test problems, and applications. In: Pardalos, P. M. and Romeijn, H. E., eds. Handbook of Global Optimization, Volume 2, pp. 515-569. Kluwer Academic Publishers, Dordrecht, 2002.  Now distributed by Springer Science +  Business Media, New York.

Pintér, J.D. and Kampas, F.J. (2005) O.R. model development and optimization with Mathematica. In: Golden, B., Raghavan, S., and Wasil, E., eds. The Next Wave in Computing, Optimization, and Decision Technologies, pp. 285-302. Springer Science + Business Media, New York.

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

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.

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

Lasdon, L.S. and Pintér, J.D. (2007) Computational global optimization. In: Klastorin, Th., ed. TutORials in Operations Research, Chapter 2. Published by the Institute for Operations Research and the Management Sciences (INFORMS), Hanover, MD.

Castillo, I., Kampas, F.J. and Pintér, J.D. (2008) Solving circle packing problems by global optimization: numerical results and industrial applications. European Journal of Operational Research 191, 786-802.

Pintér, J.D. (2008) Model development and optimization in interactive computing environments. Central European Journal of Operations Research 16, 165-178.

Pintér, J.D. (2009) Software development for global optimization. In: Pardalos, P.M. and T. F. Coleman, eds. Global Optimization: Methods and Applications, pp. 183-204. Fields Institute Communications Volume 55. Published by the American Mathematical Society, Providence, RI.

Pintér, J.D. and Kampas, F.J. (2013) Benchmarking nonlinear optimization software in technical computing environments: Global optimization in Mathematica with MathOptimizer Professional. TOP (An official journal of the Spanish Society of Statistics and Operations Research) 21, 133-162.

Fasano, G. and Pintér, J.D., eds. (2015) Optimized Packings and their Applications. Springer Science + Business Media, New York. (To appear)

Kampas, F.J. and Pintér, J.D. (2015) Globally optimized packings of non-uniform size circles and spheres: A scalable model development and customized algorithmic solution approach, with illustrative results. In: Fasano, G. and Pintér, J.D., eds., Optimized Packings and their Applications. Springer Science + Business Media, New York. (To appear)

Further Information and Technical Support

Licensed MathOptimizer Professional users: please feel free to contact the developers for additional details:

Thanks for your attention: we welcome your suggestions regarding the software and its documentation.

János D. Pintér, PhD, DSc                         	Frank J. Kampas, PhD, MBA
Pintér Consulting Services, Inc.         		Physicist at Large Consulting, LLC
Email: janos.d.pinter@gmail.com                	Email: frank@physicistatlarge.com
Website: www.pinterconsulting.com