MathOptimizer

An Advanced Nonlinear  
Optimization Program System 
For Mathematica Users

User's Guide

Frank J. Kampas and János D. Pintér

MathOptimizer Product Development and Support

MathOptimizer is developed and supported by Pintér Consulting Services, Inc.

Principal product developers: Frank J. Kampas, PhD, MBA and János D. Pintér, PhD, DSc.

Support requests and technical questions can be sent via email to fkampas@msn.com and janos.d.pinter@gmail.com.

We also offer topical consulting services, workshops and tutorials, as well as customized software development.

Trademarks, Copyright, and License Agreement

Trademarks and Copyright

Mathematica is a trademark of Wolfram Research, Inc. 
www.wolfram.com

The MathOptimizer software package and its documentation have been continuously developed and maintained since 2002.
MathOptimizer is a trademark of Pintér Consulting Services, Inc. 
www.pinterconsulting.com
MathOptimizer User's Guide Copyright © Pintér Consulting Services, Inc. 
Neither the MathOptimizer software, nor this document or parts thereof may be copied, reproduced, translated, engineered or reduced to any readable or interpretable form, without the express written consent of Pintér Consulting Services, Inc., except as permitted by the Software License Agreement.

Software License Agreement

Software License

The MathOptimizer program system and its entire documentation  (referred to below as the 'Product') is the property of Pintér Consulting Services, Inc. (abbreviated below as 'PCS').

PCS grants you, an individual or an organization (referred to below as 'User'), the right to use the MathOptimizer software that you have received from an authorized dealer, following a proper licensing procedure.

The Product may be stored and used on at most two computers, as long as you have a corresponding valid license and you are its only User.

The licensed User may not make available for use the software to any other party, in any format, without a proper licensing agreement. Multiple simultaneous (or separate) licenses are available from PCS and its authorized partners. Please contact PCS and its partners for further licensing details.

Limited Warranty

The MathOptimizer software is provided 'as is', without warranty of any kind, regarding its usability for a specific purpose. The entire risk as to the results and the performance of the Product is assumed by the User. PCS disclaims all warranties, either express or implied, including but not limited to implied warranties of merchantability, and fitness for a particular purpose, with respect to the entire Product.

In no event will PCS be liable for any business related damages – including but not limited to: loss of profits, business interruption, loss of information, and the like – arising out of the use or inability to use the Product.

PCS guarantees that all shipped Products are free from defects in materials and workmanship under normal use and service for a period of 90 days.

The functionality of the Product is also guaranteed to the specific details and examples, exactly as documented by this User's Guide.

Product Replacement and Returns

PCS offers replacement of the Product, if it arrives in unusable form (due to damage caused during shipment). Only shipping and handling charges apply in such cases of replacement.

Product returns are accepted within 30 days of purchase with a full refund, except shipping and handling charges that apply also in such cases. Return expenses are also to be covered by the User. All returned items have to be in resellable – as if brand new – condition.

A brief explanation regarding the reason for return will be much appreciated: we will make every reasonable effort to keep all MathOptimizer users satisfied.

Product Upgrades

MathOptimizer users will be offered the opportunity to purchase upgrades as these are released, at a full deduction of the price already paid. Shipping and handling charges will typically be payable also in such cases.

Acknowledgements

The MathOptimizer development work significantly benefited from advice, books and software, and a visiting scholar grant (to JDP), all generously provided by Wolfram Research, Inc to us. Several colleagues at WRI have been most helpful in answering our technical questions over the years.

The MathOptimizer (Version 1) software development project was made possible by a grant received by JDP from the National Research Council of Canada, and by a subsequent research contract from DRDC-Atlantic (Dartmouth, NS, Canada). The valuable advice and suggestions received from Dr. Christopher Purcell throughout this original development work are also gratefully acknowledged.

MathOptimizer Summary

MathOptimizer is an integrated global-local optimization software package for solving constrained nonlinear optimization problems. Since MathOptimizer is a native Mathematica package, optimization models can be formulated using a very broad range of Mathematica features and functionality.

We shall consider the following general optimization model form:

minimize f(x)

subject to the constraints

g(x) ≤ 0
xlb ≤ x ≤ xub

Here x is an n-dimensional real vector, f is a scalar objective function, and g is an m-component real-valued vector function. We shall assume (require) that the component-wise bounds xlb and xub are finite, and that f and g (the latter component-wise) are continuous. Let us point out that all well-posed, but seemingly more general finite-dimensional constrained optimization models can be brought to the above canonical form by (elementary or more advanced) transformations.

To address the very general model-type stated above, MathOptimizer can be used in a variety of operational modes, driven by option settings. We first review the technical background of MathOptimizer, including also some modeling tips. This is followed by a detailed description of MathOptimizer's options, extensively illustrated by solving simple to more advanced optimization problems.

One of the great advantages of using Mathematica – and, hence, of MathOptimizer – is that the entire application development process can be documented in an interactive Mathematica notebook, or a set of related notebooks. This approach enables high quality project documentation, and fully supports the actual software development and computational work. Mathematica notebooks – being ordinary text files – are portable across all supported hardware and operating system platforms: this feature enables convenient information exchange within the user community. The same applies, of course, to all notebook documents in which MathOptimizer is used.

We assume that the Reader is familiar with the essential concepts and usage of Mathematica, as well as with the basics of nonlinear optimization. All MathOptimizer options are explained and illustrated in sufficient detail, and we will provide topical references for technical details not discussed here.

Nonlinear Optimization: Technical Background Notes and Model Formulation Tips

MathOptimizer Scope of Application

MathOptimizer is aimed at the solution of a very general class of constrained optimization problems. It provides a combination of advanced global and local search procedures for optimizing the value of a given objective function subject to a finite collection of constraints.

Although in principle MathOptimizer can be used to solve arbitrary optimization models, its target area of applicaiton is continuous nonlinear optimization. This includes traditional local optimization, as well as the far more general case of global optimization (GO). GO problems can be defined e.g. by a multimodal objective function on a simple 'box' region, on a linearly constrained set, or on a possibly far more complicated (non-convex and perhaps even disconnected) feasible solution set. A non-convex feasible region in itself can make the numerical solution of the resulting GO model very difficult, even if the objective function is simple: an example is the solution of a system of nonlinear equations without a stated objective function.

