An Overview of the GlobSol Package for Verified Global

Optimization

by R. Baker Kearfott Department of Mathematics

University of Louisiana at Lafayette

        rbk@louisiana.edu

This talk will

• Contrast verified and non-verified global optimization

• Briefly discuss the underlying mathematics of verified global optimization

• Review capabilities of the GlobSol software package

• Review an example of how to use GlobSol

• Give an on-line demonstration.

GlobSol July, 2000 Dresden–1

Local Versus Global Optimization

[figure: x_{1}^{l}; x_{2}^{l}; x_{∗}^{l}; x_{1}^{u}; x_{2}^{u}; local minima; global minimum; critical points]

Local Optimization Versus Global Optimization

Local Optimization

• The model is steepest descent with univariate line searches (for monotone decrease of the objective function). (Start a ball on a hill and let it roll to the bottom of the nearest valley.)

• Algorithm developers speak of “globalization,” but mean only the design of algorithm variants that increase the domain of convergence. (See J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Least Squares, Prentice–Hall, 1983.)

• Algorithms contain many heuristics, and do not always work. However, many useful implementations exist.

GlobSol July, 2000 Dresden–3

Local Optimization Versus Global Optimization

Global Optimization

• is a much harder problem. Progress has accelerated with increases in computing power.

• Early milestones are L. C. W. Dixon and G. P. Szego, Towards Global Optimization (North–Holland, 1975), and Towards Global Optimization 2 (North–Holland, 1977).

• Two types of algorithms: \underline{stochastic} and \underline{deterministic}.

• Deterministic algorithms can be either \underline{rigorous} or \underline{heuristic}.

GlobSol July, 2000 Dresden–4

Global Optimization

Stochastic Algorithms

Monte-Carlo search: Random points are generated in the search space. The point with lowest objective value is taken to be the global optimum.

Simulated annealing: is similar to a local optimization method, although larger objective values are accepted with a probability that decreases as the algorithm progresses.

Genetic algorithms: Attributes, such as values of a particular coordinate, correspond to particular “genes.” “Chromosomes” of these genes are recombined randomly, and only the best results are kept. Random “mutations” are introduced.

Deterministic Optimization

• involves some kind of systematic global search over the domain.

• The various algorithms rely on estimates of the range of the objective function over subdomains.

• Some algorithms (due to Mladineo, Schubert, Wood, etc.) rely on \underline{Lipschitz constants} to obtain estimates of ranges.

• Bounds on ranges are also obtained with \underline{outwardly rounded interval arithmetic}.

Hybrid Stochastic / Deterministic Algorithms

Recently, several people have combined statistical methods with deterministic methods.

• Janos Pinter uses a statistical model to estimate local approximations to Lipschitz constants for a global search.

• Donald Jones constructs a cumulative statistical model of of the objective, and uses deterministic global optimization with a simpler objective to determine optimal placement of the next sample point.

• Kaj Madsen uses multiple starts of a local optimizer to simulate a rigorous global search with interval methods.

On the State of the Art

• Minimizing a function over a compact set in ℝ^{n} is an NP-complete problem.

• Thus, barring monumental discoveries, any general algorithm will fail for some high-dimensional problems.

• There are many practical problems that can be solved in low-dimensional spaces.

• Some low-dimensional problems are difficult.

• Advances in computer speed and algorithm construction have allowed many more practical problems to be solved, including high-dimensional ones.

Deterministic Global Optimization

Interval Methods

• Evaluation of an objective function ϕ(X) at an interval vector X gives bounds on the actual range of ϕ over X.

– If \underline{directed rounding} is used, the bounds rigorously contain the mathematical range.

– The bounds, in general, are overestimates.

• If the lower bound of ϕ(X) is greater than a previously computed objective value ϕ(X), then X can be discarded.

• \underline{Interval Newton Methods}, combined with directed rounding, can prove existence and uniqueness of critical points, as well as reduce the size of regions X.

Interval Methods

Advantages

easier to use: Obtaining bounds with interval methods involves programming the objective function, while using Lipschitz constant-based methods may require extensive preliminary analysis.

more efficient: Despite interval overestimation of ranges, the overestimation is often less than with a fixed Lipschitz constant. (But keep in mind the success of hybrid deterministic / stochastic algorithms.)

more capable: With directed roundings, interval methods \underline{cannot lie}. Also, interval Newton iteration results in \underline{quadratic convergence effects}.

Underlying Mathematics

• Classical fixed point theory implies existence.

– Contraction Mapping Theorem

– Brouwer Fixed Point Theorem

– Miranda’s Theorem

• Regularity (non-singularity) implies uniqueness.

• Fundamental property of interval arithmetic allows computational existence and uniqueness.

