GLOBAL OPTIMIZATION 9.0

GLOBAL NONLINEAR OPTIMIZATION USING MATHEMATICA

Loehle Enterprises
1258 Windemere Ave.
Naperville, IL 60564
USA
630-527-8402
630-416-9902 fax
info@loehleenterprises.com
www.loehleenterprises.com

Version 9.0

(C) Copyright 1998-2012 Loehle Enterprises.  No part of this manual or software may be copied or reverse engineered except that purchasers of this package may make back up copies of the software.

*****************************************************************
DISCLAIMER

*****************************************************************

NOTICE: The software product described in this manual is a research tool.  No warranty is offered as to operability of this product.  No guarantee can be made that any particular nonlinear equation can be solved.  No safety critical or financial decisions should be based on the operation of this software.  The manufacturer accepts no liability for losses resulting from its use.
*****************************************************************

i  CONTENTS

TABLE OF CONTENTS

  I	Getting Started                                                                                             	  3

	1  Welcome                                                                                                    	  3
	2  Registering                                                                                                 	  3
	3  Installation                                                                                                    4
	4  A Quick Start and Example							 4

II	Tips for Performance                                                                         		10

	1  Improving Performance                                                                              10
	2  Parallel Computing								11
	3  The Compile Function                                                                               	13
	4  User Functions: Working with Named Functions				14
	
 III	General Nonlinear Optimization: GlobalSearch, IntervalMin, 
 	and GlobalPenaltyFn								17

	1  Introduction								17
	2  Example									18
	3  Program Operation.								19
		A  Parameters and Default Values					19
		B  Bounds								19
		C  Constraints								19
	4  Output									24
	5  Maximization								25
	6  Limitations									25
	7  Error Messages								27
	8  Performance/Speed								29
	9  Testing and Examples							33

	
V	Nonlinear Regression: the NLRegression Function              			39

	1  Introduction								39
	2  Utilizing chi-square fit criteria						46
	3  Multiple Independent Variables Regression				46

V	Maximum Likelihood Estimation: the MaxLikelihood Function		52

	1  Introduction								52
	2  Examples and Built-In Functions						52

VI	Discrete Variable Problems: 
		InterchangeMethod & TabuSearchMin				58

	1  Introduction								58
	2  Applications								59
		A  Capital Allocation							59
		B  Vehicle Routing/Travelling Salesman				60
		C  Minimum Spanning Tree						64
	
 VII	The MaxAllocation Function							65

	1  Introduction								65
	2  Applications								68
		A Investment Allocation Problems					68
		B Derivative Hedging with Quadratic Programming			69

VIII 	Applications									70
								
	1  Zeros of a Function/Roots of a Polynomial				70
	2  Integer Programming: the Knapsack Problem				72
	3  Differential Equation Models						73
	
IX	Literature Cited								75

I  GETTING STARTED

I.1  Welcome

Welcome to Global Optimization.  This package provides a suite of tools for solving nonlinear optimization problems, as well as a variety of other applications such as finding the roots or zeros of a nonanalytic function.  The package is easier, simplier, and more robust than most optimization tools, and you will find yourself working more efficiently and more easily than before.  This package solves many problems competing products either do not solve at all, or yields a better solution.  It handles any combination of Real, Integer, and Binary (0-1) variables.  It is particularly effective for problems in finance, engineering, model estimation, biochemistry and other difficult fields where other tools are unable to perform well.  It accepts black-box models.  The user does not need to give a feasible starting point.  The package finds multiple solutions if they exist, which competing tools do not do.

I.2  Registering

Before you continue please register by emailing your name and address with a notice that you have purchased Global Optimization to info@loehleenterprises.com.  This will enable us to reach you with free updates and other news.  The developer does not receive customer information when the package is purchased from Wolfram Research or from distributors.  Please include the version of Global Optimization that you purchased and the machine you are using it on.

I.3  Installation

I.3.A  Recommended System Capabilities

This package is designed to work with Mathematica, which must be installed on your computer.  The following are recommended system capabilities for installation:

	•	Personal computer, 600 MHz or faster
	•	200Mb free RAM or more
	•	Mathematica 7.0 or higher installed

I.3.B  Installation on a Hard Drive

To install the package on your hard drive on any machine, copy the contents of the disk to the Mathematica directory.  The file may need to be unzipped.  On Mac the file may be a self-extracting archive; double-click it to extract the contents.  The go90v.mx file (where v denotes the version) is an encoded Mathematica package.  This encoded file can not be read or opened and will not work on other platforms.  It needs to be placed in the main Mathematica directory or in the AddOns/Applications directory.  An example notebook file .nb is included for illustration.  The .nb file(s) may be placed anywhere.

I.4  A Quick Start and Example

Version 9.0 contains nine functions: GlobalSearch, IntervalMin, GlobalPenaltyFn, GlobalMinima, NLRegression, MaxLikelihood, InterchangeMethodMin, TabuSearchMin, and MaxAllocation.  The package has been designed for easy use with Mathematica.  It uses all valid Mathematica functions and syntax.  The user must have Mathematica installed.  This manual assumes basic familiarity with Mathematica.  Input files should be in the form of Mathematica notebooks, with a file designation "filename.nb". This manual is executable.  Timing examples in this manual are based on a 3.4 Ghz Pentium IV machine using Mathematica 9.0.

	To begin the execution of a notebook, follow this procedure:

1.	Start Mathematica.

2. 	Use the File pulldown menu to Open one of the example notebook files supplied with the disk.  

3.	Once this notebook is open, pull down the Kernel menu and select Evaluate Notebook.

	In the examples, the package is installed with the Get (Get["go90v.mx"] or <<go90v.mx) command, where "v" is the current version (a, b, etc) sent to the user. This notebook should execute with no errors and should produce output that resembles closely what was in the notebook when you opened it.  If not, some error may exist.  The basic formats are similar to the Mathematica function NMinimize.  The following definitions define the functions in the package:

<<go90v.mx

?GlobalSearch

?GlobalPenaltyFn

?IntervalMin

?MaxAllocation

?MaxLikelihood

?NLRegression

?InterchangeMethodMin

?TabuSearchMin

The following options are available for one or more functions.  PenaltyMethod is used for GlobalPenaltyFn.

?ParallelOption

?CompileOption

?EvaluateObj

?PenaltyMethod

?TabuListLength

?MaxIterations

?SensitivityPlots

?FunctionCalls

?ShowProgress

?SimplifyOption

?Starts

?StartsList

?UserResiduals

?Weights