For illustration, consider the box-constrained global optimization model

minimize sin((Subscript[x, 1]) – (Subscript[x, 1])(Subscript[x, 2])) + 0.3 ((Subscript[x, 1])^2) 
–2 ≤ (Subscript[x, 1]) ≤ 3 
–5 ≤ (Subscript[x, 2]) ≤ 3

Although the feasible set is just a two-dimensional interval, the objective function is multimodal as shown by the figure below.

Evidently, in lack of an initial point that is 'sufficiently close' to the (unknown) global optimum, a local scope optimization method in itself may not find the best possible solution.

Optimization Methods in MathOptimizer

MathOptimizer combines several local and global scope algorithmic approaches.

The global solver approach is based on an adaptive stochastic search algorithm. Theoretically, this method is globally convergent with probabilty 1, as the sample size grows (to infinity).

The current local solvers are based on Lagrangian duality theory in convex optimization; they also use Mathematica's FindMinimum function. Theoretically, the methods used are locally convergent for sufficiently smooth models defined by twice continuously differentiable model functions.

We wish to emphasize that although all search approaches implemented in MathOptimizer are theoretically convergent, in numerical practice it could occur that the global or even local optimum of some difficult model is not found directly, applying default option settings and a preset (limited) numerical effort. The same limitation applies to all numerical optimization algorithm implementations, even in much simpler contexts such as linear optimizaiton. Many of the examples discussed later on suggest option settings to overcome this inherent difficulty. In 'easy' models, one can expect that all reasonable option settings will generate the same global optimum estimate (except perhaps small numerical differences). In more difficult models, MathOptimizer option changes can lead to a range of numerical solutions, and the user can eventually select the best one.

Due to the objectives and intended scope of this document, we will not discuss the underlying theory that leads to the algorithms implemented in MathOptimizer. Illustrative general references are Bertsekas (1999), Edgar, Himmelblau and Lasdon (2001), Horst and Pardalos (1995), and Pardalos and Romeijn (2002). More specific topical works are e.g. Pintér (1996, 2001, 2002, 2006), as well as Pintér and Kampas (2005), Kampas and Pintér (2006).

Model Formulation, Input and Output Information

Input Information

x	vector of decision variables

xlb	finite lower bound vector of x

xub	finite upper bound vector of x

xnom	(or xinit) nominal (or initial) value of x

f(x)	scalar objective function

g(x)	vector of constraint functions

Model Formulation

The following standard model form is considered throughout this document

minimize f(x) 	
subject to the constraints
g(x) ≤ 0
xlb ≤ x ≤ xub  (xlb < xub component-wise)

All model functions f, g are assumed to be at least computable over the domain [xlb,xub]. If the functions are continuous, then this guarantees the theoretical convergence of the global search algorithm used. In the context of the built-in local scope search methods, at least local smoothness of the model functions is also postulated.

The explicit consideration of lower and upper variable bounds is necessary, since the existence of the global solution set can not be guaranteed without such bounds.

Output of Results

- the optimal solution vector found (i.e., numerically estimated) by the algorithm
- the numerical optimum value
- the maximal constraint function violation at the optimum estimate

Explicit Bound Settings and Variable Range Scaling

By results from classical analysis, the existence of finite bounds xlb ≤ x ≤ xub is needed to assure that the general optimization model form considered indeed has a finite (globally optimal) solution. One can think of minimizing or maximizing the trivially simple univariate function f(x)=x over the interval (-∞, 0) to see that the closedness of the variable range [xlb, xub] is also essential.

It is always good practice to select 'safe', but not 'exaggerated' bounds. In general, global and local search methods will benefit from well-chosen, reasonable bounds and a good quality initial (nominal) solution vector if such is available. If the bounds are possibly binding, then one should add the corresponding bound constraints also to the explicit list of inequality constraints, at least when also applying local search methods. The global search methods embedded in MathOptimizer will never leave the interval [xlb, xub], but in general they will run slower when larger intervals and/or more model constraints are considered.

It is often beneficial and hence advisable to formulate well-scaled models, to improve algorithmic search efficiency. The interval [xlb,xub] can be simply scaled to the unit 'box' [0,1] of corresponding dimensionality, by the transformation

xi = xlbi + (Subscript[x, s,i])(xubi- xlbi)  i=1,...,n

Here n is the model dimensionality, i denotes the indices of the vector components, and (Subscript[x, s])ϵ[0,1] is the scaled vector. If all variable ranges are at least in the same or close order of magnitude, then we may not need the scaling shown above, but in other cases scaling is advisable.

Model Function Scaling

Let us assume that the model variables are already properly scaled, if this was necessary. Next, we shall briefly address the equally important issue of model function scaling.

To this end, first one can introduce a simple equivalent model transformation: we shall minimize z, and add f(x) – z ≤ 0 to the set of constraints. Next, all constraints will be considered individually. Two cases will be distinguished regarding whether the constraint considered gj(x) ≤ 0 has a non-zero constant term in it, or not. For those constraints that do not have such a constant term, a given (preset) tolerance level t > 0 is introduced: this parameter should correspond to the acceptable constraint violation by the solution. All equality constraints will be replaced by the componentwise pair of inequalities –t ≤ gj(x) ≤ t; a similar but only 'one-sided' transformation is applied to all constraints that do not have a non-zero constant term. After the needed such transformations are completed, each model function has the form c(x) ≤ k. Now constraint c(x) ≤ k can be simply replaced by (c(x) – k)/k ≤ 0, if k is sufficently large, that is at least k ≥ 100*t holds. For those constraints where k < 100*t (this will specifically include all equality relations), one can use instead of the original constraint the approximation (c(x) - k)/(1+k) ≤ 0.

Let us point out that there are also other model scaling options, and that there is no guaranteed recipe regarding the 'best possible' model scaling. In many cases some experimentation with scaling and solver options will help to find good quality numerical solutions even in difficult models.

