nlpdemo.html

How to Use the NonlinearProgramming Package

Algorithms for Constrained Nonlinear Optimization Problems

Introduction

This is a package for solving nonlinear optimization problems. The NonlinearProgramming package contains iterative algorithms for finding local minima and maxima of nonlinear functions subject to nonlinear inequality constraints. The algorithms used are Newton's Method for unconstrained problems, and a primal-dual log-barrier method for constrained problems. Both algorithms are implemented with special adaptations for nonconvex problems.

You may run either of these two iterative algorithms ( UnconstrainedNewton(...) or PrimalDualLogBarrier(...) ) on your optimization problem directly, or you may call a top-level procedure Optimize(...) that checks your optimization problem for special structure before passing it on to one of the iterative algorithms. In particular, Optimize(...) checks (1) whether your problem is a linear program, in which case Maple's built-in simplex package is used; (2) whether your problem is an unconstrained convex quadratic program, in which case the problem is solved explicity using first-order optimality conditions; and (3) whether your problem is a constrained convex quadratic program, in which case the adaptations intended for nonconvex programs are turned off to save computation time.

The NonlinearProgramming package was developed by Jason Schattman of Waterloo Maple, Inc. The author will cheerfully answer any questions about the package, or about optimization modeling in general. Please direct questions to jschattman@maplesoft.com.

Installation Instructions

1. Download the Maple worksheet NLPcode.mws from the NonlinearProgramming Package homepage.

2. Go to the last section of the worksheet NLPcode.mws and follow the instructions there.

3. Execute the worksheet NLPcode.mws. This saves the package as a Maple library on your machine.

4. From now on, in any Maple worksheet where you want to use the NonlinearProgramming package, load the packages LinearAlgebra and linalg by typing:

> with( LinearAlgebra ): with( linalg ):

5. Load the NonlinearProgramming package by typing

> with( NonlinearProgramming );

Examples

For all examples, we'll need Maple's built-in LinearAlgebra and linalg packages.

> restart;

> with( LinearAlgebra ): with( linalg ):

Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value

Warning, the protected names norm and trace have been redefined and unprotected

We now load the NonlinearProgramming package. You will need to set your libname path to where you have saved the NonlinearProgramming package.

> libname := "C:/mylib/nlp",libname;

> with( NonlinearProgramming );

We would like to see the individual iterates, so we set the user feedback level on all procedures to 2. To see the final answer only, use level 1.

> infolevel['UnconstrainedNewton'] := 2: infolevel['Optimize'] := 2: infolevel['PrimalDualLogBarrier'] := 2:

Example 1: Unconstrained nonlinear program

minimize

The NonlinearProgramming package requires the decision variables to be input in terms of , ..., . We define the function to be minimized as

> f := x[1]^2 + x[2]^2 + x[3]^2;

We have 3 decision variables.

> numDecVars := 3;

Let's start our search at the point (4, 4, 4).

> x_start:=<4, 4, 4>;

Because this is an unconstrained problem, we'll use UnconstrainedNewton(...).

> UnconstrainedNewton( f, numDecVars, x_start, 'convex', 'float[8]' );

UnconstrainedNewton:   "Starting point"

UnconstrainedNewton:

UnconstrainedNewton:

UnconstrainedNewton:

UnconstrainedNewton:

UnconstrainedNewton:

UnconstrainedNewton:   Global optimum found at

The optimal solution is (0, 0, 0), and Newton's Method found it, albeit with some roundoff error. We note that the objective function is quadratic. The top-level procedure Optimize(...) discovers this and solves the problem explicity without having to call UnconstrainedNewton(...) .

> Optimize( f, min, numDecVars, {}, 'free', x_start );

OptimizeExplicitFONOC:   Convex quadratic problem: solving explicitly with first-order conditions

Example 2: Constrained nonlinear program

minimize

subject to

We can think of this optimization problem as finding the point inside the sphere of radius 9 centered at (10, 0, 0) that is closest to the origin. Because this is a constrained problem, we use PrimalDualLogBarrier(...) . As with the objective function, constraints must be entered in terms of , ..., . The constraints are entered as a set, in this case with only one element.

> constraintSet := { (x[1]-10)^2 + x[2]^2 + x[3]^2 <= 81 };

$constraintSet := {(x[1]-10)^2+x[2]^2+x[3]^2 <= 81}$

Both the objective function (a parabola) and the feasible region (a sphere) are convex, so we bypass the adaptations for nonconvex problems by setting the convexity parameter to 'convex'.

> PrimalDualLogBarrier( f, numDecVars, constraintSet, 'free', x_start, 'convex' );

PrimalDualLogBarrier:   Minimize

PrimalDualLogBarrier:   subject to   [0 <= -(x[1]-10)^2-x[2]^2-x[3]^2+81]

PrimalDualLogBarrier:   Starting point

PrimalDualLogBarrier:

InitialBarrierParameter:   Initial barrier parameter   2600.

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:   Global optimum found at