GlobalSearch[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

GlobalPenaltyFn[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

IntervalMin[expression,inequalities,Null,{InitialPoint__Interval},varlist,tolerance,options]

MaxAllocation[f[var1name,...],rhs,{varlist}..,tolerance,options]

NLRegression[data,expression,independent_variables,inequalities,{{var1name,lowbound,highbound}..},tolerance,options]

MaxLikelihood[data,expression,independent_variables,inequalities,{{var1name,lowbound,highbound}..},tolerance,options]

InterChangeMethodMin[expression,inequalities,varlist,tolerance,options]

TabuSearchMin[expression,inequalities,varlist,tolerance,options]

Where "expression" is the equation or function name for the function to be minimized, "inequalities" is a list of inequality constraints, and "equalities" is a list of equalities.  Speed of execution is enhanced with Mathematica versions 7 and higher.  The Mathematica kernel and the executable program need not reside on the same machine.

The functions in this package are designed for nonlinear problems with local minima, multiple minima, and/or linear or nonlinear constraints.

The function GlobalSearch is a multiple-start generalized hill-climbing algorithm designed to work with or without constraints.  It is robust to noisy functions and local minima.  It can handle large problems (200+ variables).  It can handle linear and nonlinear equality and inequality constraints, which are assumed to be analytic.   GlobalSearch is the general purpose optimizer in the package, and is the engine used for the NLRegression and MaxLikelihood functions.  Variables can be Real or Binary (0-1).  If more than 10 binary variables are included in a problem, GlobalPenaltyFn with PenaltyMethod->True or the discrete functions InterchangeMethodMin or TabuSearchMin should be used.

	The function IntervalMin is a general minimizer using Interval methods.  It can handle inequality constraints, which must be analytic.  It is robust to local minima, but is slower than other methods.

The function GlobalPenaltyFn is a multiple-start generalized hill-climbing algorithm designed to work with or without constraints.  It handles the special case of constraints that are not Solvable for any of the variables, or that are black box.  It is robust to noisy functions and local minima.  It can handle large problems (200+ variables).  For equality constraints, it has two options.  The default attempts to work with equality constraints analytically.  PenaltyMethod when True uses a penalty method.  This option is useful when equality constraints are extremely complicated or non-analytic.  Variables can be Real, Integer or Binary (0-1).  If more than 10 binary variables are included in a problem, GlobalPenaltyFn with PenaltyMethod->True or the discrete functions InterchangeMethodMin or TabuSearchMin should be used.

The function MaxAllocation is designed for allocation problems, such as arise in finance.  In an allocation problem, a fixed quantity, such as an investment sum, is to be distributed to a portfolio.  All variables must be nonnegative and the sum of all investments must equal the total available for investing.  This creates a nonlinear problem with a single equality constraint.

The function NLRegression solves nonlinear regression problems.  It is particularly designed for noisy problems or those requiring constraints to achieve a good fit.  L1 and L2 norms are allowed.  Least-squares and chi-square options are also available.  Sensitivity plots and confidence intervals on the parameters are provided as output.

The function MaxLikelihood solves maximum likelihood estimation problems using Log-likelihood estimation.  A library of common distributions is provided and fit statistics are computed.  Summary statistics are computed.

	The function InterchangeMethodMin uses the interchange method to search for a solution to binary problems.
	
	The function TabuSearchMin uses the tabu method to search for a solution to binary problems.

A basic example is:

GlobalSearch[x^2+y^2,,,{{x,-5,5},{y,-5,5}},.000001,Starts->1]

With a constraint:

GlobalSearch[x^2+y^2,{-x+1,-y+1},,{{x,-5,5},{y,-5,5}},.000001,Starts->1]

With an equality constraint:

GlobalSearch[(x-1.3)^2+(y-1.5)^2,{-x+1,-y+1},{x-y},{{x,-5,5},{y,-5,5}},.000001,Starts->1]

With an equality constraint and x binary:

GlobalSearch[(x-1.3)^2+(y-1.5)^2,{-x+1,-y+1},{x-y},{{x,-5,5,Binary},{y,-5,5}},.000001,Starts->1]

II  TIPS FOR PERFORMANCE

A key factor for numerical optimization is the time taken to compute the objective function.  This is critical because the objective function may be called thousands of times to hundreds of thousands of times.  This section discusses tricks to achieve better performance.

II.1  Improving Performance

Some general tips include the removal of the reading of external files from the user function, computing constant expressions once and then using the result, and avoiding logical operations that will prevent Compile from being used.  A key is to generate a function once as an explicit expression, and then pass this in to the package.  If many steps are involved in generating a user function and these must be performed every time the function is called, this can be very expensive.  It often happens that parts of an expression involve constants.  If these can be evaluated up front, the savings can be enormous.  In the following example, the function "f" evaluates numerically when Evaluate is executed during Function creation.  In function "g" functions like Sin and Tan are evaluated each time the function is called:

ClearAll[f,g]

x={x1,x2,x3,x4,x5,x6,x7,x8,x9,x10};

f=Function[Evaluate[x],Evaluate[Sum[Sin[Pi/2.1]*x[[i]]+Tan[19.4*Pi]/x[[i]]+Cos[.1]+Sqrt[99],{i,1,10}]]]

g[p_]:=Module[{d},{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10}=p;r=Sum[Sin[Pi/2.1]*x[[i]]+Tan[19.4*Pi]/x[[i]]+Cos[.1]+Sqrt[99],{i,1,10}];Return[r]]

dat=Table[i,{i,1,10}]

Apply[f,dat]

g[{1,2,3,4,5,6,7,8,9,10}]

Do[Apply[f,dat],{i,1,500}]//Timing

Do[g[{1,2,3,4,5,6,7,8,9,10}],{i,1,500}]//Timing

Because the constants in f were Evaluated before execution (during Function definition), f was very fast to execute, even without Compile.  In contrast, g, which computes Sin etc. every time it is called, took 46 times longer to compute.

II.2  Parallel Computing

Mathematica 7.0 introduces the use of parallel computing.  Version 9.0 of Global Optimization takes advantage of this capability.  In the functions GlobalSearch and GlobalPenaltyFn, the multiple starts can be solved with multiple processors if available. These multiple starts help to ensure an optimal solution or to find multiple solutions if they exist. For problems that are very costly to evaluate, this can provide a substantial speed-up in execution.  If Starts is set to 1, the job will not run in parallel.  To run in parallel mode, use ParallelOption -> True.  ShowProgress can still be set True, but multiple outputs will result (1 for each processor).

You can check how many processors your computer has with :

n=$ProcessorCount

On machines with single processors, Mathematica can still launch up to 4 Kernels, but performance will not be improved.  The Global Optimization package can launch kernels from within a run, but recognizes if kernels are already open. In the example below

LaunchKernels[n]

It is necessary to load the Global Optimization package for each kernel, using :

ParallelEvaluate[<<go90v.m]

where the "v" in the file name is replaced with the file version you have (such as "a" etc).  The example below has multiple solutions, nicely solved using the parallel capabilities:

b=GlobalPenaltyFn[(y-x^2*5.1/(4.*π^2)+5.*x/π-6.)^2+10.*(1.-1./(8.*π))*Cos[x]+10.,{-40-x,x-40,-40-y,y-40},{},{{x,-20,30},{y,0,40}},.000001,Starts->100,ParallelOption->True,ShowProgress->False];

Length[b]

c={};Do[AppendTo[c,{x,y}/.b[[ii,1]]],{ii,1,Length[b]}]

ListPlot[c,PlotRange->{{-20,30},{0,40}}]

For this particular case, the amount of message passing leads to only a small improvement in performance with two processors.  In general, the longer it takes to compute the function and the more variables, the more advantage one obtains from using the parallel computing option.

If one is using a named function such as "func" for the objective (see following sections), particularly if it is a Module, one also needs to distribute the definitions:

DistributeDefinitions[func]

When finished with parallel computing, the extra kernels can be closed.

CloseKernels[]

II.3  The Compile Function

The Mathematica function Compile is used to speed up program execution.  The user input function and constraints are both Compiled.  The advantage of the Compile function is greatest with more complex expressions and larger problems.  For small problems, it can be difficult to even detect a benefit, but there is little cost to the Compile.  The range of speed improvement due to Compile can range from 30% to 30x.  Compilation does not work for all Mathematica expressions.  For example, if "Apply" or "If" are used in the function definition, the function will not Compile.  It is suggested that if a complex expression involving Mathematica special functions is to be used, the user attempt to Compile the expression to test it.  The CompileOption can be set to False (CompileOption->False) in the function call if the user function is not compilable.  It may also be set to False if the user Compiles the expression before passing it in.  It is possible to get errors during the Compile step.  These will be highlighed in blue and will indicate that they are compilation errors.  It is also possible that execution errors can result if non-machine numbers result from some computation of the compiled function.  Mathematica usually reverts to the uncompiled expression in such cases, but the user may wish to run with CompileOption->False in this case to obtain the most accuracy and best speed.  

	We can see the benefit of the Compile in the following example:

ClearAll[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10]

x={x1,x2,x3,x4,x5,x6,x7,x8,x9,x10};

f=Sum[(1+x[[i]])^2,{i,1,10}]

g1=Function[Evaluate[x],Evaluate[f]];

g2=Compile[Evaluate[x],Evaluate[f]];

Do[Apply[g1,{.1,2.,3.,4.,3.,.3,.3,.3,.5,.5}],{i,1,10000}]//Timing

Do[Apply[g2,{.1,2.,3.,4.,3.,.3,.3,.3,.5,.5}],{i,1,10000}]//Timing

We see that Compile speeded up the execution by a factor of more than 4.  On some functions, the speedup can be a factor of 20 or more.  Setting up a function as a Module prevents Compile from working on the objective function, and should be avoided.  Some Mathematica functions, particularly those that are math functions, can be included in a Compiled function.  Examples include Abs, Sqrt, Max, Min, Sin, etc.  Others, such as FindRoot, can not.

II.4  User Functions: Working with Named Functions

The objective (goal) function in all the functions in this package can be defined in several ways.  The function must in some cases be passed in differently, depending on the function format.  In the simplest case, the expression is itself passed in.  Here we find the minimum of the given function, which is zero.

ClearAll[x]

GlobalSearch[x^2+y^2,,,{{x,38,40},{y,38,40}},.00000001,Starts->1]

We can also pass in a function by name:

f:=x^2+y^2

GlobalSearch[f,,,{{x,38,40},{y,38,40}},.00000001,Starts->1]

In order to Compile or use the objective function, the programs in this package will attempt to Evaluate the expression.  If there are certain operations involved in the user function, such as obtaining the Inverse or Determinant of a matrix, Mathematica will attempt to Evaluate these symbolically when the function is used in the program.  It will then work with the Evaluated function.  If the user expression is complex, however, this may become impossible.  For example, a symbolic inverse of a matrix with 10 variables may take a very long time.  The user can test this by trying to Evaluate the objective function before passing it in to any of the functions of this package.  If it fails to Evaluate, then use the parameter EvaluateObj->False to prevent evaluation.  Another time when EvaluateObj should be set to False is when any logical operations are performed in the objective function, as in the following example:

f=Function[{x,y},If[x+y<0,tem=0,tem=y+x];tem];

Apply[f,{-1,-1}]

Apply[f,{1,1}]

Clearly, the function is minimized for any set with x + y negative or for x==-y.  Because of the conditional logic in the expression, it should not be symbolically evaluated.  This means that it should not be simplified (SimplifyOption->False) or Evaluated (EvaluateObj->False):

GlobalSearch[f,{},,{{x,38,40},{y,38,40}},.000001,Starts->1,EvaluateObj->False,SimplifyOption->False]//Timing

GlobalPenaltyFn[f,{},,{{x,38,40},{y,38,40}},.000001,Starts->1,EvaluateObj->False,SimplifyOption->False]//Timing

If, in contrast, we set up a function as a Module, the function can not be Compiled or Evaluated without giving an incorrect answer.  However, this problem is checked for internally and does not cause a problem, though it does cause a problem in Mathematica functions such as NMinimize.  In this case, it is also critical to use ?NumericQ in the Module definition to prevent symbolic evaluation of the expression, which is not valid for many functions such as those using logic or in which numerical quantities are computed.  Global Optimization will return an error if this has not been done.

ClearAll[x,y,g]

u4[x_]:=Module[{y},If[x>-2,y=3+x^2,y=0];Return[y]]

GlobalSearch[u4[x],,,{{x,38,40}},.000001,Starts->1]

NMinimize does not trap this error.

NMinimize[u4[x],{{x,38,40}}]

ClearAll[x,y,u4]

u4[x_?NumericQ]:=Module[{y},If[x>-2,y=3+x^2,y=0];Return[y]]

GlobalSearch[u4[x],,,{{x,38,40}},.000001,Starts->1]

The valid answer for u4 is any value less than -2, but NMinimize fails to work properly.

NMinimize[u4[x],{{x,38,40}}]

Note that some functions if evaluated without the ?NumericQ, will be evaluated symbolically, and will use a huge amount of memory and take a very long time, possibly even crashing the kernel.

III  GENERAL NONLINEAR OPTIMIZATION: GlobalSearch, GlobalPenaltyFn, and IntervalMin

III.1  Introduction

Three general nonlinear solvers are included in this package.  They each are designed for a particular class of problems.  These functions provide tools for global optimization.  Traditional gradient (local) approaches require the user to know how many optima are being sought for the function to be solved, and roughly where the optima are, so that a good initial guess can be made.  The user, however, rarely can make a good initial guess and usually has no information about the existence of multiple solutions.  In the absence of such information, existing algorithms will generally converge to an answer, but this may be only a local solution that may not be globally optimal.  Furthermore, even global optima may not be unique.  The functions in this package can solve black-box objective functions, objective functions that are nondifferentiable, differential equation models, and functions with discrete steps in them.  The initial bounds on the parameters given on input do not need to bound the true solution, but only help get the algorithms started in a good region.   The solvers work for problems where NMinimize and competing products fail.

GlobalSearch approaches the difficult problem of finding a global optimum with several techniques.  A generalized hill climbing technique is used that is based on Newton's method but using a generalized gradient rather than a derivative, and allowing for constraints.  Constraints must be analytic, but can be linear or nonlinear.  Multiple starts are used to test for the existence of multiple solutions.  The multiple starts are generated randomly from the region defined by the range of parameter values input by the user.  Feasible starting regions are not needed.  When a step can not be made that improves the solution by at least "tolerance", the program terminates.   Discrete 0-1 variables can be mixed with Real variables (assumed if not specified).  GlobalSearch is the solver used by the regression and maximum likelihood functions.  The function is defined by

GlobalSearch[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

GlobalPenaltyFn addresses the problem of constraints that are nonanalytic. If GlobalSearch fails on the constraints, use GlobalPenaltyFn.  A generalized hill climbing technique is used that is based on Newton's method but using a generalized gradient rather than a derivative, and allowing for constraints.  An adaptive penalty method is used that adjusts the penalties for constraint violation depending on the degree of violation of constraints.  Multiple starts are used to test for the existence of multiple solutions.  The multiple starts are generated randomly from the region defined by the range of parameter values input by the user.  Feasible starting regions are not needed.  When a step can not be made that improves the solution by at least "tolerance", the program terminates. Two methods are available for solving equality constraints.  An analytic method is default.  To override this method, use PenaltyMethod->True.  The alternate method is a penalty method for both inequalities and equalities. An optional 4th parameter in the variable bounds list can be used to define a variable to be Integer or Binary (as in {{var1name,lowbound,highbound, Integer}..}  ).  The function is defined by

GlobalPenaltyFn[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

IntervalMin uses Interval methods to find a solution.  Intervals are computed in Mathematica based on concepts of limits.  An initial Interval is input by the user.  It is beneficial for this Interval to cover the true solution point, but this is not necessary as long as it is not orders of magnitude distant.  Both the input objective function and the constraints must be analytic.  IntervalMin is robust to local minima, but is not as fast as GlobalSearch because Compile can not be used with Interval variables.  IntervalMin can not solve problems with equality constraints at this time.  We illustrate Intervals below:

x^2/.x->Interval[{-1,2}]

In this case, the true minimum of the function is given by the lower bound of the Interval result.  In simple or separable cases, this will pertain.  In other cases it will not:

1.`-2.` x+(x^2)+100.` (x^4)-200.` (x^2) y+100.` (y^2)/.{x->Interval[{0,2}],y->Interval[{0,2}]}

The true minimum of this function is 0 at {1,1}.  The true minimum is within the Interval of the result in this case.  Note that when IntervalMin is used, both the estimated parameters and the function value are given in terms of an Interval.  The format for the function is given by:

IntervalMin[expression,inequalities,Null,{InitialPoint__Interval},varlist,tolerance,options]

III.2  Example

A simple example is to minimize a square function subject to a constraint, as follows.  The constraint

x>1&&y>1

must be converted to standard form as a list:

-x<-1&&-y<-1

-x+1<0&&-y+1<0

c={-x+1,-y+1};

GlobalSearch[x^2+y^2,c,,{{x,-5,5},{y,-5,5}},.000001]//Timing

The output shows four solutions from random initial points, the default, with the {x,y} value and the function value in each list element.  Initial feasible points are not needed.  In this case, the function is smooth, so all three solutions are the same, and there is really no need to use multiple starts.

III.3  Program Operation

III.3.A  Parameters and Default Values

The tolerance (T) defines the amount of improvement (absolute, not relative) in the function value required at each iteration.  If at least this much improvement is not found, the program stops.  The tolerance can be set to very small values, such as (10^(-15)), to achieve a highly accurate estimate of the minimum.  If the user sets tolerance = 0, this is an error.

The functions have seven options.  The defaults are MaxIterations->10000, CompileOption->True, ShowProgress->True, StartsList->{}, EvaluateObj->True, SimplifyOption->True,  ParallelOption->False, FunctionCalls->False, and Starts->4.  MaxIterations and other checks prevent the program from running away to infinity when a mistake is made in the function (Min[x] for example).  For very computationally expensive problems (like those with many parameters), it is useful to use Starts->1 to evaluate the time a run takes.  For general usage, Starts->4 (the default) is good. The results of this short run can then be used to define further, more directed runs based on whether all 3 starts found the same solution.  A list of starting values can be input instead of letting the program find them with random search (e.g., StartsList->{{1.,2.},{2.,1.}}).  The CompileOption determines whether the user function and constraints will be compiled.  While Compile reduces execution time, some functions can not be compiled, in which case CompileOption->False should be used.  If the objective function should not be Evaluated, use EvaluateObj->False.  SimplifyOption attempts to simplify the objective function.  If this should not be done or will not improve the solution, this should be set to False.

III.3.B  Bounds

Whereas most algorithms require input of an initial guess or starting point (which must be feasible) to initiate the search, GlobalSearch and GlobalPenaltyFn require only approximate (not strict) bounds on the variables.  Bounds are generally easier to define than are feasible starting points.  Physical, logical, or economic considerations often provide guidance on realistic bounds for a variable.  If the optimization problem is wing design and one variable is wing thickness, a certain range of thicknesses is sensible to consider.  For a financial problem, the range of values for an investment mix across a product line is clearly bounded between zero and the total budget, and the upper limit can generally be bounded below a much smaller value than this.  If narrower bounds are known, then convergence will require fewer function calls, but even quite wide bounds are acceptable.  If erroneously narrow bounds are input, then initial step size will be too small, and the function will act more like a local solver.  It is better to err on the side of wider bounds.  Upper and lower bounds must be input by the user, as illustrated below.  These bounds do not restrict the searching of GlobalSearch,n, IntervalMin, or GlobalPenaltyFn but merely provide an initial guide for generating feasible starts.  Initial bounds do not need to produce feasible solutions, the program can find initial feasible points.  If hard bounds on variables are necessary (such as positivity restrictions), they can be entered as constraints.

III.3.C  Constraints

Constraints are entered as Mathematica functions, which may be linear or nonlinear.  Equality constraints are entered in the third position in standard form.  For example, positivity restrictions on x and y would be represented by the list:

{-x,-y}

The following problem is a typical LP problem, with the solution at the intersection of two inequality constraints.  The solution is {2/3,10/3,0}.

GlobalSearch[-x1-2x2+x3,{x1+x2+x3-4,-x1+2x2-2x3-6,2x1+x2-5,-x1,-x2,-x3},,{{x1,0,1},{x2,0,1},{x3,0,1}},.000000001,Starts->1]//Timing

IntervalMin[-x1-2x2+x3,{x1+x2+x3-4,-x1+2x2-2x3-6,2x1+x2-5,-x1,-x2,-x3},,{Interval[{0.,4.}],Interval[{0.,4.}],Interval[{0.,1.}]},{x1,x2,x3},.0000001]//Timing

We see that IntervalMin is much slower, but there are problems where it performs better.  We next test GlobalPenaltyFn.

GlobalPenaltyFn[-x1-2x2+x3,{x1+x2+x3-4,-x1+2x2-2x3-6,2x1+x2-5,-x1,-x2,-x3},,{{x1,0,3},{x2,0,3},{x3,0,1}},.000001,Starts->2]//Timing

GlobalPenaltyFn can be slower than GlobalSearch because it can take nonanalytic constraints.

It is useful to illustrate a function in which the constraints are nonanalytic.  Any function for which Solve will not work falls in this category.  We take as the example:

Min[y^2+x]    s.t.   x==Sin[y]

GlobalSearch  and  GlobalPenaltyFn nevertheless succeed.

res=GlobalSearch[y^2+Sin[y],,,{{y,0,1}},.000000001,Starts->1]//Timing

Where x is then

x=N[Sin[y]]/.res[[2,1,1]]

ClearAll[x]

res=GlobalSearch[y^2+x,,{x-Sin[y]},{{x,0,1},{y,0,1}},.00000001,Starts->1,ShowProgress->False]//Timing

In this problem, GlobalPenaltyFn can solve it.

res=GlobalPenaltyFn[y^2+x,,{x-Sin[y]},{{x,-1,1},{y,0,1}},.000000001,Starts->1,ShowProgress->False]//Timing

In the next example, we see how a nonlinear constraint can be defined that creates a whole set of solution points, an approximation to which can be obtained with enough starts.  In this example, the solution must lie outside the circle of radius Sqrt[2], but the unconstrained solution is at the origin.  This means that all points on this circle are optimal.  This problem is a good candidate for parallel computing.  Some of the random starts do not give valid solution, and the program tries another point.  This example is run in parallel mode.

ClearAll[x]

b=GlobalPenaltyFn[x^2+y^2,{2-(x^2+y^2)},,{{x,-1,1},{y,-1,1}},.0001,Starts->100,ParallelOption->True];

b[[1]]

c={};Do[AppendTo[c,{x,y}/.b[[ii,1]]],{ii,1,100}]

ListPlot[c]

Equality constraints are input as a list of equality statements in the third parameter position.  In the following problem, the solution is {1/3,1/3,1/3}, with obj=1/3.

GlobalSearch[x^2+y^2+z^2,{},{x+y+z-1},{{x,0,1},{y,0,1},{z,0,1}},.00000001,Starts->1]//Timing

GlobalPenaltyFn[x^2+y^2+z^2,{},{x+y+z-1},{{x,0,1},{y,0,1},{z,0,1}},.00000001,Starts->1]//Timing

The option PenaltyMethod->True is useful for complex constraints.  The following problem needs these options.

vars={x1,x2,x3};nvars=Length[vars];varlist={{x1,0,1},{x2,0,1},{x3,0,1}};(* randomly generated solution to test problem *)sol=Table[Random[],{nvars}];

objf=(x1-sol[[1]])^2+2*(x2-sol[[2]])^2+3*(x3-sol[[3]])^2;

eqs={Sin[5*(x1-sol[[1]])]-12*(x3-sol[[3]])*(x2-sol[[2]]),Exp[(x1-sol[[1]])*(x2-sol[[2]])*(x3-sol[[3]])]-1};

ineqs={((x1-sol[[1]])*(x2-sol[[2]])*(x3-sol[[3]]))^2+Sin[15*(x1*x2*x3-sol[[1]]*sol[[2]]*sol[[3]])]-0.05};

sol

res=GlobalSearch[objf,ineqs,eqs,varlist,.0000000001,Starts->1,ShowProgress->False]//Timing

res=GlobalPenaltyFn[objf,ineqs,eqs,varlist,.00000000001,Starts->3,PenaltyMethod->True]//Timing

In some cases, the solution to a problem occurs at the intersection of constraints.  These problems are difficult but can be solved.  In the following problem, the solution is exactly at the intersection of the line and the circle:

r=Solve[{x1+x2-1.2==0,x1^2+x2^2==1}]

(x1-.1)^2+(x2-.1)^2+x3^2/.r[[1]]

In this problem, some of the solutions are not optimal.

GlobalSearch[(x1-.1)^2+(x2-.1)^2+x3^2,{x1+x2-1.2},{x1^2+x2^2-1},{{x1,0,1},{x2,0,1},{x3,5,10}},.000000001,Starts->1,ShowProgress->False]//Timing

We see that GlobalSearch found one of the solutions with 1 starts.  GlobalPenaltyFn can also find a solution but is slower.

GlobalPenaltyFn[(x1-.1)^2+(x2-.1)^2+x3^2,{-x1,-x2,x1+x2-1.2},{x1^2+x2^2-1},{{x1,0,1},{x2,0,1},{x3,5,10}},.0000000001,Starts->1]//Timing

In the following example, the constraint forces a solution outside the unit circle.  The solution is 0.737157 at {.707,.707}, which is found by GlobalSearch and GlobalPenaltyFn.  IntervalMin gets pretty close.

GlobalSearch[(x-.1)^2+(y-.1)^2,{-(x^2+y^2-1)},{},{{x,-4,5},{y,-4,5}},.000000001,Starts->1]//Timing

IntervalMin[(x-.1)^2+(y-.1)^2,{-(x^2+y^2-1)},,{Interval[{0.,1.}],Interval[{0.,1.}]},{x,y},.0001]//Timing

GlobalPenaltyFn[(x-.1)^2+(y-.1)^2,{-(x^2+y^2-1)},,{{x,0,1},{y,0,1}},.000000001,Starts->3,ShowProgress->False]//Timing

WARNING:  Constraints can cause difficulties.  Constraints can be mistakenly formulated such that no feasible space exists.  In this case, GlobalSearch will stop with an error message:

ClearAll[x,y,z]

GlobalSearch[x^2+y^2+z^2,{x-1,-x+2},{},{{x,0,1},{y,0,1},{z,0,1}},.00000001,Starts->1]//Timing

III.4  Output

Following solution, a list is returned that contains the solution.  For example, two points in a 2 parameter problem would give the list {{{1.1,1.2},-5},{{1.2,1.21},-4.9}}.  Intermediate output can be printed out by setting ShowProgress->True.

GlobalSearch[x^2+y^2,{-x},,{{x,-5,5},{y,-5,5}},.000000001,Starts->1,ShowProgress->True]

Following solution, a list is returned that contains the solution as a list of replacement rules with corresponding solutions.

III.5  Maximization

To find a maximum instead of a minimum, simply call GlobalSearch with the negative of the function, keeping in mind that a negative function value will be printed out.

III.6  Limitations

No optimization algorithm can solve all conceivable problems.  Consider a response surface in two dimensions that is flat everywhere except for a very small but very deep pit at one location.  Any hill-climbing algorithm will fail completely to solve this problem, because for almost all starting points there is no gradient.

Some functions have infinitely many solutions.  Consider Sin(1/x) between 0 and 1.  As the function approaches zero, the minima (with z = -1) get closer and closer together, as we can see from the plot following:

Plot [Sin[1/x], {x, 0, 1}]

Although it is fundamentally impossible to list all solutions when there are infinitely many, GlobalSolve can show that there are very many solutions.  If we run the program with 40 starts on the interval [0, 0.6], we obtain many points with values close to -1.0 (shown below).  This result demonstrates that there are many solutions and that the density of solutions increases with an approach to zero.  This is thus an approximate solution to the true situation.  x is bounded away from 0 to prevent underflows.

GlobalSearch[Sin[1/x],{-x+.00001,x-.6},,{{x,0.0000001,.6}},.000001,CompileOption->False,Starts->20]

By restricting the search region, we find more solutions:

GlobalSearch[Sin[1/x],{-x+.000000000001,x-.2},,{{x,0.000001,.2}},.00001,Starts->20]

These examples illustrate that not all problems can be solved.  Some are not solvable by any algorithm.  Some may be ill-posed.  However, for many problems that cause other algorithms to fail, GlobalSolve either succeeds, has a probability of succeeding, or provides a partial or approximate solution. We may thus say that the set of unsolvable problems is much smaller than it is for other solution methods.  This is particularly so with respect to the inclusion of constraints.

The optimization algorithm fundamentally assumes real valued functions.  It is possible for some parameter values that are tested to return values that are complex or otherwise not Real.  The program assumes that these are illegal values and treats them like values that fail a constraint test.  This means that it can solve a problem where values can be complex.  The problem below works even if the input range does not include positive numbers.

GlobalSearch[x^.5,{},{},{{x,-2,-1}},.000000000001,Starts->1]

III.7  Error Messages

An effort has been made to trap as many errors as possible.  Many user errors will prompt Mathematica error messages.  On input, only the order of parameters is evaluated in the passing function call.  If a valid value for tolerance (> 0.0) is not input, the following error message is printed, and the program stops:

	Error: tolerance must be > 0

If the upper bound of the variable definitions is smaller than the lower bound, the program will also stop with an error message.

If a parameter is omitted from the calling sequence, Mathematica cannot continue.  To indicate the problem, it echoes back the names for functions and values for constants in the function call, and then it stops.  This type of echoed-back output means that a format error exists in the calling sequence.  Some syntax errors are trapped within the Global Optimization package.

After input, the program compiles the user input function and the constraints, speeding execution.  If these functions are improperly formatted, this step may produce a compilation error message, which will terminate the run.  If a function is defined such that it does not return numeric values, this will cause error messages from Mathematica, followed by termination of GlobalSearch.

When constraints are used, there may be no feasible region, particularly if the problem has many constraints.  If this occurs, the following error message is printed, and the program terminates:

	Error: no valid initial points found

In this case, check the constraints.  Constraints can be in conflict if care is not taken.

A common mistake is to maximize when the intention was to minimize.  This error can be recognized because the solution will tend to run into one of the parameter bounds and will have an illogical value.  For example, the optimal airplane wing thickness for a design problem might come back as zero.

An unbounded solution may result from an improperly formulated problem (e.g., min(-1/x) over {-1, 1} which becomes -inf at x=0).  Because the function z continues to get larger the closer one gets to zero at an increasing rate, the program will never report convergence.  ShowProgress->True can be used to check for such problems during a run.

If the user defines variables in his program by mistake that should be parameters of the function, this will cause the program to stop or malfunction.  For example, if pp=5 is defined by mistake, this will cause GlobalSearch to detect an error which can be understood from looking at the error output which shows an invalid parameter list:

pp=5;

GlobalSearch[(pp-30)^2 (y-60)^2,,,{{pp,0.,100.},{y,1,5}},0.001]

Note how the program echoes back the parameter input, showing the mistake.

If the function to be solved is by mistake unbounded, then the parameter MaxIterations will stop the program from running forever.

GlobalSearch[x+y,{},{},{{x,-5,5},{y,-5,5}},.0000001,Starts->1,MaxIterations->4]

MaxIterations can also be used to make a preliminary run to test out execution.

r=GlobalSearch[100.*(x^2-y)^2+ (1.-x)^2,{},{},{{x,-5,5},{y,-10.,10.}},.0000001,Starts->1,MaxIterations->1]

We now wish to restart the program with the given ending value.

rr=r[[1]]

s={x,y}/.rr[[1]]

GlobalSearch[100.*(x^2-y)^2+ (1.-x)^2,{},{},{{x,-5,5},{y,-10.,10.}},.0000001,StartsList->{s},ShowProgress->False]//Timing

III.8  Performance/Speed

For larger problems, speed can become a limiting factor.  It is suggested that the objective function and constraints be formulated such that they can be Compiled (see sec. II) which can speed up operation by 2x to 20x.  If a function is very complex, it can be helpful to pull out subexpressions that do not need to be computed each time, and compute them once in advance.  Simplify or FullSimplify may also speed up the function.  Some large problems are run next to illustrate performance.  Either GlobalSearch or GlobalPenaltyFn can be run for large unconstrained problems, for which they use the same basic algorithm.

vars=Table[ToExpression[StringJoin["x",ToString[i]]],{i,1,100}];

fg=Evaluate[Sum[vars[[i]]^2,{i,1,100}]];

vardefs=Table[{vars[[i]],-10,11},{i,1,100}];

GlobalSearch[fg,,,vardefs,.00000001,Starts->1,ShowProgress->False]//Timing

GlobalPenaltyFn[fg,,,vardefs,.0000001,Starts->1,ShowProgress->False]//Timing

vars=Table[ToExpression[StringJoin["x",ToString[i]]],{i,1,400}];

fg=Evaluate[Sum[vars[[i]]^2,{i,1,400}]];

vardefs=Table[{vars[[i]],-10,11},{i,1,400}];

GlobalSearch[fg,,,vardefs,.0000001,Starts->1]//Timing

Note that for large problems with constraints or lots of data, execution will be slower.

III.9  Testing and Examples

In this section, the general nonlinear functions are tested.  The test criteria are the ability to obtain accurate answers, time required to execute, size of problem, ability to handle noisy functions, and ability to find multiple solutions.  In the first example, we test for the ability to solve to arbitrary precision:

GlobalSearch[(x-30)^2 +(y-60)^2,,,{{x,0.,10.},{y,1,5}},0.000000000000000001,Starts->1]

In this example, a square root function is used.  This function is undefined for x<0.  Min is 0 at 0.

f=x^.5

GlobalSearch finds the answer for any starting range.

GlobalSearch[f,,,{{x,-3,5}},.00000000001,Starts->1,ShowProgress->False]//Timing

NMinimize can not solve it even with when all options are tried.  The save is true for power functions.

NMinimize[f,{{x,-3,5}}]

NMinimize[f,x]

NMinimize[x^2.21,{{x,-3,5}}]

The following problem is the difficult Rosenbrock function, but it is easily solved.

GlobalSearch[100.*(x^2-y)^2+ (1.-x)^2,{},{},{{x,-5,5},{y,-10.,10.}},.00000000001,Starts->1,ShowProgress->False]//Timing

Of particular interest is how the function performs on highly irregular functions.  The following irregular function is tested:

f:=Abs[2(x-24)+(x-24)*Sin[x-24]]

Plot [f, {x,-20,60}]

GlobalSearch[f,,,{{x,-40,40}},.0000001,Starts->1]//Timing

NMinimize does not solve this properly, instead finding a local minimum.  It will work with particular values for the initial range.

NMinimize[f,{x}]

Problems can be solved that cause FindMinimum and NMinimize to fail.  In the following, any number below 3 has the minimum function value of 1.

Plot[Max[1,x-2],{x,0,4}]

GlobalSearch[Max[1,x-2],,,{{x,5,9}},.1,Starts->1,CompileOption->False]

FindMinimum[Max[1,x-2],{{x,5}}]

NMinimize[Max[1,x-2],{{x,5}}]

NMinimize[Max[1,x-2],{x}]

In the next problem, a more irregular function is tested:

Plot [4Sin[x]+Sin[4x]+Sin[8x]+Sin[16x]+Sin[32x]+Sin[64x], {x,0,6.28}]

GlobalSearch[4Sin[x]+Sin[4x]+Sin[8x]+Sin[16x]+Sin[32x]+Sin[64x],{-x,x-6},,{{x,0,6.28}},.000001,Starts->10]//Timing

NMinimize finds a local minimum.

NMinimize[{4Sin[x]+Sin[4x]+Sin[8x]+Sin[16x]+Sin[32x]+Sin[64x],{-x≤0,x-6≤0}},{{x,0,6.28}}]//Timing

All of the starts found the minimum (-6.98759).  We see next that IntervalMin fails on this problem because the way terms are handled separately does not allow the algorithm to narrow down the range.

IntervalMin[4Sin[x]+Sin[4x]+Sin[8x]+Sin[16x]+Sin[32x]+Sin[64x],{},,{Interval[{0.,6.28}]},{x},.000001]//Timing

In the next test, the Csendes function is tested, which has a highly dissected but flat region around the optimum, shaped like an ashtray (see next section for illustrative figures).

GlobalSearch[x^6 * (2.0 + Sin[1.0/x]) + y^6 * (2.0 + Sin[1.0/y]), ,,{{x, -1., 1.}, {y, -1., 1.}},0.0000000000000000001]//Timing

We can see that the minimum point {0,0} is being approached with arbitrary accuracy.

The following function has many local minima:

Plot3D[20Sin[Pi/2(x-2Pi)]+20Sin[Pi/2(y-2Pi)]+(x-2Pi)^2+(y-2Pi)^2,{x,0,10},{y,0,10}]

We can see that this could be difficult to find.  It is known that the solution is -38.0779 at {5.32216,5.32216}.

GlobalSearch[20Sin[Pi/2(x-2Pi)]+20Sin[Pi/2(y-2Pi)]+(x-2Pi)^2+(y-2Pi)^2,{},,{{x,0,10},{y,0,10}},.000000001,Starts->4]//Timing

GlobalSearch is able to find the true minimum when several starts are used.  IntervalMin solves this problem reliably:

IntervalMin[20Sin[Pi/2(x-2Pi)]+20Sin[Pi/2(y-2Pi)]+(x-2Pi)^2+(y-2Pi)^2,,,{Interval[{0.,10.}],Interval[{0.,10.}]},{x,y},.0001]//Timing

IV  NONLINEAR REGRESSION: THE NLRegression FUNCTION

IV.1  Introduction

NLRegression is a function that performs nonlinear regression using nonlinear least-squares.  It is able to incorporate constraints which may represent physical limits on the parameters being estimated, such as that a growth rate can not be negative.

The optimizer underlying NLRegression is GlobalSearch.  All options used by GlobalSearch are accepted by NLRegression and passed down to it except StartsList.  See the GlobalSearch section for its options.  The only other option accepted is Weights, which is the weight to be assigned to each data point.  Weighting might be employed when some of the data have a greater reliability than others.  The format for the function call is

NLRegression[data,expression,independent_variables,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

The output of the NLRegression procedure is a table of confidence intervals on the parameters, sensitivity plots, the fit parameter values, and the fit statistics.  The sensitivity plot can be very useful for detecting parameters that are redundant.

The typical regression problem involves fitting data to a nonlinear model.  If the model is nonanalytic (black box), the problem can not be solved by NLRegression, but instead must be solved by setting up the least-squares error function which is then passed in to GlobalSearch.  If the function returns values that are below machine accuracy or can return complex results, CompileOption->False should be used.  There are four options: FitStatistic, SensitivityPlots, Weights and Norm:

?FitStatistic

?Weights

?Norm

?SensitivityPlots

?UserResiduals

Let's consider a simple regression problem.

p=Table[i,{i,0,20}]

ClearAll[a,b,x]

eqn=a+b*x^2

Where a and b are parameters, and x is the independent variable.

data={};Do[AppendTo[data,{x,eqn/.{a->1,b->.1}}],{x,0,20}]

data

ListPlot[data]

r=NLRegression[data,eqn,{x},{},{},{{a,0,3},{b,0,1}},.0000001,SensitivityPlots->False,StartsList->{{2,.5}}]//Timing

We recover exactly our original parameters.  It is possible that the best fit to data violates some constraint such as physical feasibility.  In this case, it is necessary to constrain the optimization problem.  We introduce some measurement error into the data.

data[[1,2]]=data[[1,2]]-1

data[[2,2]]=data[[2,2]]-1

data

ListPlot[data]

r=NLRegression[data,eqn,{x},{},{},{{a,0,3},{b,0,1}},.0000001]//Timing

We see above that the error added to the data causes "a" to fall down to .792.  If we know that parameter a can not fall below .9 for some reason, we may add this constraint to the problem.  a>.9 is converted to standard form as {-a+.9}.  We add this to the function call.

r=NLRegression[data,eqn,{x},{-a+.9},{},{{a,0,3},{b,0,1}},.0000001]//Timing

The addition of the constraint satisfies the restrictions on the problem. Note that with constraints the confidence intervals and sensitivity plots are not defined and can not be printed. Another way to approach this problem is to use the L1 norm, which reduces the influence of outliers.  We try this next, without constraints.

r=NLRegression[data,eqn,{x},{},{},{{a,0,3},{b,0,1}},.0000001,Norm->L1]//Timing

We note two things immediately.  First, the sensitivity plots are now linear rather than quadratic, because of the different norm.  Second, the L1 norm recovers the true parameters in spite of the error in the data and without the use of constraints.

IV.2  Utilizing Chi-Square fit Criteria

It is sometimes useful to use a chi-square fit statistic:

FitStatistic -> ChiSquare

IV.3  Multiple Independent Variables in Regression Problems

In the next example,the fitting of a function to data is illustrated when there are multiple independent variables.  A common example of such a problem is when several experiments are done under different conditions and all data from all experiments are to be used to estimate certain parameters.

p=Table[i,{i,1,20}]

ClearAll[x,a,b,c,g]

y1=d(a*x^g+b)+c*(2*b+a*x^g)

c=1;d=0;
dat1={};curve1={};curve2={};Do[k=y1/.{x->p[[i]],a->.1,b->10.,g->2};AppendTo[dat1,{c,d,p[[i]],k}];AppendTo[curve1,{p[[i]],k}];,{i,1,20}];
c=0;d=1;
Do[k=y1/.{x->p[[i]],a->.1,b->10.,g->2};AppendTo[dat1,{c,d,p[[i]],k}];AppendTo[curve2,{p[[i]],k}];,{i,1,20}]

p1=ListPlot[curve1]

p2=ListPlot[curve2]

Show[p1,p2]

ClearAll[c,d,x]

indepvars={c,d,x}

r=NLRegression[dat1,y1,indepvars,{},{},{{a,0,1},{b,0,20},{g,1,3}},.0000001]//Timing

rr=r[[2]]

rrr=rr[[3]]

pars=rrr[[2]]

y1

pars2=Join[pars,{c->1,d->0}]

p3=Plot[y1/.pars2,{x,1,20}]

pars3=Join[pars,{c->0,d->1}]

p4=Plot[y1/.pars3,{x,1,20}]

Show[p1,p2,p3,p4]

We see that we have simultaneously fit both curves.  Problems with multiple independent variables can be fit even if they can not be graphed in the manner above.

V  MAXIMUM LIKELIHOOD ESTIMATION: THE MaxLikelihood FUNCTION

V.1  Introduction

In this section, the MaxLikelihood function is described.  This function is an estimation procedure which attempts to maximize the log-likelihood of a function with respect to a data set.  It is a useful alternative to nonlinear regression (the NLRegression function).  The statistics used are based on Crooke et al. (1999).  The function uses GlobalSearch as the minimizing function, which is needed because some functions require constraints to be solved properly.  It thus accepts all GlobalSearch options, except CompileOption and SimplifyOption.  This is because Compile does not work well with Log[], and because a sum of Log terms will generally not simplify.  It is assumed that the function to be solved is algebraic, and black box functions are not allowed.  Weighting is not yet implemented.  Special functions are pre-programmed for efficiency and ease of use.  The pre-programmed functions can be ascertained as follows:

?LikelihoodModels

The syntax for the function is as follows:

MaxLikelihood[data,expression,independent_variables,constraints,{{var1name,lowbound,highbound}..},tolerance,options]

V.2  Examples and Built-In Functions

We first illustrate use of the MaxLikelihood function for determining the mean μ  and standard deviation σ  from a random sample that has been generated by the Normal distribution:

			f(x;μ,σ)=((e^(-((((x+μ))^2)/(2(σ^2)))))/(σ(Sqrt[2π]))), x∈ℝ.
			
This probability distribution density is defined as NormalDistribution.

For demonstration purposes, we draw a set of pseudo random numbers from the Normal distribution with μ=10 and σ=1.

T=100;
exact={μ->10,σ->1};
data=Table[Random[NormalDistribution[μ,σ]/.exact],{T}];

Next, we define the probability density that we would like to fit as a pure function.

Clear[f];
f=PDF[NormalDistribution[μ,σ],#]&;

g=f[x]

mle=MaxLikelihood[data,g,{x},{},{{μ,8,12},{σ,.8,1.5}},.000001,CompileOption->True]//Timing

We see that we have recovered the true parameters pretty well.  A large sample size would be needed for a better estimate.  We can compare the results using the built-in function NormalModel:

mle=MaxLikelihood[data,NormalModel,{x},{},{{μ,8,12},{σ,.8,1.5}},.000001,CompileOption->True]//Timing

We see that the exact same result is obtained, but the run is about 4 times faster.  Note that the order of the parameters, but not their names, is assumed to match the {mean,variance} parameters of the normal model.  In this example, since the data is a sample from the distribution, we do not get back exactly the parameters we used to generate the data unless we draw a very large sample data set.

T=200;
exact={μ->10,σ->1};
data=Table[Random[NormalDistribution[μ,σ]/.exact],{T}];

mle=MaxLikelihood[data,NormalModel,{x},{},{{μ,8,12},{σ,.8,1.5}},.000000001]//Timing

We see above that a larger sample size gave a much better estimate of the parameters.  Note also that the execution time is nearly unaffected by the size of the data set for the built in model, although time required goes up linearly with data set size for the standard function.

Next we consider the estimation of the mean λ for a sample data drawn from the Poisson distribution:

		f(x;λ)=(((e^(-λ))(λ^x))/(x!)),x=0,1,2,… .
		
We use the same approach to illustrate the method, generating a sample with a known mean and then estimating the mean from the sample.  We draw sample data from a distribution with λ=10.

Clear[f];
T=100;
exact=λ->10;
data=Table[Random[PoissonDistribution[λ]/.exact],{T}];
f=PDF[PoissonDistribution[λ],#]&;

g=f[x]

mle=MaxLikelihood[data,PoissonModel,{x},{},{{λ,6,9}},.000001]//Timing

Once again, the preprogrammed model is faster.  Next we evaluate the Beta function.

Our next example is the Beta Distribution

		f(x;α,β)=((Γ(α+β)(((1-x))^(q-1))(x^(p-1)))/(Γ(α)Γ(β))),0≤x≤1
		
where α,β>-1 are parameters and Γ denotes the gamma function.   As above, we draw sample data from the specified distribution.

Clear[f];
T=200;
exact={α->2,β->5};
data=Table[Random[BetaDistribution[α,β]/.exact],{T}];
f=PDF[BetaDistribution[α,β],#]&;

g=f[x]

In this case, we must add constraints or the function goes to negative infinity during the solution.

mle=MaxLikelihood[data,g,{x},{-α+1,-β},{{α,1,5},{β,3,8}},.000001]//Timing

mle=MaxLikelihood[data,BetaModel,{x},{-α+1,-β+1},{{α,1,5},{β,3,8}},.000001]//Timing

Next we illustrate the Gamma distribution.  Note that the order of parameters is assumed to match that in the Gamma built in function.

Clear[f];
T=100;
exact={α->1,λ->1};
data=Table[Random[GammaDistribution[α,λ]/.exact],{T}];
f=PDF[GammaDistribution[α,λ],#]&;

g=f[x]

Note that we must constrain the parameters to be positive to prevent the function running away to -Infinity.  In addition, lambda must be constrained away from 0 to prevent underflow during computations, even though the final answer is not near zero.

mle=MaxLikelihood[data,g,{x},{-α+.0001,α-5,-λ+.0001,λ-8},{{α,.2,5},{λ,3,8}},.000001]//Timing

mle=MaxLikelihood[data,GammaModel,{x},{-α+.0001,α-5,-λ+.0001,λ-8},{{α,.2,5},{λ,3,8}},.000001]//Timing

In this case we get a huge speedup (10 times faster).  Next, we test the LogNormal predefined model.

Clear[f];
T=100;
exact={μ->20,σ->5};
data=Table[Random[LogNormalDistribution[μ,σ]/.exact],{T}];
f=PDF[LogNormalDistribution[μ,σ],#]&;

g=f[x]

mle=MaxLikelihood[data,g,{x},{},{{μ,8,12},{σ,.8,10}},.000001]//Timing

mle=MaxLikelihood[data,LogNormalModel,{x},{},{{μ,8,12},{σ,.8,10}},.000001]//Timing

VI  DISCRETE VARIABLE PROBLEMS: THE 
 InterchangeMethodMin & TabuSearchMin FUNCTIONS

VI.1  Introduction

InterchangeMethodMin and TabuSearchMin are functions which maximize a linear or nonlinear function of integer 0-1 variables either with or without constraints.  This type of problem comes up often in operations research (e.g., network, scheduling, and allocation applications).  Integer 0-1 problems are not generally solvable with an analytic approach such as Linear Programming.  The method used here is the Interchange method (Densham and Rushton, 1992; Goldberg and Paz, 1991; Lin and Kernighan, 1973; Teitz and Bart, 1968), which is similar to the author's SWAP algorithm (Loehle, 2000).  It has been shown that the Interchange method is guaranteed to get close to the optimum solution, but for complex problems may not find the exact solution.    On the other hand, it is much faster than Simulated Annealing (Murray and Church, 1995).  The Interchange algorithm begins with a feasible start generated randomly.  It then makes changes in the configuration at each iteration until no more progress is possible.  Note that while it is possible to input this type of problem into MultiStartMin because it accepts discrete variables, MultiStartMin assumes that only a few discrete variables are involved, and may generate a large set of possible moves and run slowly for large problems.  TabuSearchMin operates with a basic Interchange framework, but then adds a tabu list feature.  The tabu list is the list of n previous moves that will not be revisited in looking for the next good move.  This feature adds efficiency because bad moves are not checked repeatedly.  While GlobalSearch and GlobalPenaltyFn can accept Binary variables, they should not be used for large problems because memory requirements can become excessive.

The functions are defined by

InterchangeMethodMin[expr,ineqs,{varlist},tolerance,options]

TabuSearchMin[expr,ineqs,{varlist},tolerance,options]

where expr is the function to be minimized, which need not be linear, and ineqs is the inequalities, which are optional and are in standard form. There are three options for the program.  Starts defines the number of random starting points to test (default 5).  Input starting points can be input with the StartsList option.  CompileOption can be set to False if the function is not compilable.  ShowProgress shows intermediate stages of the computation.  Program operation is shown next in the following examples of applications.  Note that the range of applications is not exhausted by these examples.

VI.2  Applications

VI.2.A Capital Allocation

In many investment problems, the goal is to allocate a fixed amount of money across a series of investments.  When the items to be allocated are discrete, and must be allocated to one or another use in their entirety, then the problem must be solved with integer programming.  In the following example, the variables {x1,x2,x3,x4} are 0-1.  There are several constraints.  The objective function is given by

Max(z)=.2x1+.3x2+.5x3+.1x4

This problem may be solved as follows.  We minimize the negative of z to maximize, and enter all the constraints, which are less than inequalities:

ClearAll[x1,x2,x3,x4]

InterchangeMethodMin[-(.2x1+.3x2+.5x3+.1x4),{.5x1+x2+1.5x3+.1x4-3.1,.3x1+.8x2+1.5x3+.4x4-2.5,.2x1+.2x2+.3x3+.1x4-.4},{x1,x2,x3,x4},.001,Starts->4]//Timing

Seven random starts were used, but only 2 unique starts were used due to duplications.  Two starts found the correct answer {0,0,1,1}.  On larger problems, performance might not be quite this good.  TabuSearch also can solve this problem correctly:

TabuSearchMin[-(.2x1+.3x2+.5x3+.1x4),{.5x1+x2+1.5x3+.1x4-3.1,.3x1+.8x2+1.5x3+.4x4-2.5,.2x1+.2x2+.3x3+.1x4-.4},{x1,x2,x3,x4},.001,Starts->4,TabuListLength->1]//Timing

Note that Tabu Search may perform worse than Interchange on small problems.  Note also that a small tabu list size is needed for this problem.  If we specify a tabu list too long, a warning results:

TabuSearchMin[-(.2x1+.3x2+.5x3+.1x4),{.5x1+x2+1.5x3+.1x4-3.1,.3x1+.8x2+1.5x3+.4x4-2.5,.2x1+.2x2+.3x3+.1x4-.4},{x1,x2,x3,x4},.001,Starts->4,TabuListLength->100]//Timing

VI.2.B Vehicle Routing/Travelling Salesman

A very common problem in operations research is vehicle routing.  This problem arises in delivery service, airline, and other transportation contexts.  To illustrate, define a set of 8 cities with coordinates:

c={{1,1},{2,0.1},{3,0.1},{4,1},{1,2},{2,2},{3,2},{4,2}}

Where the objective is to minimize the sum of the links between cities.  We define the objective function as the sum of the distance*existence of the link over all possible links between cities.  Self-loops (e.g., x11) are removed from the lists.  Constraints are added that force every city to have at least 2 links (for entry and exit).  Finally, the sum of all edges must be at least 8.  We really want exactly 8 links, but the algorithm does not accept equality constraints.  The minimization of the sum of link distances will force the solution to exactly 8. Links are defined between cities such that there are no repeated links (e.g., x13 but no x31) to make it faster.

lis={x12,x13,x23,x14,x24,x34,x15,x25,x35,x45,x16,x26,x36,x46,x56,x17,x27,x37,x47,x57,x67,x18,x28,x38,x48,x58,x68,x78}

InterchangeMethodMin[2(1.3453624047073711` x12+2.193171219946131` x13+3 x14+x15+(Sqrt[2]) x16+(Sqrt[5]) x17+(Sqrt[10]) x18+1.` x23+2.193171219946131` x24+2.1470910553583886` x25+1.9` x26+2.1470910553583886` x27+2.758622844826744` x28+1.3453624047073711` x34+2.758622844826744` x35+2.1470910553583886` x36+1.9` x37+2.1470910553583886` x38+(Sqrt[10]) x45+(Sqrt[5]) x46+(Sqrt[2]) x47+x48+x56+2 x57+3 x58+x67+2 x68+x78)+Abs[x12+x13+x14+x15+x16+x17+x18-2]+Abs[x12+x23+x24+x25+x26+x27+x28-2]+Abs[x13+x23+x34+x35+x36+x37+x38-2]+Abs[x14+x24+x34+x45+x46+x47+x48-2]+Abs[x15+x25+x35+x45+x56+x57+x58-2]+Abs[x16+x26+x36+x46+x56+x67+x68-2]+Abs[x17+x27+x37+x47+x57+x67+x78-2]+Abs[x18+x28+x38+x48+x58+x68+x78-2],{-(x12+x13+x14+x15+x16+x17+x18+x23+x24+x25+x26+x27+x28+x34+x35+x36+x37+x38+x45+x46+x47+x48+x56+x57+x58+x67+x68+x78)+8},lis,.01,Starts->5,CompileOption->True,ShowProgress->False]//Timing

We see that the function found the solution in 4 out of 5 starts, which is a loop around the points in a rough circle.  Even the starts that failed got quite close to the best solution and had the proper number of links defined.  The timing is quite good.  The problem of subgraphs that are not linked to form a complete vehicle route can occur in this problem but is unlikely unless the subgraphs produce exactly the same total travel distance as the linked route.

TabuSearchMin[2(1.3453624047073711` x12+2.193171219946131` x13+3 x14+x15+(Sqrt[2]) x16+(Sqrt[5]) x17+(Sqrt[10]) x18+1.` x23+2.193171219946131` x24+2.1470910553583886` x25+1.9` x26+2.1470910553583886` x27+2.758622844826744` x28+1.3453624047073711` x34+2.758622844826744` x35+2.1470910553583886` x36+1.9` x37+2.1470910553583886` x38+(Sqrt[10]) x45+(Sqrt[5]) x46+(Sqrt[2]) x47+x48+x56+2 x57+3 x58+x67+2 x68+x78)+Abs[x12+x13+x14+x15+x16+x17+x18-2]+Abs[x12+x23+x24+x25+x26+x27+x28-2]+Abs[x13+x23+x34+x35+x36+x37+x38-2]+Abs[x14+x24+x34+x45+x46+x47+x48-2]+Abs[x15+x25+x35+x45+x56+x57+x58-2]+Abs[x16+x26+x36+x46+x56+x67+x68-2]+Abs[x17+x27+x37+x47+x57+x67+x78-2]+Abs[x18+x28+x38+x48+x58+x68+x78-2],{-(x12+x13+x14+x15+x16+x17+x18+x23+x24+x25+x26+x27+x28+x34+x35+x36+x37+x38+x45+x46+x47+x48+x56+x57+x58+x67+x68+x78)+8},lis,.01,Starts->5,TabuListLength->60]//Timing

We see that the TabuSearchMin function found the solution in 1 out of 5 starts, which is a loop around the points in a rough circle, but in 1/3 the time of the InterchangeMethod.  Even the starts that failed got quite close to the best solution and had the proper number of links defined.  The timing is quite good.  The problem of subgraphs that are not linked to form a complete vehicle route can occur in this problem but is unlikely unless the subgraphs produce exactly the same total travel distance as the linked route.

VI.2.C Minimum Spanning Tree

The construction of a minimum spanning tree is a problem that comes up in problems of delivery from fixed warehouses, utility network design, and other applications.  We begin with the weighted graph, with nodes (a,b,c,d,e,f), and with edges and weights defined by:

{xab,xac,xae,xbc,xbd,xcd,xce,xcf,xdf,xef}

{7,2,2,1,3,3,4,4,5,6}

This can be solved with the interchange method.  The objective function is to minimize the sum of the weighted edges.  This produces the minimum tree we desire.  We know that we want 5 edges as the final solution (#nodes -1), so we add an inequality that the sum of the edges must be greater than 5.  This allows more than 5 initially, but then edges are trimmed as the algorithm proceeds.  If an equality constraint is put in, the program will fail.  Additional constraints force every node to have at least one edge touching it.

InterchangeMethodMin[7xab+2xac+2xae+1xbc+3xbd+3xcd+4xce+4xcf+5xdf+6xef,{-(xab+xac+xae+xbc+xbd+xcd+xce+xcf+xdf+xef)+5,-(xab+xac+xae)+1,-(xab+xbc+xbd)+1,-(xac+xbc+xcd+xce+xcf)+1,-(xbd+xcd+xdf)+1,-(xae+xce+xef)+1,-(xcf+xdf+xef)+1},{xab,xac,xae,xbc,xbd,xcd,xce,xcf,xdf,xef},.1,Starts->1,CompileOption->True]//Timing

We see that the function returns the correct solution {ac,ae,bc,bd,cf} with a very fast execution time.  This reflects the fast computational time for binary variables.  Large problems can thus be solved.

TabuSearchMin[7xab+2xac+2xae+1xbc+3xbd+3xcd+4xce+4xcf+5xdf+6xef,{-(xab+xac+xae+xbc+xbd+xcd+xce+xcf+xdf+xef)+5,-(xab+xac+xae)+1,-(xab+xbc+xbd)+1,-(xac+xbc+xcd+xce+xcf)+1,-(xbd+xcd+xdf)+1,-(xae+xce+xef)+1,-(xcf+xdf+xef)+1},{xab,xac,xae,xbc,xbd,xcd,xce,xcf,xdf,xef},.1,Starts->1,TabuListLength->10,CompileOption->True]//Timing

We see that the function returns the correct solution {ac,ae,bc,bd,cf} with a very fast execution time. For such smaller problems, the Tabu Search approach does not provide superior speed.

VII  THE MaxAllocation FUNCTION

VII.1  Introduction

MaxAllocation is a function to maximize a nonlinear function subject to a linear constraint.

Max(f((Subscript[x, i]))) s.t.0<(Subscript[x, i])<rhs, (Underoverscript[∑(Subscript[x, i]), (i=1), n])=rhs

This type of problem comes up often in finance, economics, and investment.  Traditional algorithms have a great deal of trouble with such a problem, though apparently simple.  A new algorithm enables a very efficient solution of this type of problem.  The algorithm is based on the idea of path following.  If one can stay on the constraint line, then the solution stays feasible at all times.  An approach to path following is based on the SWAP algorithm (Loehle, 2001). We first discretize the interval of the rhs.  This discretization is made proportional to the number of variables, N, for convenience.  The initial 2N pieces are placed onto all the variables equally (each variable is given a value 2*rhs/N).  At iteration 1, each variable is tested to see if swapping a single piece from that variable would improve the solution.  This performs a rectangular move in N space.  After no more pieces can be moved from variable 1 (either all have been moved or no improvement is obtained), testing proceeds to variable 2 to see if any pieces can be moved, and so on.  Note that some of the pieces moved from 1 to 2 might be redistributed at this point.  After the algorithm has tested all the variables in order for the possibility of making a swap, the number of pieces on each variable is doubled (with half the length) for the next iteration.  This enables a finer resolution of the function.

The function is defined by

MaxAllocation[function,varlist,rhs,tolerance,options]

where rhs is the right hand side value (from eqn (1)) and tolerance defines the stopping criterion.  There are three options for the program.  MaxIterations will stop the program at the specified number of iterations if desired.  CompileOption is useful if the function is not compilable.  ShowProgress shows intermediate stages of the computation.

Operation of the function is shown in the following example.  This example is a sum of diminishing return terms.  This is a function with no nonlinear interaction terms.  It can be shown that the algorithm converges to the analytic optimum in this case.

ClearAll[x2]

varlist={x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32,x33,x34,x35,x36,x37,x38,x39,x40,x41,x42,x43,x44,x45,x46,x47,x48,x49,x50,x51,x52,x53,x54,x55,x56,x57,x58,x59,x60};

v[x_]:=Sum[i*(1.0-E^(-3.*x[[i]])),{i,1,60}]

Note the technique used here to take advantage of indexed variables.  The function v is defined with a tensor as an argument.  Then function vvv is defined to accept a list and turn it into a tensor by passing it down to v.  This function can be compiled, and is in the following example.

vvv=Function[{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32,x33,x34,x35,x36,x37,x38,x39,x40,x41,x42,x43,x44,x45,x46,x47,x48,x49,x50,x51,x52,x53,x54,x55,x56,x57,x58,x59,x60},v[{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32,x33,x34,x35,x36,x37,x38,x39,x40,x41,x42,x43,x44,x45,x46,x47,x48,x49,x50,x51,x52,x53,x54,x55,x56,x57,x58,x59,x60}]];

t=Table[1,{i,1,60}]

v[t]

Apply[vvv,t]

MaxAllocation[vvv[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32,x33,x34,x35,x36,x37,x38,x39,x40,x41,x42,x43,x44,x45,x46,x47,x48,x49,x50,x51,x52,x53,x54,x55,x56,x57,x58,x59,x60],varlist,1.0,.1,CompileOption->True,ShowProgress->False]//Timing

Next, we can compare the value of compiling the function in terms of timing.

MaxAllocation[vvv[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32,x33,x34,x35,x36,x37,x38,x39,x40,x41,x42,x43,x44,x45,x46,x47,x48,x49,x50,x51,x52,x53,x54,x55,x56,x57,x58,x59,x60],varlist,1.0,.1,CompileOption->False]//Timing

We can see that the compiled version runs 5 times faster.

It would be useful to have a proof that this algorithm works. At this time, no proof is available.  The algorithm has been tested on a suite of problems with known solutions, and gives good results.   The algorithm has been designed to guard against infinite loops.   It is possible that a problem could be devised with local minima that would prevent a global optimum from being reached, but no such problem has yet been encountered. The initial equitable distribution of pieces across the variables guards against interaction (e.g., x1*x2) type terms which can not be improved by rectangular moves if they are initialized to zero.

The above problem with 250 variables took 10 minutes to run on a Pentium III 600 MHz machine.  Because the execution time goes up nonlinearly with the number of variables, it is likely that the upper limit for an overnight run is 350-1000 variables, depending on the machine used.

VII.2  Applications

VII.2.A Investment Allocation Problems

In many investment problems, the goal is to allocate a fixed amount of money across a series of investments.  We need not invest in all the options available.  A typical problem to solve is the following, (from A.K. Dixit, Optimization in Economic Theory):  A capital sum C is available for allocation among n investment projects.  The expected return from a portfolio of (Subscript[x, j]) projects is

(∑_(j=1)^n)[(Subscript[α, j])(Subscript[x, j])-0.5(Subscript[β, j])(Subscript[x, j]^2)]

subject to:

(∑_(j=1)^(Subscript[nx, j]))=C

This problem may be solved as follows:

a={1000,2000,3000};

rhs=10.0;

b={100,500,1000.};

h[x_]:=Sum[a[[i]]*x[[i]]-.5*b[[i]]*x[[i]]^2,{i,1,3}];

fff=Function[{a1,a2,a3},h[{a1,a2,a3}]];

MaxAllocation[fff[a1,a2,a3],{a1,a2,a3},rhs,.1,ShowProgress->True]

Here we find the solution which is an interesting mix of investments.  This problem can be extended to 350 to 1000 variable cases for overnight runs, depending on machine speed.

VII.2.B Derivative Hedging with Quadratic Programming

A common problem in finance is derivative construction to hedge an investment.  One technique for derivative construction is to construct a target payoff value (T) for the instrument, and then to allocate investment to a series of instruments with return R to get as close to the payoff value as possible using a quadratic objective function (for 3 instruments here):

(((∑_(i=1)^3)[(Subscript[T, i])-(∑_(j=1)^3)(Subscript[x, j])(Subscript[R, i,j])])^2)

(∑_(j=1)^3)(Subscript[x, j])(Subscript[p, j])≤C

This can be solved with MaxAllocation.  First, for a problem with 3 possible investments, we make a dummy variable (index 4 here) to take up the slack in the investment (to represent the part of C NOT invested).  Also, the equation is rearranged so that x represents the dollar amount invested in each instrument, so x is divided by the price p in eqn. (3).  Then C represents our rhs and the problem is given by:

p={1,1,1};

targetpayoff={1000,2000,3000};

payoff={{1000,0,0},{0,2000,0},{0,0,3000}};

fy[x_]:=-Sum[(targetpayoff[[i]]-Sum[(x[[j]]/p[[j]])*payoff[[i,j]],
						{j,1,3}])^2,{i,1,3}];

fff=Function[{a1,a2,a3,a4},fy[{a1,a2,a3}]];

MaxAllocation[fff[a1,a2,a3,a4],{a1,a2,a3,a4},6.,.1]

In this example, we know the exact answer, and the algorithm converged to the exact answer quickly.  Such derivative construction problems can be solved with hundreds of variables.

VII  APPLICATIONS

VII.1  Zeros of a Function/Roots of a Polynomial

Although the Mathematica function RootFind can find the zeros of a function, it can occur that a function is not analytic, in which case RootFind will not work.  For example, a function might need to be calculated recursively, algorithmically, or by simulation.  FindRoot also will only find a single root at a time from a given starting point.  In addition, we may invoke the need for identifying the indifference zone: we may want all the solutions that are within some indifference region of the zero or root.  It is a simple matter to use either GlobalMinima or MultiStartMin to find roots of a function.  Consider the function

(x^2)-25

where the roots or zeros are obviously +5  and -5.  As formulated, the function looks like

Plot[x^2-25, {x, -10,10}]

and the minimum is at -25.  If, however, we evaluate the absolute value of the function

Plot[Abs[x^2-25], {x, -10,10}]

then we can simply minimize the function as in all the above examples and we will find the roots or zeros of the function, which are now the minima.  In cases where there is a region of solutions with values near zero, this procedure will likewise allow us to define these regions.

GlobalSearch[Abs[x^2-25],{},{},{ {x, -10,10}},.0000001,Starts->4]//Timing

We can see that all three functions find the two solutions but on this small problem GlobalMinima is faster.  The other 2 functions also need multiple starts to find both solutions.

VIII.2  Integer Programming: the Knapsack Problem

As discussed above, MultiStartMin allows integer variables.  This implies that it can solve integer programming problems.  Optimization problems with integer variables are difficult because they violate the assumptions of continuous methods such as LP or gradient descent.  Various methods have been applied to integer programming problems, including heuristic search.  The method used here for integer variables is a generalized descent with discrete step sizes, combined with limited interchange.  This method is illustrated for the knapsack problem.  In the knapsack problem, the optimum return on packing of discrete variables is desired.  This could occur when packing a truck, where different packages have different shipping values.  On the other hand, the size of each package may differ, and the total space available is constrained.  We next test such a problem.  Because we wish to maximize, we put a negative sign in front of the objective function.  To use GlobalPenaltyFn on this problem, it is necessary to multiply 100 times the positivity constraints to prevent it from going negative. Only the best solutions are shown in the output.  The GlobalPenaltyFn solution can be quite close.

GlobalPenaltyFn[-(x1+2x2+4x3+15x4+80x5+100x6),{-100x1,-100x2,-100x3,-100x4,-100x5,-100x6,1x1+3x2+15x3+29x4+70x5+59x6-300},,{{x1,0,5,Integer},{x2,0,5,Integer},{x3,0,5,Integer},{x4,0,5,Integer},{x5,0,5,Integer},{x6,0,5,Integer}},.00001,ShowProgress->False,Starts->4]//Timing

We see that the correct solution {{5,0,0,0,0,5},-505.} is found.

NMinimize[{-(x1+2x2+4x3+15x4+80x5+100x6),-100x1≤0&&-100x2≤0&&-100x3≤0&&-100x4≤0&&-100x5≤0&&-100x6≤0&&1x1+3x2+15x3+29x4+70x5+59x6-300≤0&&{x1,x2,x3,x4,x5,x6}∈Integers},{{x1,0,5},{x2,0,5},{x3,0,5},{x4,0,5},{x5,0,5},{x6,0,5}}]

VIII.3  Differential Equation Models

The following example shows how to fit parameters in a differential equation model.

sol=NDSolve[{y'[x]==y[x],y[0]==1},y,{x,0,2}]

y[1.5]/.sol

Plot[y[t]/.sol,{t,0,2}]

fn=Function[{a},so=NDSolve[{y'[x]==y[x],y[0]==a},y,{x,0,2}];j=y[1.5]/.so;j[[1]]]

fn[1]

gg[x_?NumericQ]:=Abs[fn[x]-fn[1]]

GlobalSearch[gg[x],{},,{{x,0,2}},.00001,CompileOption->False,SimplifyOption->False,EvaluateObj->False,Starts->1]//Timing

NMinimize fails

NMinimize[gg[x],{{x,0,2}}]

In this example we see that we recovered the original parameter a=1.  Note that CompileOption->False was necessary because a Mathematica function was passed in, which prohibits compilation.  SimplifyOption->False was also needed to prevent error messages (although the correct answer is still obtained without it).  It is also possible to fit models in cases where we want to fit the entire time trajectory to the data, in which case the function is evaluated at many points and some measure of fit computed.

IX.  Literature Cited

Cadzow, J. A. and Martens, H.R., 1970. Discrete-Time and Computer Control Systems, Prentice-Hall, Inc.

Courrieu, P.  1997.  The hyperbell algorithm for global optimization: A random walk using Cauchy densities.  J. Global Optimization 10:37-55.

Crooke, P., L. Froeb, and S. Tschantz.  1999.  Maximum likelihood estimation.  Mathematica in Education and Research 8:17-23.

Csendes, T.  1985.  A simple but hard-to-solve global optimization test problem.  IIASA Workshop on Global Optimization, Sopron, Hungary.

Cvijovic, D. and J. Klinowski.  1995.  Taboo search: an approach to the multiple minima problem.  Science 267:664-666.

Densham, P., and G. Rushton.  1992.  A more efficient heuristic for solving large p-nedian problems.  Pap. Reg. Sci. 71:307-329.

Dixon, L.C.W., and G.P. Szego.  1978.  The global optimization problem: An introduction.  Towards Global Optimization 2, North-Holland, Amsterdam.

Goldberg, J. and L. Paz.  1991.  Locating emergency vehicle bases when service time depends on call location.  Transportation Science  25:264-280.

Griewank, A.O.,  19 ___.  Generalized descent for global optimization.

Ingber, L., and B. Rosen.  1992.  Genetic algorithms and very fast simulated annealing: A comparison.  Mathematical and Computer Modelling 16:87-100.

Jones, D.R., C.D. Perttunen, and B.E. Stuckman.  1993.  Lipschitzian optimization without the Lipschitz constant.  J. Optimization Theory and Applications 79:157-181.

Lin, S. and B. Kernighan.  1973.  An effective heuristic algorithm for the traveling salesman problem.  Operations Research 21:498-516.

Loehle, C. 2000.  Optimal control of spatially distributed process models.  Ecological Modelling 131:79-95. 

Murray, A.T. and R.L. Church.  1995.  Heuristic solution approaches to operational forest planning problems.  OR Spektrum 17:193-203.

Nelson, J. and J.D. Brodie.  1990.  Comparison of random search algorithm and mixed integer programming for solving area-based forest plans.  Canadian J. Forest Research 20:934-942.

Pinter, J.D.  1996.  Global Optimization in Action.  Kluwer Academic Publishers, Boston.

Phillips, C. L. and Nagle, H. T., 1995. Digital Control System Analysis and Design, 3rd edition, Prentice-Hall, Inc.

Shahian, B. and Hassul, M., 1993. Control System Design Using MATLAB, Prentice-Hall, Inc.

Teitz, M. and P. Bart.  1968.  Heuristic methods for estimating the generalized vertex median of a weighted graph.  Operations Research 16:955-961.