MathOptimizer  Installation

System Requirements

MathOptimizer will run on all hardware and operating system platforms for which a current, complete Mathematica implementation and technical support is available from Wolfram Research.

MathOptimizer has been originally developed and tested using Mathematica 4.0 (in 2002). Since then it has been consistently updated and tested for all new Mathematica versions: in September 2009 the current version is 7.0.1.

We recommend to use Mathematica version 5.2 and higher. In fact, for 5.2, we need to ship a somewhat earlier MathOptimizer version, since the current MathOptimizer implementation uses some of Mathematica's latest features. (This may be another reason to upgrade at least to Mathematica version 6.0.)

The test results and optional timings reported below were obtained using several 'average capability' personal computers (as of 2009). For example, one of these machines has an Intel Core Duo CPU T5450 2.26 Ghz processor., running under the 32-bit Windows XP SP3 operating system. Using Mathematica version 7.0.1, the machine's overall benchmark score is 1.94 based on the MathematicaMark7 Benchmark Report. For comparison, similar category computers (shown in the benchmark study) score in the range of 1 to 4.7. Most of the examples included in this document typically can be solved at most within a few seconds, unless explicitly noted otherwise.

$Version

Date[]

Install and Launch MathOptimizer

The current MathOptimizer installation file system can be obtained on a CD, via e-mail or some other download method. The package files (possibly in operating system specific executable only versions) are intended for loading into Mathematica using the Needs or Get command. The documentation consists of standard notebook (plain text) files, and it is the same for all operating systems.

The installer should input $UserAddOnsDirectory into Mathematica and put the MathOptimizer directory into the directory that is returned. 
The $UserAddOnsDirectory may be hidden, requiring the installer to check "Show Hidden Files" or something similar when looking for it. 
For example, in Suse Linux, $UserAddOnsDirectory is /home/<UserName>/.Mathematica.

Next, place the MathOptimizer documentation into the ...\MathOptimizer\Documentation\English subdirectory. Additional delivery files (technical notes etc.) may also be present, and then can also be placed in the documentation subdirectory.

Following this setup structure, MathOptimizer can be directly invoked for use by the following statement.

Needs["MathOptimizer`MO`"]

The next statement serves only to verify the availability of MathOptimizer's MO package.

$ContextPath

MO.m contains the Optimize package that integrates all MathOptimizer solvers.

?Optimize

As shown above, the Optimize function requires three key arguments. The first argument is the model objective function, the second is the list of model constraints, and the third is list of decision variables with corresponding bounds. Notice that this structure is a bit different from that of NMinimize, which has the objective function as the first argument of a list containing the objective function and the constraints; furthermore, variable bounds are not required.

A number of options can be selected and added to the input list: we will discuss each of these options later on.

Getting Started: Illustrative Examples

A Box-Constrained Model

If there are no constraints, then the second argument of Optimize is an empy list as shown below.

Optimize[(((x-1))^2)+(((y-2))^2),{},{{x,-5,5},{y,-5,5}}]

Optimize returns a list containing the objective function value at the numerical optimum found, a list of rules giving the optimized values for the decision variables, and the maximum constraint violation at the numerical optimum. (The latter should be 0 here, since there are no constraints.)

Explicit Variable Bounds

The stated variable bounds are strictly enforced by MathOptimizer in its default mode, as shown by the next example. This is another difference with NMinimize, in which variable bounds are not required, and thus the search range is set only when explicitly defined.

Optimize[(((x-6))^2)+(((y-2))^2),{},{{x,-5,5},{y,-5,5}}]

NMinimize[(((x-6))^2)+(((y-2))^2),{{x,-5,5},{y,-5,5}}]

Let us remark here again that – in a mathematically precise setting – one does need variable bounds and continuous model functions to guarantee that the global optimization problem itself is well-posed.  The continuity assumptions can be replaced by somewhat milder semicontinuity conditions that will not be discussed here, since continuity is most typically satisfied in the prospective domain of MathOptimizer applications.

It is also noted that, in general, tighter variable bounds lead to better numerical performance: therefore it is worthwile to give due consideration to all variable bound settings.

Setting Initial Variable Values

It is possible to define initial values for some or all of the variables. If such a value is not given, then the midpoint (component-wise central point) of the bounds is used by default. The initial values of the variables are used as one set of starting values for the local search – quite possibly together with further initial points generated by a preceding global search phase, as it will be discussed later on.

Optimize[(((x-1))^2)+(((y-2))^2),{},{{x,-5,4,5},{y,-5,3,5}}]

Equality and Inequality Constraints

If the constraints are all inequalities, then it is possible (but not guaranteed) that the maximum constraint violation will be 0.

Optimize[(((x-1))^2)+(((y-2))^2),{x+y≤ 3},{{x,-2,5},{y,-1,3}}]

If equality constraints are also present, then one can definitely expect corresponding small constraint violations, especially if the exact solution vector has non-rational number components. To illustrate, in the following example, the optimal values for x and y are -(Sqrt[(1/2)])~ -0.707107 which cannot be represented as a machine number.

Optimize[x+y,{(x^2)+(y^2)≤ 1},{{x,-1,1},{y,-1,1}}]

Optimize[x+y,{(x^2)+(y^2)==1},{{x,-1,1},{y,-1,1}}]

Model Development Step-by-Step

It is easy to modify the objective function or a constraint, to add new constraints, or to change the model functions otherwise in a stepwise fashion. This feature makes easy modular development, as illustrated by the sequence of models shown below. Model changes can also be made simultaneously, of course.

Optimize[(((x-1))^2)+(((y-2))^2),{x-y≥1},{{x,-5,4,5},{y,-5,3,5}}]

Optimize[(((x-1))^2)+(((y-2))^2),{x-y≥1, Sin[2*y]-Cos[x^2]≥1 },{{x,-5,4,5},{y,-5,3,5}}]

Optimize[(((x-1))^2)+(((y-2))^2),{x-y≥1, Sin[2*y]-Cos[x^2]==1 },{{x,-5,4,5},{y,-5,3,5}}]