The algorithm found the global optimum at (1, 0, 0).

Example 3: Linear Program

Maximize x + y

subject to

x + 2y >= 2

x - y <= 2

4x - y >= 0

x <= 4

y <= 6.5

> linearConstraints := {x[1]+2*x[2]>=2, x[1]-x[2]<=2, x[2]<=4*x[1], x[1]<=4, x[2]<=6.5};

For two-dimensional linear programs, we can graph the feasible region.

> with(plots): inequal(linearConstraints, x[1]=-10..10, x[2]=-10..10, optionsexcluded=(color=green,thickness=4), optionsfeasible=(color=blue,thickness=4));

Warning, the name changecoords has been redefined

> f_linear := x[1] + x[2];

Let us compare the performance of Optimize(...) and PrimalDualLogBarrier(...) on this linear program.

> Optimize( f_linear, max, 2, linearConstraints, free, <0,0> );

OptimizeSimplex:   Linear problem: solving with simplex method

Optimize(...) recognized the linear structure of the problem and called Maple's built-in simplex method. PrimalDualLogBarrier(...) , however, does no structure checking. We must also remember to enter -1*f_linear because we wish to maximize, and PrimalDualLogBarrier(...) takes minimization problems as input.

> PrimalDualLogBarrier( -1 * f_linear, 2, linearConstraints, 'free', <0, 0>, 'convex' );

PrimalDualLogBarrier:   Minimize

PrimalDualLogBarrier:   subject to   [0 <= -x[1]+4, 0 <= -x[2]+6.5, 0 <= -x[2]+4*x[1], 0 <= -x[1]+x[2]+2, 0 <= -2+x[1]+2*x[2]]

PrimalDualLogBarrier:   Starting point

PrimalDualLogBarrier:

InitialBarrierParameter:   Initial barrier parameter   325.5764119

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:   Global optimum found at

After some rounding, we see that both methods arrive at the same answer.

Example 4: Nonconvex, higher dimensional problem

minimize

subject to

, , ..., all nonnegative

Until now, our examples have been quite tame. Let's turn up the screws a bit and try a nonconvex, 10-dimensional problem with 11 constraints. We'll also throw in an infeasible starting point.

> f_tenD:=add(x[i]^4, i=1..10) + sin( product(x[i],i=1..10 ) );

> numDecVars := 10;

> tenDSphere:={(x[1]-10)^2+ add(x[i]^2, i=2..10)<=81};