Miranda’s Theorem

Theorem 1 Suppose X ∈ 𝕀ℝ^{n}, and let the faces of X be denoted by

X\underline{i} = (x_{1}, …, x_{i−1}, \underline{x}_{i}, x_{i+1}, …, x_{n})^{T}
X\overline{i} = (x_{1}, …, x_{i−1}, \overline{x}_{i}, x_{i+1}, …, x_{n})^{T} .

Let F = (f_{1}, …, f_{n})^{T} be a continuous function defined on X. If

f_{i}^{u}(X\underline{i})f_{i}^{u}(X\overline{i}) ≤ 0  (1)

for each i between 1 and n, then there is an X ∈ X such that F (X) = 0.

Regularity

Lemma 2 Suppose F : X → ℝ^{n} and A is a Lipschitz matrix, such as F′(X). If A is regular, then any root of F in X is unique.

Proof: Suppose X^{∗} ∈ X and X ∈ X have F (X^{∗}) = 0 and F (X) = 0. If A is a Lipschitz set, then there is an A ∈ A such that

F (X^{∗}) − F (X) = 0
= A(X^{∗}) − A(X)
= A(X^{∗} − X).

If X^{∗} ≠ X, then A would have a null vector, contradicting the regularity of A.

A Computational Existence Tool

Interval Gauss–Seidel Method

Assume F (X) = (f_{1}(X), …, f_{n}(X)) is continuously differentiable. Then the mean value theorem gives

f_{i}(x_{i}) = f_{i}(x_{¬i}, x̌_{i}) + \frac{∂f_{i}}{∂x_{i}}(x_{¬i}, ξ_{i})(x_{i} − x̌_{i}),

whence f_{i}(X) = 0 provided

x_{i} = x̌_{i} − f_{i}(x_{¬i}, x̌_{i})/\frac{∂f_{i}}{∂x_{i}}(x_{¬i}, ξ_{i})
⊆ x̌_{i} − f_{i}(x_{¬i}, x̌_{i})/\frac{∂f_{i}}{∂x_{i}}(x_{¬i}, x_{i})

The above inclusion forms the Interval Gauss–Seidel iteration equation.

Interval Gauss–Seidel Method

Verification Properties

Do sequentially for 1 = 1 to n:

1. x̃_{i} ← x̌_{i} − f_{i}(x_{¬i}, x̌_{i})/\frac{∂f_{i}}{∂x_{i}}(x_{¬i}, x_{i})

2. x_{i} ← x̃_{i}

Then:

1. If x̃_{i} ⊆ x_{i} after each step 1, then the hypotheses of Miranda’s theorem are satisfied, and, hence, there exists a solution to F (X) = 0 within X.

2. If x̃_{i} ⊆ x_{i} after each step 1, then F′(X) must be regular, and hence, the solution in X is also unique.

Note: The system F (X) = 0 is usually preconditioned by a point matrix Y to make the Jacobi matrix Y F′(X) approximately diagonal.

GlobSol July, 2000 Dresden–15

Preconditioned Interval Gauss–Seidel

An Example

f_{1}(x_{1}, x_{2}) = x_{1}^{2} − 4x_{2},
f_{2}(x_{1}, x_{2}) = x_{2}^{2} − 2x_{1} + 4x_{2},