Optimize[(((x-1))^2)+(((y-2))^2) + Sin[2*y]-Cos[x^2],{x-y≥1, Sin[2*y]-Cos[x^2]==1 },{{x,-5,4,5},{y,-5,3,5}}]

A Model Development Template

The next example shows an easy-to-use model development template. According to the input information requirements of MathOptimizer, we define the following model components:

varswithbounds: the decision variables with given lower and upper bounds, nominal values can be optionally added

objective: the model objective function to minimize

constraints: the model constraints.

varswithbounds = {{x1, -5, 3}, {x2, -5, 4}}; 	
objective = (x1^2 - 2*x2)^2 + (x1 - 2)^2; 
constraints = {Sin[x1^2 - 3*x1*x2] == 0, x1^4 + 5*x2^3 - 40 ≤ 0};

objective 
constraints 
varswithbounds

Optimize[objective, constraints, varswithbounds]

In spite of its small size, this model is not trivial. For comparison, let us solve the same model using the built-in function NMinimize. To this end, we introduce the auxiliary symbols

vars = {x1,x2}
Flatten[{objective, constraints}]

NMinimize[Flatten[{objective, constraints}], vars]

The solution found by NMinimize is not as good as MathOptimizer's solution. This is not always the case, of course. However, the example illustrates an important point: using different approaches (algorithms and corresponding software implementations) to solve difficult numerical problems can be useful.

MathOptimizer Options

List of Options

Options[Optimize]

Options[Optimize]

All default option settings can be overruled (changed) using standard Mathematica conventions.  This point will be illustrated by a set of detailed examples in the following sub-sections.

The combined global-local search approach in MathOptimizer theoretically guarantees stochastic convergence to the global solution (set).  The default option settings are frequently – but not always – sufficient to find the global solution. Generally speaking, in "easy" optimization problems one should always obtain the same solution, except small numerical differences, even when changing MathOptimizer options.  In difficult models – that are typical in the global optimization context – it often makes sense to change one or several of the options, to possibly explore alternative solutions.  Although a general recipe can not be provided, trying a selection of options can lead to improving numerical solutions, with a moderate amount of experimentation.

PenaltyMultiplier

?PenaltyMultiplier

As shown above, the default setting of PenaltyMultiplier is 1: this gives equal weights to the objective and all constraint functions.

In the following example we attempt to calculate a Pythagorean triplet, enforcing integrality using (non-trivial) equality constraints.  A penalty multiplier of 1 is not high enough to give the correct result.

Optimize[Abs[(x^2)+(y^2)-(z^2)],{Sin[π x]==0,Sin[π y]==0,Sin[π z]==0},{{x,1,5},{y,1,5},{z,1,5}},PenaltyMultiplier -> 1]

Raising the penalty multiplier to 2 gives a feasible solution, and 3 gives the expected global solution. However, very high penalties may result in "neglecting" the objective function while aiming "only" at feasiblity. As always, a little experimentation with option settings could help, especially in poorly-scaled or otherwise difficult models – like this unusual formulation to find a Pythagorean triplet.

Optimize[Abs[(x^2)+(y^2)-(z^2)],{Sin[π x]==0,Sin[π y]==0,Sin[π z]==0},{{x,1,5},{y,1,5},{z,1,5}},PenaltyMultiplier -> 2]

Optimize[Abs[(x^2)+(y^2)-(z^2)],{Sin[π x]==0,Sin[π y]==0,Sin[π z]==0},{{x,1,5},{y,1,5},{z,1,5}},PenaltyMultiplier -> 3]

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

Optimize[Abs[(x^2)+(y^2)-(z^2)],{Sin[π x]==0,Sin[π y]==0,Sin[π z]==0},{{x,1,5},{y,1,5},{z,1,5}},PenaltyMultiplier -> 1000]

MultiStarts

?MultiStarts

If used in its global and local search modes, then MathOptimizer first performs MultiStarts (number of) global searches in a stochastic multistart framework, and then it performs local searches starting from the best result or several of the best results found by the global searches. The default setting of MultiStarts is 10; it may be advisable to increase this number for higher-dimensional or otherwise more difficult models. A suggested heuristic setting is m+1+n, where m+1 is number of model functions, and n is the number of decision variables.

It is important to point out that MultiStarts -> 0 leads to local search only, without a preceding global search phase. This local search is started from a given initial point or – if such information is absent, then –  from the midpoint of the search range.

The next example illustrates that in difficult global optimization problems the quality of solution could improve by increasing the number of global scope multistarts.  The example chosen is Trefethen's Problem 4: for details, consult e.g. Weisstein, Eric W. "Hundred-Dollar, Hundred-Digit Challenge Problems." From MathWorld – A Wolfram Web Resource, http://mathworld.wolfram.com/Hundred-DollarHundred-DigitChallengeProblems.html.

We define the objective function

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

The first three calculations below illustrate the multiextremality of this objective function: different starting points chosen from a given box region lead to different solutions.

Optimize[expr,{},{{x,-3,1,3},{y,-3,2,3}},MultiStarts-> 0]

Optimize[expr,{},{{x,-3,-1,3},{y,-3,-2,3}},MultiStarts-> 0]

Optimize[expr,{},{{x,-3,-2,3},{y,-3,-1,3}},MultiStarts-> 0]

The next sequence of calculations shows that

i) without a useful starting point, local search on its own can lead to much poorer solutions than even a quick global sceope search,

ii) the quality of solution typically improves when increasing the number of global scope multistarts,

iii) once the results seem to stabilize, (heuristically) there is no point to increase the search effort further.

Table[Optimize[expr,{},{{x,-3,3},{y,-3,3}},MultiStarts-> i],{i,0,10}]//TableForm

Samples

?Samples

As shown above, the default setting of Samples is 1000.  Again, it may be advisable to increase this number for higher-dimensional or otherwise more difficult models. This point is illustrated by the following example in which we use again Trefethen's Problem 4.

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

Table[{Round[(10^i)],Optimize[expr,{},{{x,-3,3},{y,-3,3}},Samples -> Round[(10^i)]]},{i,1,5,0.5}]//TableForm