$tenDSphere := {(x[1]-10)^2+x[2]^2+x[3]^2+x[4]^2+x[5...$

> Optimize( f_tenD, min, numDecVars, tenDSphere, nonnegative, <5,1,5,1,5,0,5,2,1,1> );

OptimizeInteriorPointMethod:   Constrained   nonconvex   problem: solving with   nonconvex   primal-dual log-barrier algorithm

PrimalDualLogBarrier:   Minimize

PrimalDualLogBarrier:   subject to   [0 <= x[1], 0 <= x[2], 0 <= x[3], 0 <= -(x[1]-10)^2-x[2]^2-x[3]^2-x[4]^2-x[5]^2-x[6]^2-x[7]^2-x[8]^2-x[9]^2-x[10]^2+81, 0 <= x[4], 0 <= x[5], 0 <= x[6], 0 <= x[7], 0 <= x[8], 0 <= x[9], 0 <= x[10]]

PrimalDualLogBarrier:   Starting point

PrimalDualLogBarrier:

InitialBarrierParameter:   Initial barrier parameter   526.0173142

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:

PrimalDualLogBarrier:   Local optimum found at

Some reality checking confirms that (1, 0, 0, 0, 0, 0, 0, 0, 0, 0) is the optimal solution.

Documentation

Entering unconstrained nonlinear problems

If you want to minimize a function with no constraints and don't want the algorithm to check for special structure, call:

UnconstrainedNewton (<function f to be minimized>, <number of decision variables>, <starting point>, <convexity of f ('convex' or 'nonconvex')>, <data type of the range of f ('float[8]' or 'complex[8]'));

The function f must be in terms of , ..., (as opposed to x, y, z etc.). If you want to maximize f , enter -1* f .

The starting point must be a vector entered in form < v1, v2, ... >.

If you know f is convex, enter ' convex'. This turns off extra computation applied to nonconvex problems. Otherwise, enter ' nonconvex' .

If you know f will never take the square root of a negative number, use ' float[8]' as the data type. Otherwise, enter ' complex[8]' , and the algorithm will take the real part of all function evaluations.

If you want to optimize a function with no constraints and want the algorithm to check for special structure (e.g. convex quadratic), call:

Optimize( <function f to be minimized OR maximized>, <objective ('max' or 'min')>, <number of decision variables>, {}, free, <starting point>);

Note that Optimize(...) assumes f is always real-valued. If f can evaluate to a complex number, use UnconstrainedNewton(...) with the datatype argument

' complex[8]' .

Entering constrained nonlinear problems

If you want to minimize a function subject to inequality constraints and don't want the algorithm to check for special structure, call:

PrimalDualLogBarrier (<function f to be minimized>, <number of decision variables>, <set of inequality constraints>, <variable bounding ('nonnegative' or 'free')>, <starting point>, <convexity of the problem ('convex' or 'nonconvex')>);

The function f and the constraint functions must be in terms of , ..., (as opposed to x, y, z etc.).

If you want to maximize f , enter -1* f .

The inequality constraints must be in the form of a set {<constraint 1>, ..., <constraint m>}.

If your decision variables must be nonnegative, use the parameter 'nonnegative' . Otherwise, use ' free' .

The starting point must be a vector < v1, v2, ... >.

If you know that both f and the feasible region defined by the constraints are convex, use ' convex' as the convexity parameter . This turns off extra computation applied to nonconvex problems. Otherwise, enter ' nonconvex' .

If you want to optimize a function subject to inequality constraints and want the algorithm to check for special structure (e.g. linear or quadratic), call:

Optimize( <function f to be minimized OR maximized>, <objective ('max' or 'min')>, <number of decision variables>, <set of inequality constraints>, <variable bounding ('nonnegative' or 'free')> , <starting point>);

Output options

If all you want to see is the final answer, simply call the procedures as described above. However, if you want to monitor the progress of the algorithm in any way, set the infolevel for the algorithm you're using to the appropriate level. For example,

> infolevel['PrimalDualLogBarrier'] := 2: infolevel['Optimize']:=1:

Infolevel 1 Echo the objective function and constraints and report which algorithm is being used (e.g. simplex method, unconstrained newton)

Infolevel 2 Show the values of the decision variables in every iteration, plus Infolevel 1 information.
Infolevel 3 Show all diagnostic information: the symbolic gradient vector and hessian matrix of the objective function and the constraint functions, dual variable iterates, progress on primal and dual feasibility, current barrier parameter and stepsizes, reports on augmentations of the hessian's diagonal due to nonconvexity, plus Infolevel 2 information.

Hints and Troubleshooting

> restart;

Problem

> with(NonlinearProgramming);

Error, (in pacman:-pexports) NonlinearProgramming is not a package

Solution : You didn't set libname to point to the directory where you've saved the NonlinearProgramming package. Execute first

> #libname:="C:/mylib/nlp",libname:

Problem

> Optimize(x[1]^2+x[2]^2, min, 2, {}, 'free', <0,0>);

Error, (in Matrix) dimension parameters are required for this form of initializer

Solution : You didn't load the LinearAlgebra and/or linalg packages. The NonlinearProgramming package uses these packages. Execute first

> #with(LinearAlgebra): with(linalg):

Problem

> UnconstrainedNewton(-x[1]^2, 1, <5>, 'convex', 'float[8]');

Error, (in LinearAlgebra:-LA_Main:-LinearSolve) Matrix is not positive-definite

Solution : The convexity parameter was set to 'convex', but the objective function given, -x^2, is not convex. This causes an error because the algorithm solves linear systems using a Cholesky decomposition of the objective function's Hessian matrix. This only works when this matrix is positive-definite (PD). Hessians of nonconvex objective functions are, in general, not PD. To fix this, set the convexity parameter to 'nonconvex' instead. Then the algorithm augments the Hessian's diagonal until it is PD. Instead of the above, execute

> #UnconstrainedNewton(-x[1]^2, 1, <5>, 'nonconvex', 'float[8]');

Problem: The algorithm keeps repeating the same iterate over and over and never terminates.
Solutions :
(1) Try a different starting point, in particular, one that's close to the iterate that is being repeated. In general, nonlinear programming algorithms are quite sensitive to the starting point.
(2) If you didn't specify nonnegative variable bounding, and one or more of your decision variables will never be negative, try adding some nonnegativity constraints, either by using the 'nonnegative' parameter, or by adding individually to the constraints, for example, {x[17] >= 0}. Adding nonnegativity constraints can sometimes provide stability to the primal-dual log-barrier algorithm on difficult problems.

Problem : The iterates increase in magnitude without bound.
Solution : Your problem may be unbounded. For example, suppose you want to find the local maximum of the cubic (x[1] - 2)(x[1]-4)(x[1]-6). If you enter the starting point <7> and specify a maximization objective, the algorithm will increase the value of x[1] without bound because the cubic is increasing on the interval [7, +infinity). To find the local maximum of the cubic, you would need a starting point between the first two roots, i.e., in the interval [2,4], or add constraints on x[1], e.g. {x[1] >= -5, x[1] <= 10}.

References

The general idea behind the primal-dual log-barrier algorithm implemented in this package is from

Vanderbei, R. and Shanno, D., An Interior Point Algorithm for Nonconvex Nonlinear Programming , Technical Report SOR 97-21, Statistics and Operations Research, Princeton University, 1997 .

For a general mathematical treatment of nonlinear programming theory, see

Bertsekas, D., Nonlinear Programming, Athena Scientific, 1995.