and X = (x_{1}, x_{2})^{T} = ([−0.1, 0.1], [−0.1, 0.3])^{T}. X^{∗} = (0, 0)^{T}, of F (X^{∗} = 0 is unique in X. Take X̌ = (0, .1)^{T}, so F (X̌) = (−.4, .41)^{T}, and

F′(X) = ( \begin{matrix}2x_{1} & −4 \\ −2 & 2x_{2} + 4 \end{matrix} )
= ( \begin{matrix}[−.2, .2] & −4 \\ −2 & [3.8, 4.6] \end{matrix} ) .

An Example

(Continued)

Take Y to be the inverse of the midpoint matrix of F′(X):

Y = \{m(F′(X))\}^{−1} = ( \begin{matrix}−.525 & −.5 \\ −.25 & 0 \end{matrix} ) .

Then Y F′(X) = ( \begin{matrix}[0.895, 1.105] & [−.2, .2] \\ [−.05, .05] & 1 \end{matrix} ),

Y F (X̌) = (.005, .1)^{T}, and

x̃_{1} = 0 − \frac{−.005 + [−.2, .2][−.2, .2]}{[.895, 1.105]}
⊆ \frac{[−.035, .045]}{[.895, 1.105]} ⊆ [−.0392, .0503]
⊂ x_{1} = [−.1, .1].

Similarly, x̃_{2} ⊂ x_{2}, so there is a unique root of F in X.

What is GlobSol?

• A Fortran 90 package – well-tested.

– self-contained.

• Solves constrained and unconstrained global optimization problems

• Separate program solves square algebraic systems of equations.

• Utility programs for interval and point evaluation, etc.

• Subroutine / module libraries for interval arithmetic, automatic differentiation, etc.

• Publicly available free of charge

  http://interval.louisiana.edu/GlobSol/download GlobSol.html

GlobSol

Special Features

• Objective function and constraints are coded as Fortran 90 programs.

• Can use \underline{constraint propagation} (substitution/iteration) on the intermediate quantities in objective function, equality, and inequality constraint evaluation.

• Can use an overestimation-reducing “peeling” process for bound-constraints.

• Uses an effective point method to find approximate feasible points.

• Has a special augmented system mode for least squares problems.

GlobSol Features

(continued)

• Has extensive error-checking (user input, internal errors, etc.)

• Has on-line web page documentation.

• The algorithm is configurable.

• Has various levels of printing, for various algorithm aspects.

• Source code and libraries for components are available.

– Automatic differentiation access.

– Interval arithmetic access.

– User-modifiable, with adequate study.

• Gives performance statistics, both in report form and for input to spreadsheets.

Use of GlobSol

An Example

The following Fortran 90 program defines the objective function

minimize ϕ(X) = −2 ∗ x_{1}^{2} − x_{2}^{2} subject to constraints

x_{1}^{2} + x_{2}^{2} − 1 ≤ 0
x_{1}^{2} − x_{2} ≤ 0
x_{1}^{2} − x_{2}^{2} = 0

Use of GlobSol

An Example, continued

PROGRAM SIMPLE_MIXED_CONSTRAINTS
USE CODELIST_CREATION
PARAMETER (NN = 2)
TYPE(CDLVAR), DIMENSION(NN) :: X
TYPE(CDLLHS), DIMENSION(1):: PHI
TYPE(CDLINEQ), DIMENSION(2) :: G
TYPE(CDLEQ), DIMENSION(1)   :: C
CALL INITIALIZE_CODELIST(X)
PHI(1) = -2*X(1)**2 - X(2)**2
G(1) = X(1)**2 + X(2)**2 - 1
G(2) = X(1)**2 - X(2)
C(1) = X(1)**2 - X(2)**2
CALL FINISH_CODELIST
END PROGRAM SIMPLE_MIXED_CONSTRAINTS

GlobSol Example

(continued)

1. Running the above program produces an internal representation, or \underline{code list}.

2. The code list is then symbolically differentiated.

3. The optimization code interprets the derivative code list at run time to produce floating point and interval evaluations of the objective function, gradient, and Hessian matrix.

4. A separate data file defines the initial search box, the bound constraints, and the initial guess, if any.

5. A configuration file GlobSol.CFG supplies algorithm options, such as which interval Newton method to use and how to precondition the linear systems.

GlobSol July, 2000 Dresden–23

The Data File

1D-5       ! General domain tolerance
 0   1     ! Bounds on the first variable
 0   1     ! Bounds on the second variable
F F        ! X(1) has no bound constraints
F F        ! X(2) has no bound constraints

Subsequent optional lines can give an initial guess point.

Output File – abridged first part

Output from FIND_GLOBAL_MIN on  04/06/1999  at  08:03:52.
Version for the system is: March 20, 1999
Codelist file name is: MIXEDG.CDL
Box data file name is: MIXED.DT1
Initial box:
[    0.0000E+00,   0.1000E+01 ] [    0.0000E+00,   0.1000E+01 ]
BOUND_CONSTRAINT:
  F F   F F
---------------------------------------
CONFIGURATION VALUES:
EPS_DOMAIN:   0.1000E-04   MAXITR:   60000
DO_INTERVAL_NEWTON: T  QUADRATIC: T  FULL_SPACE: F
VERY_GOOD_INITIAL_GUESS:F
USE_SUBSIT:T
OUTPUT UNIT:7 PRINT_LENGTH:1
Default point optimizer was used.

Output File – abridged second part

THERE WERE NO BOXES IN COMPLETED_LIST.
LIST OF BOXES CONTAINING VERIFIED FEASIBLE POINTS:
Box no.:1
Box coordinates:
[    0.7071E+00,   0.7071E+00 ] [    0.7071E+00,   0.7071E+00 ]
PHI:
[   -0.1500E+01,  -0.1500E+01 ]
Level: 3
Box contains the following approximate root:
   0.7071E+00   0.7071E+00
OBJECTIVE ENCLOSURE AT APPROXIMATE ROOT:
[   -0.1500E+01,  -0.1500E+01 ]
U0:
[    0.3852E+00,   0.3852E+00 ]
U:
[    0.5777E+00,   0.5777E+00 ] [    0.0000E+00,   0.1000E+01 ]
V:
[    0.1926E+00,   0.1926E+00 ]
INEQ_CERT_FEASIBLE:
 F T
NIN_POSS_BINDING:1
Number of bisections: 1
BEST_ESTIMATE:   -0.1500E+01
Total number of boxes processed in loop:            4
Overall CPU time:    0.5000D-01

Simple Example of “Thick” Constants

minimize (x_{1} − [1, 2])^{2} + (x_{2} − [3, 4])^{2}

PROGRAM THICK_PARAMETER_EXAMPLE
USE CODELIST_CREATION
IMPLICIT NONE
INTEGER, PARAMETER:: NN=2
TYPE(CDLVAR), DIMENSION(NN):: X
TYPE(CDLLHS) :: PHI
OUTPUT_FILE_NAME='thick_parameter_example.CDL'
CALL INITIALIZE_CODELIST(X)
PHI = (X(1) - INTERVAL(1,2))**2 &
      + (X(2)-INTERVAL(3,4))**2
CALL FINISH_CODELIST
END PROGRAM THICK_PARAMETER_EXAMPLE

Example with Thick Constants

(continued)

The “solution” is the set of all possible minima with constants in [1, 2] and [3, 4]. Thus, a minimum minimum and maximum minimum are obtained. The solution is

x_{1} ∈ [1, 2], x_{2} ∈ [3, 4], and ϕ = 0. With initial box ([-10,10], [-10,10]), GlobSol gives

Box no.:           1
Box coordinates:
[ 0.1000D+01, 0.1500D+01 ] [ 0.3000D+01, 0.4000D+01 ]
PHI:
[ 0.0000D+00, 0.2000D+01 ]
Box no.:           2
Box coordinates:
[ 0.1500D+01, 0.2000D+01 ] [ 0.3000D+01, 0.4000D+01 ]
PHI:
[ 0.0000D+00, 0.2000D+01 ]
BEST_ESTIMATE:    0.5000D+00
Total number of boxes processed in loop:            4
Overall CPU time:    0.0000D+00

GlobSol July, 2000 Dresden–28

Thick Constants

A Practical Application

• Work by Claudio Rocco

• Maintenance scheduling optimization based on a model by Dekker et al

• A multi-component system has n components. Each component can be scheduled for preventive maintenance every k_{i}T time units, where T is a fixed interval and the k_{i} are integers.

• The total cost of a particular schedule is written as an unconstrained minimization problem, with integer variables k_{i} and the base interval T.

• The expected deterioration cost of each component is not known precisely; the cost of the i-th component depends on a parameter c_{i} known to lie in an interval.

GlobSol July, 2000 Dresden–29

Maintenance Scheduling Application

GlobSol’s Possibilities

1. For an objective of the form f(x, p), where p is a set of parameters subject to interval uncertainty, GlobSol can compute a lower bound on

min_{p∈p} {min_{x} f(x, p)}

and an upper bound on

max_{p∈p} {min_{x} f(x, p)} .

2. This lower bound and upper bound can be made sharper by subdividing the parameter intervals p, solving the subproblems, then taking the union of the solutions. Subdividing can also decrease the overall execution time.

3. Details can be found in a preprint at

http://interval.usl.edu/preprints/TOMS thick.ps

GlobSol July, 2000 Dresden–30

Maintenance Scheduling Application

GlobSol’s Performance

1. GlobSol solves an 8-component model with fixed c_{i} very effectively.

2. If the c_{i} are intervals, GlobSol gives satisfactory results if only one of the c_{i} is assumed to be uncertain, and if the variation is within the range of 10% or so.

3. GlobSol presently takes excessive time when more than one c_{i} is assumed to be uncertain.

4. Certain algorithmic improvements appear promising for increasing the practicality of GlobSol on this problem when more than one parameter is uncertain.

GlobSol References

• For the source, installation instructions, user guide, etc.:

http://www.mscs.mu.edu/~globsol/

• Rigorous Global Search: Continuous Problems, R. B. Kearfott, Kluwer Academic Publishers, 1996. Contains

– Most of the basic ideas underlying GlobSol.

– Structure of the research code that eventually became GlobSol.

• For various preprints related to techniques in GlobSol and applications:

http://interval.louisiana.edu/preprints.html

• For these transparencies:

http://interval.louisiana.edu/preprints/2000 Dresden.ps

(Postscript)

http://interval.louisiana.edu/preprints/2000 Dresden.dvi

(TEX DVI)

GlobSol July, 2000 Dresden–32