RandomSeed

?RandomSeed

The default setting for RandomSeed is 0. Changing the random seed to another (positive) integer may lead to a numerically different solution, and thus can serve as an easy mechanism to generate a sequence of alternative solutions in difficult multimodal problems.

To illustrate this point, consider the following optimization problem that has two global solutions; the corresponding optimized decision variables have opposite signs.  By using different random seeds, both global solutions can be obtained.

Table[Optimize[(x^2)-(y^2),{Cos[x-y]≥ 0.5 },{{x,-5,5},{y,-5,5}},RandomSeed -> i],{i,4}]//TableForm

SamplingMethod

?SamplingMethod

The default setting leads to generating pseudo-random global search points. The alternative low discrepancy generator may produce a more uniform coverage of the box search region, especially when the sample size is relatively limited. Hence, changing the method may lead to a numerically different solution in difficult multimodal problems.

The following example illustrates the basic usage of this option.

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}};

Timing @ Optimize[obj,cons,varswithbounds,SamplingMethod ->"Random"]

Timing @ Optimize[obj,cons,varswithbounds, SamplingMethod -> "QuasiRandom"]

ListableExpressions

?ListableExpressions

Many functions in Mathematica are listable expressions. For example, in the next example Cos is a listable function. The default setting True in general leads to somewhat faster program runs, whenever applicable.

Timing @ Optimize[Cos[x*y*z-(x^2)-Cos[(y^2)+z-x^3]],{Cos[3*x-y+z]≥ 0.5, Cos[x-y+z]≤ 0.5 },{{x,-5,5},{y,-5,5},{z,-15,5}}]

Timing @ Optimize[Cos[x*y*z-(x^2)-Cos[(y^2)+z-x^3]],{Cos[3*x-y+z]≥ 0.5,Cos[x-y+z]≤ 0.5 },{{x,-5,5},{y,-5,5},{z,-15,5}}, ListableExpressions -> False]

If some of the model functions are not listable, then ListableExpressions should be set to False.  For example, the function Max is not listable.

Optimize[Max[x,y],{x+y≤ 3,-x-y≥-3},{{x,1,2},{y,2,6}}, ListableExpressions -> False]

LocalSearches

?LocalSearches

In general, setting larger values in this option could improve the solution quality in difficult optimization problems: for illustration, see the following example.  (You will also notice the corresponding runtime increase.)

Optimize[Max[x,y,z],{Sin[x]+Cos[y-z]≤ 1,Sin[z^2-x^2-y]≤ 0.2},{{x,0,2},{y,0,6},{z,0,10}}, ListableExpressions -> False, LocalSearches->1]

Optimize[Max[x,y,z],{Sin[x]+Cos[y-z]≤ 1,Sin[z^2-x^2-y]≤ 0.2},{{x,0,2},{y,0,6},{z,0,10}}, ListableExpressions -> False, LocalSearches->5]

ShowGlobalResults

?ShowGlobalResults

MathOptimizer performs a number of global searches in a multistart framework, and then uses a local search to improve the result of the best global search.  Although theoretically this approach guarantees global convergence, in numerical practice the best result from the global search phase may not give the best final result after a local search is performed.  ShowGlobalResults enables the user to see the merit function values returned by all the global searches, in order to decide if local searches should be performed starting from other global search results. In other words, setting this parameter to True (default is False) can help to explore a range of initial solution guesses.

Let us remark that since we report here merit function values (as opposed to the best feasible objective function value), it is possible that the reported values are better (for infeasible points) than at the final result.  This is illustrated by the following example.

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}};

Optimize[obj,cons,varswithbounds, ShowGlobalResults -> True]

Consider also the next example: here almost all global search results are already fairly close to each other.  This – together with the coherent final result – indicates the feasibility (or near feasibility) of the global search results.

Optimize[Cos[x*y*z-(x^2)-Cos[(y^2)+z-x^3]],{Cos[3*x-y+z]≥ 0.5,Cos[x-y+z]≤ 0.5 },{{x,-5,5},{y,-5,5},{z,-15,5}}, ShowGlobalResults -> True]

ShowDetailedGlobalResults

?ShowDetailedGlobalResults

Let us point out that this option can be useful, when the model function evaluations are expensive, and we would like to get some useful information about the best solutions found in a resource-limited run.

Another good reason to use it can be when the local search could run into numerical difficulties, due e.g. to local non-differentiability of some model functions. In such cases, one may just do a global search and omit the local search option, by setting LocalSearches -> 0, and RepeatedLocalSearches->False.

The next examples illustrate the ShowDetailedGlobalResults option.

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

Optimize[expr,{},{{x,-3,3},{y,-3,3}}, MultiStarts->1,Samples->100,ShowDetailedGlobalResults->True, LocalSearches->0, RepeatedLocalSearches->False]

Optimize[expr,{},{{x,-3,3},{y,-3,3}}, MultiStarts->10,Samples->10,ShowDetailedGlobalResults->True, LocalSearches->0, RepeatedLocalSearches->False]

Optimize[expr,{},{{x,-3,3},{y,-3,3}}, MultiStarts->15,Samples->20,ShowDetailedGlobalResults->True, LocalSearches->0, RepeatedLocalSearches->False]

The following two runs show that the global search on its own (without using local search) obtains a reasonably good approximate solution, at least in this specific example, in a small fraction of the runtime for the same model solved with local searches.

Optimize[expr,{},{{x,-3,3},{y,-3,3}}, MultiStarts->20,Samples->100,ShowDetailedGlobalResults->True, LocalSearches->0, RepeatedLocalSearches->False]//Timing

The combined global and local search based solution shown below is the true global solution with high precision.

Optimize[expr,{},{{x,-3,3},{y,-3,3}}, MultiStarts->20,Samples->100,ShowDetailedGlobalResults->True, LocalSearches->10, RepeatedLocalSearches->True]//Timing

ShowLocalResults

?ShowLocalResults

It is also possible to inspect all local search results corresponding to the global search results. See the following example and notice that one of the local results returned is inferior to the coinciding others.

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}};

Optimize[obj,cons,varswithbounds, ShowGlobalResults -> True, ShowLocalResults -> True]

Here is another example that shows only the local results. If these are identical, then that serves as a (better) guarantee that the numerical global optimum has been found inded.

Optimize[Cos[x*y*z-(x^2)-Cos[(y^2)+z-x^3]],{Cos[3*x-y+z]≥ 0.5,Cos[x-y+z]≤ 0.5 },{{x,-5,5},{y,-5,5},{z,-15,5}}, ShowLocalResults -> True]

ShowDetailedLocalResults

?ShowDetailedLocalResults

Consider the following example where selecting this option is useful, since it helps to discover several global solutions numerically (when this fact could be overlooked otherwise).

Optimize[x^2-y^2,{Cos[x-y]≥0.5},{{x,-5,5},{y,-5,5}},ShowDetailedLocalResults->True]

LocalSearchMethod

?LocalSearchMethod

In some situations, a better or faster result can be obtained by using the Augmented Lagrangian local search method, rather than using the default FindMinimum saetting: two illustrative examples are shown below.  Note however that this is not always the case, and there may be a need to try both local searches in solving difficult models.

Optimize[Max[x,y,z],{Sin[x]+Cos[y-z]≤ 1,Sin[z^2-x^2-y]≤ 0.2},{{x,0,2},{y,0,6},{z,0,10}}, ListableExpressions -> False, LocalSearchMethod->"FindMinimum",LocalSearches->1]

Optimize[Max[x,y,z],{Sin[x]+Cos[y-z]≤ 1,Sin[z^2-x^2-y]≤ 0.2},{{x,0,2},{y,0,6},{z,0,10}}, ListableExpressions -> False, LocalSearchMethod->"AugmentedLagrangian",LocalSearches->1]

Timing @ Optimize[x+y+z,{(x^2)+(y^2)+(z^2)==1,x(y+1)==(1/2)},{{x,-2,2},{y,-2,2},{z,-2,2}}, LocalSearchMethod -> "FindMinimum"]

Timing @ Optimize[x+y+z,{(x^2)+(y^2)+(z^2)==1,x(y+1)==(1/2)},{{x,-2,2},{y,-2,2},{z,-2,2}}, LocalSearchMethod -> "AugmentedLagrangian"]

Timing @ NMinimize[{x+y+z,(x^2)+(y^2)+(z^2)==1,x(y+1)==(1/2)},{x,y,z} ]

BoundedLocalSearch

?BoundedLocalSearch

The global search phase (by its algorithm design) is always bounded by the stated variable bounds.  However, the local search could leave the feasible region – by violating the set variable bounds – if BoundedLocalSearch is set to False.  Applying this option may lead to better results, if the bounds were (perhaps) mis-specified. The next example illustrates this possibility.

Optimize[(((x-6))^2)+(((y-2))^2),{},{{x,-5,5},{y,-5,5}}]

Optimize[(((x-6))^2)+(((y-2))^2),{},{{x,-5,5},{y,-5,5}},BoundedLocalSearch -> False]

However, if the originally set variable bounds were meaningful but are omitted, then it is possible to obtain infeasible solutions by setting BoundedLocalSearch to False.  Please take another look at the example shown above in this light: if we have to constrain x by 5 indeed, then the second solution is infeasible, while the first one is the valid optimal solution.

RepeatedLocalSearches

?RepeatedLocalSearches

RepeatedLocalSearches is another useful option.  In the next example, setting it to True enables finding the global minimum with less – adaptively chosen, hence a priori unknown – sampling effort, than using the default False setting. As always, it could make sense to test both options in difficult models.

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

Table[{Round[(10^i)], Optimize[expr,{},{{x,-3,3},{y,-3,3}},Samples -> Round[(10^i)],RepeatedLocalSearches -> False]},{i,1,4,0.5}]//TableForm

Table[{Round[(10^i)], Optimize[expr,{},{{x,-3,3},{y,-3,3}},Samples -> Round[(10^i)],RepeatedLocalSearches -> True]},{i,1,4,0.5}]//TableForm

Further Illustrative Applications

A University Math Challenge

The following problem has been posted by the Department of Mathematics at Bilkent University, Ankara, Turkey for July-August 2009.

Find the maximal value of the expression a1/a2+a2/a3+ · · · +a2008/a2009 where a1, a2, . . . , a2009 are positive numbers satisfying ai ≥ a1 + a2 + · · · + ai−1 for each 2 ≤ i ≤ 2009.

We conjecture that the solution is a1=1, and then ai=(2^(i-1)) for 2 ≤ i ≤ 2009, with the corresponding optimum value 1004.5.

Here are a few instances solved numerically by MathOptimizer, to support our guess. In order to obtain the maximum, we need to introduce a change of sign in the objective function, and then compensate for it by multiplying the result by -1.

a1=1 ;
result3= -Optimize[-a1/a2-a2/a3,{a2≥ a1,a3≥a1+a2},{{a2,1,5},{a3,1,5} }]

a1=1 ;
result5= Optimize[-a1/a2-a2/a3-a3/a4-a4/a5,{a2≥ a1,a3≥a1+a2, a4≥ a1+a2+a3,a5≥ a1+a2+a3+a4},{{a2,1,5},{a3,2,4} ,{a4,4,8} ,{a5,5,16}}]

One can make the model development more elegant and compact, by using (automatically generating) indexed variables and a list of constraints, using the model development template discussed earlier.

However, in this particular model this leads to increased runtime and inThe next example shows the 10-variable instance of the same model.

nmax = 10;
a1=1;
varswithbounds = Table[{a[i],i,2^i},{i,1,nmax}]
constraints = Table[a[i]≥Sum[a[j],{j,1,i-1}],{i,2,nmax}]
objective = -Sum[a[i-1]/a[i],{i,2,nmax}] 
Timing[Optimize[objective,constraints,varswithbounds]]

It is emphasized that this model is not really meant to be solved numerically, and we use it only as an illustrative example. (The proper theoretical approach is to use induction.)

Optimization Models with Arbitrary (Continuous) Mathematica Functions

MathOptimizer can handle a broad range of user-defined model functions, including Mathematica's own functions and their suitable customized extensions.  In principle, one could use all continuous Mathematica functions as model components. A technical consideration is to define the model functions so that they are only evaluated numerically, and for general safety (if in doubt, then) to set ListableExpressions to False.

For illustration, in the next example we want to find the unknown parameter 0 ≤ a ≤ 3 in the parametric integral (∫_0^Π)cos(a (x^2))sin (a x)dx so that the corresponding integral is minimal. We then compare our result to the result obtained by NMinimize. Since solving this model takes more than a few seconds, we also time the runs.

f[a_?NumberQ]:= NIntegrate[Cos[a x^2]*Sin[a x],{x,0,π}]

Timing @ NMinimize[{f[a],0 ≤ a ≤ 3},{a}]

Timing @ Optimize[f[a],{},{{a,0,3}},ListableExpressions -> False, MultiStarts -> 1]

Notice that in this example MathOptimizer finds the correct global solution (unlike NMinimize) as shown by the figure below. This is not always the case, of course. However, the example illustrates an important point: using different approaches (algorithms and corresponding software implementations) to solve difficult numerical problems can be useful – especially so in the area of nonlinear optimization.

The runtime is about 2.7 times higher for MathOptimizer: this is not always the case as we will illustrate on circle packing test examples in the next section.

Plot[f[a],{a,0,3}]

Non-uniform Circle Packing Problems

Introduction and Model Formulation

Optimized object packings (configurations or arrangements) are important in various engineering and scientific fields such as numerical integration, potential energy models, experimental design and others.  Next, we will illustrate MathOptimizer performance by solving non-trivial circle packing problems.  Specifically, we will study the following general problem-type: given an arbitrary collection of circles, find the smallest circumscribing circle that includes all circles in a non-overlapping arrangement. This and similar problems have been tackled recently by a number of researchers. In general, the global solution is unknown, hence only putative optima are published (except in special cases and/or for small model instances).

The function pack shown below generates the objective function, constraints and variables with bounds for this problem. The first argument of the function is the number of circles; the second argument is a function determining the size-sequence of the packed circles.  The third argument is the target (minimized) radius for the circumscribing circle. See also the comments embedded in the code. Notice also that we have to clear the symbol obj that was used earlier.

Clear[obj];

pack[n_Integer,func_,rmax_]:=
Module[
{r,cons,cons1,cons2,vars,lb,ub,soln,rcirc,pd, bds,noms, vwb},
Do[r[i]=func[i],{i,n}];
(* The circles don't overlap *)
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}]];
(* The objective function is ≥ the furthest point on each circle from the origin *)
cons2 = Table[obj -r[i] ≥(Sqrt[ ((x[i])^2)+((y[i])^2)]),{i,n}];
(* The variables *)
vars=Join[Flatten[Table[{x[i],y[i]},{i,n}]],{obj}];
(* All the constraints *)
cons = Join[cons1,cons2];
(* Generate the bounds on the variables *)
bds = (rmax-#)& /@Table[r[i],{i,n}];
ub = Join[bds,bds,{rmax}]; (* Upper Bounds *)
lb= Join[-bds,-bds,{0}]; (* Lower Bounds *)
noms = RandomReal[{-0.5,0.5},Length[lb]];
vwb=Transpose[{vars,lb,ub}];
{obj,cons,vwb}
]

The following function serves to plot the result obtained.

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}]
]

3-circle Example

In all numerical examples below, the n-circle problem instances are defined by circles with radii given as ((Range[n])^((-1)/2)). That is, the radii are 1, 1/(Sqrt[2 ]), 1/(Sqrt[3]), 1/2, 1/(Sqrt[5]), and so on.. This choice guarantees that the total area to cover (by the circumscribing circle) increases unbounded as the number of circles grows.

The first example involves packing 3 circles.  There are a number of solutions with rotational symmetry, but all have the same objective value, approximately 1.70997.  First we define the 3-circle model instance, then we solve it using MathOptimizer.

prob3 = pack[3,((#)^(-1/2))&,2.];

Timing  @ (res=Optimize[Sequence @@ prob3])

plt[res,3,((#)^(-1/2))&,2.]

Since all the local minima have the same objective function value, the result can be found more quickly using only one global search and one local search phase.

Timing  @ Optimize[Sequence @@ prob3, MultiStarts -> 1, LocalSearches -> 1]

NMinimize finds a solution with the same objective function.

Timing @ NMinimize[{prob3[[1]],Sequence @@ prob3[[2]]},prob3[[3,All,1]]]

5-circle Example

When packing 4 or more circles, there are a number of local minima.  This provides a good opportunity to check out some of the various MathOptimizer options.

prob5= pack[5,((#)^(-1/2))&,2.];

Using a value of 1 for MultiStarts gives an inferior solution, as shown in the next two examples.

Timing @(res= Optimize[Sequence @@ prob5, MultiStarts -> 1, LocalSearches -> 1])

plt[res,5,((#)^(-1/2))&,2.]

The default usage of MathOptimizer gives a better result:

Timing @ (res=Optimize[Sequence @@ prob5])

plt[res,5,((#)^(-1/2))&,2.]

After some experimentation, our best result has been obtained by setting BoundedLocalSearch to False (enabling the circles to move out and then back in again in the course of the optimization) and using RepeatedLocalSearches. See the results below.

Timing @ (res=Optimize[Sequence @@ prob5, BoundedLocalSearch -> False, RepeatedLocalSearches -> True])

plt[res,5,((#)^(-1/2))&,2.]

For comparison, we also ran NMinimize: it finds the same objective value.

Timing @ NMinimize[{prob5[[1]],Sequence @@ prob5[[2]]},prob5[[3,All,1]]]

10-circle Example

Packing more and more circles quickly becomes a real computational challenge, since the number of constraints grows quadratically, and the number of local minima grows exponentially. Next, we look at the 10-circle model instance

prob10 = pack[10,((#)^(-1/2))&,2.];

Timing @ (res=Optimize[Sequence @@ prob10])

plt[res,10,((#)^(-1/2))&,2.]

For comparison, we ran also NMinimize: for this problem (using default settings), it finds a somewhat inferior solution to the result obtained by MathOptimizer.

Timing @ NMinimize[{prob10[[1]],Sequence @@ prob10[[2]]},prob10[[3,All,1]]]

15-circle Example

prob15 = pack[15,((#)^(-1/2))&,2.];

Timing @ (res=Optimize[Sequence @@ prob15,BoundedLocalSearch -> False, RepeatedLocalSearches -> True])

plt[res,15,((#)^(-1/2))&,2.5]

For comparison, we ran NMinimize for this problem as well, using default settings. Again, it finds a somewhat inferior solution to the result obtained by MathOptimizer.

Timing @ NMinimize[{prob15[[1]],Sequence @@ prob15[[2]]},prob15[[3,All,1]]]

20-circle Example

As a final illustrative example shown here, we look at the 20-circle model instance

prob20 = pack[20,((#)^(-1/2))&,2.5];

Timing @ (res=Optimize[Sequence @@ prob20,BoundedLocalSearch -> False, RepeatedLocalSearches -> True])

plt[res,20,((#)^(-1/2))&,2.5]

Again, MathOptimizer has found a somewhat better solution than NMinimize, see below.

Timing @ NMinimize[{prob20[[1]],Sequence @@ prob20[[2]]},prob20[[3,All,1]]]

References

This list includes a few classical textbooks, more books from recent years, as well as a few of our topical articles. Several of the works cited focus on mathematical concepts, theory and algorithms, or software issues pertinent to the subject, while others discuss also a number of interesting applications. In addition to the Mathematica Book, several well-received other books on Mathematica are also listed.

Bazaraa, M.S., Sherali, H.D. and Shetty, C.M. (1993) Nonlinear Programming: Theory and Algorithms. Wiley, New York.

Bhatti, M.A. (2000) Practical Optimization Methods With Mathematica Applications. Springer-Verlag, New York / Berlin / Heidelberg.

Bertsekas, D.P. (1999) Nonlinear Programming. (2nd edn.) Athena Scientific, Cambridge, MA.

Bomze, I.M., Csendes, T., Horst, R., and Pardalos, P.M., eds. (1997) Developments in Global Optimization. Kluwer Academic Publishers, Dordrecht / Boston / London.

De Leone, R., Murli, A., Pardalos, P.M. and Toraldo, G., eds. (1998) High Performance Software for Nonlinear Optimization. Kluwer Academic Publishers, Dordrecht / Boston / London.

Dixon, L.C.W. and Szegö, G.P., eds. (1975, 1978) Towards Global Optimisation. Vols. 1-2. North-Holland, Amsterdam.

Edgar, T.F., Himmelblau, D.M, and Lasdon, L.S. (2001) Optimization of Chemical Processes. McGraw-Hill, New York.

Floudas et al. (1999) Handbook of Test Problems in Local and Global Optimization. Kluwer Academic Publishers, Dordrecht / Boston / London.

Fourer, R. (2002) Nonlinear Programming FAQ. http://www-unix.mcs.anl.gov/otc/ Guide/faq/nonlinear-programming-faq.html.

Gaylord, R.J., Kamin, S.N. and Wellin, P.R. (1996) An Introduction to Programming with Mathematica. (2nd edn.) Springer-Verlag, New York / Berlin / Heidelberg.

Grossmann, I.E., ed. (1996) Global Optimization in Engineering Design. Kluwer Academic Publishers, Dordrecht / Boston / London.

Hendrix, E.M.T. (1998) Global Optimization at Work. Ph.D. Dissertation, Agricultural University, Wageningen.

Hillier, F.S. and Lieberman, G.J. (1986) Introduction to Operations Research. (4th edn.) Holden Day, CA.

Hock, W. and Schittkowski, K. (1981) Test Examples for Nonlinear Programming Codes. Springer-Verlag, Berlin / Heidelberg / New York.

Horst, R. and Pardalos, P.M., eds. (1995) Handbook of Global Optimization, Vol. 1. Kluwer Academic Publishers, Dordrecht / Boston / London.

Journal of Global Optimization (published since 1991). Kluwer Academic Publishers, Dordrecht / Boston / London.

Journal of Heuristics (published since 1995). Kluwer Academic Publishers, Dordrecht / Boston / London.

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.

Kearfott, R.B. (1996) Rigorous Global Search: Continuous Problems. Kluwer Academic Publishers, Dordrecht / Boston / London.

Maeder, R.E. (1997) Programming in Mathematica. (3rd edn.) Addison-Wesley, Reading, MA.

Mittelmann, H.D. and Spelucci, P. (2002) Decision Tree for Optimization Software. http://plato.la.asu.edu/guide.html.

Papalambros, P.Y. and Wilde, D.J. (2000) Principles of Optimal Design, Cambridge University Press, UK.

Pardalos, P.M. and Romeijn, H.E., eds. (2002) Handbook of Global Optimization, Vol. 2. Kluwer Academic Publishers, Dordrecht / Boston / London.

Pintér, J.D. (1996) Global Optimization in Action (Continuous and Lipschitz Optimization: Algorithms, Implementations and Applications). Kluwer Academic Publishers, Dordrecht.

Pintér, J.D. (2001) Computational Global Optimization in Nonlinear Systems – An Interactive Tutorial. Lionheart Publishing, Inc., Atlanta, GA.

Pintér, J.D. (2002) Global Optimization: Software, Tests and Applications. Chapter 15 (pp. 515-569) in: Pardalos and Romeijn, eds.

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

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.

Schrage, L. (2001) Optimization Modeling with LINGO. LINDO Systems, Inc. Chicago, IL.

Wagner, D.B. (1996) Power Programming with Mathematica – The Kernel, McGraw-Hill, New York.

Wolfram, S. (2005) The Mathematica Book. (5th edn.) Wolfram Media, Champaign, IL, and Cambridge University Press, Cambridge, UK.