Nonlinear Constrained Optimization

M. Asghar Bhatti
University of Iowa, Iowa City, IA

Optimization Problems

Find x to 	Minimize f(x)

Subject to	 (Subscript[h, i])(x)=0,i = 1,2,…,p 	and	 (Subscript[g, i])(x)≤0,i = 1,2,…,m.

Building design example

To save energy costs for heating and cooling an architect is considering designing a partially buried rectangular building. The total floor space needed is 20,000 (m^2). Lot size limits the building plan dimension to 50 m. It has already been decided that the ratio between the plan dimensions must be equal to the golden ratio (1.618) and that each story must be 3.5 m high. The heating and cooling costs are estimated at $100 per (m^2)of the exposed surface area of the building. The owner has specified that the annual energy costs should not exceed $225, 000. Determining building dimensions to minimize cost of excavation.

Optimization variables

n = Number of stories
d = Depth of building below ground
h = Height of building above ground
ℓ = Length of building in plan
w = Width of building in plan

Objective function

Minimize excavation cost

Constraints

Number of stories related to the building height

Plan dimensions

Floor space requirement

Lot size

Energy cost

Physical limitations

The complete optimization problem can be stated as follows.

Find (n, ℓ, w, h, d) in order to

Minimize  d ℓ w

Subject to (({{((d + h)/n) = 3.5}, {ℓ = 1.618  w}, {n ℓ w ≥ 20, 000}, {ℓ ≤ 50}, {w ≤ 50}, {\(100 \((2\ h\ ℓ\  +                         2\ h\ w\  + \ ℓ\ w)\) ≤ 225,
                000\)}, {n ≥ 1}, {ℓ, w, h, d ≥ 0}}))

Built-in Support

Unconstrained problems — FindMinimum & FindMaximum

Information["FindMinimum", LongForm -> False]

f = ((((x^2) + (y^2) - 3 x))^2) - ((4 (x^2))/(2 + x)); 
vars = {x, y}; x0 = {1, 1};

FindMinimum[f, {x, 1}, {y, 1}]

Linear programming problems — ConstrainedMin, ConstrainedMax, LinearProgramming

Information["LinearProgramming", LongForm -> False]

Information["ConstrainedMin", LongForm -> False]

f = 2 (Subscript[x, 1]) + (Subscript[x, 2]) + 3 (Subscript[x, 4]); 
g = {9 (Subscript[x, 1]) + 14 (Subscript[x, 2]) - 6 (Subscript[x, 3]) - 6 (Subscript[x, 4]) ≤ 2,
        (Subscript[x, 1]) + (Subscript[x, 2]) - (Subscript[x, 3]) - (Subscript[x, 4]) ≥ -1, -20 (Subscript[x, 1]) -
            5 (Subscript[x, 2]) + 5 (Subscript[x, 3]) + 13 (Subscript[x, 4]) == 11}; 
vars = Table[(Subscript[x, i]), {i, 1, 4}];

ConstrainedMin[f, g, vars]

Nonlinear Constrained Problems

Large number of methods

Direct Methods

Sequential linear programming

Sequential quadratic programming

Feasible directions

Indirect Methods — Penalty methods

Interior and exterior penalty functions

Augmented Lagrangian Penalty Function method

Penalty Methods

Exterior Penalty Function

ϕ(x,μ)=f(x)+μ [(∑_(i = 1)^m)(((Max[0,(Subscript[g, i])(x)]))^2)+(∑_(i = 1)^p)((((Subscript[h, i])(x)))^(2 ))]

where μ is a large positive number, known as the penalty number. If an inequality constraint is satisfied, then (Subscript[g, i]) < 0 and Max[0, (Subscript[g, i])] will return 0 and therefore that constraint will not contribute anything to the function ϕ. If a constraint is violated, i.e. (Subscript[g, i]) > 0 or (Subscript[h, i]) ≠ 0, a large term will get added to ϕ. Minimizing ϕ should force the violation to go away. For the method to work the penalty number μ must be very large. In fact, theoretically, the minimum of ϕ corresponds to the solution of the original problem only as μ -> ∞.

Example

Minimize	f = (1/2) (Subscript[x, 1]^2) + (Subscript[x, 1]) (Subscript[x, 2]^2)

Subject to	(({{(Subscript[x, 1]) (Subscript[x, 2]^2) = 10}, {-2 (Subscript[x, 1]^2) + (Subscript[x, 2])/3 ≤ 0}}))

Optimum solution: \(TraditionalForm\`\(f\^*\) = 10.2995359848948342`, \((x\_1 = 0.774,
        x\_2 = 3.594)\)\).

The unconstrained function ϕ for the problem is as follows.

ϕ =
      ((Subscript[x, 1]^2)/2) +
        (Subscript[x, 1]) (Subscript[x, 2]^2) + μ (((((Subscript[x, 1]) (Subscript[x, 2]^2) - 10))^2) +
              Max[0, -2 (Subscript[x, 1]^2) + ((Subscript[x, 2])/3)]^2);

Contour plots of this function for two different μ values. The plots demonstrate that the minimum of the unconstrained problem is close to the optimum solution. The figure also illustrates numerical difficulty anticipated in minimizing ϕ for large μ values. As the value of μ is increased the contours become closer together with a long narrow shape.

Interior Penalty Function

In the interior penalty function method, the penalty term is defined to keep the solution from leaving the feasible region. It is not possible to directly handle equality constraints in this method. Thus the original problem for this method is defined as follows.

Find x to 	Minimize f(x)

Subject to	 (Subscript[g, i])(x)≤0,i = 1,2,…,m.

Two commonly used interior penalty functions are defined as follows.

ϕ(x,μ)=f(x)-μ [(∑_(i = 1)^m)Log[-(Subscript[g, i])(x)]] — Logarithmic barrier function

ϕ(x,μ)=f(x)+μ [(∑_(i = 1)^m)((-1)/((Subscript[g, i])(x)))] — Inverse barrier function

where μ is a positive number, known as the penalty number. It is assumed that the inequality constraints are satisfied at all time and therefore (Subscript[g, i]) < 0 and a positive term is added to the objective function. As we move closer to a constraint boundary, (Subscript[g, i]) -> 0 causing a large term to be added to the objective function. Thus the method keeps the solution away from the constraint boundaries and hence is also known as the barrier function method. Minimizing ϕ should force the violation to go away. If the solution is at a constraint boundary, for the method to work, the penalty number μ must be very small. In fact, theoretically, the minimum of ϕ corresponds to the solution of the original problem only as μ -> 0.

Example

Minimize	f =  x y  +  (y^2)

Subject to	(({{ (1/4) (y^2)  -  1 ≤ 0}, {x/5  +  y/2  -  1 ≤ 0}}))

Optimum solution: \(TraditionalForm\`\(f\^*\) = \(-16\), \((x = 10, \ y = \(-2\))\)\).

ϕ =
      x y + (y^2) + μ (-(1/(((y^2)/4) - 1)) -
              (1/((x/5) + (y/2) - 1)));

The necessary condition for minimum of this unconstrained function is that the gradient must be zero. This gives the following two equations.

eqns = Thread[Grad[ϕ, {x, y}] == 0]; MatrixForm[eqns]

Several solutions of this system of equations, for different values of penalty parameter μ, are as follows.

({{μ = 1, \({x -> 7.39260175769682348`,
                y -> \(-1.65276800315081278`\)}, {x ->
                  9.29953147342448005`,
                y -> \(-2.30849435651477907`\)}, 
	{x ->
                  10.9855165408315414`,
                y -> \(-1.71027630051842898`\)}, {x ->
                  12.1669615584025336`,
                y -> \(-2.27360301812060416`\)}\)}, {μ = 0.1, \({x -> 9.23956364740620372`,
                y -> \(-1.9009688334254065`\)}, {x ->
                  9.75210163977435584`,
                y -> \(-2.09619761314683028`\)}, 
	{x ->
                  10.2763753098701715`,
                y -> \(-1.90565937222678876`\)}, {x ->
                  10.7192859604458123`,
                y -> \(-2.0921694519929499`\)}\)}, {μ = 0.0001, \({x -> 9.99188115110397312`,
                y -> \(-1.99698356402931658`\)}, {x ->
                  9.99303706130981517`,
                y -> \(-1.99698372265408234`\)}, 
	{x ->
                  10.007165725564957`,
                y -> \(-2.00301392779849507`\)}, {x ->
                  10.0079036645412444`,
                y -> \(-2.00301382679499528`\)}\)}})

Clearly as the value of μ is approaching 0, the solutions of the unconstrained problem are all approaching the solution of the original constrained problem. It is also interesting to note that setting μ = 0 gives x = y = 0 which is obviously not the correct solution. Thus the correct solution can only be obtained through the limiting process as μ -> 0.

ϕ =
      x y + (y^2) - μ (Log[-(((y^2)/4) - 1)] +
              Log[-((x/5) + (y/2) - 1)]);

The necessary conditions for minimum of this unconstrained function are that the gradient must be zero. This gives the following two equations.

eqns = Thread[Grad[ϕ, {x, y}] == 0]; MatrixForm[eqns]

Solutions of this system of equations for different values of penalty parameter μ are as follows.

({{μ = 1, {x -> 0.000775589040265045692`,
                y -> -0.183239978857939114`}, {x ->
                  9.23548061115235441`,
                y -> -1.90424882657261136`}}, {μ = 0.1, {x -> -0.00889613832064528864`,
                y -> -0.0197694109132707005`},  {x ->
                  9.92693416567897024`,
                y -> -1.99086543107678828`}}, {μ = 0.0001, \({x -> \(-9.99888001687922844`*^-6\),
                y -> \(-0.0000199997600097596297`\)}, {x ->
                  0.000300083802230132246`,
                y -> 1.99989996747948346`}, 
	{x ->
                  0.83311597601912517`,
                y -> 1.66677760799490108`}, {x ->
                  9.99992727239199652`,
                y -> \(-1.99999090904770824`\)}\)}})

Clearly as the value of μ is approaching 0, the last solution of the unconstrained problem is approaching the solution of the original constrained problem. By checking the second order conditions, it is easy to see that the other solutions are simply inflection points for function ϕ.

Augmented Lagrangian Penalty Function — Equality constraints

The unconstrained augmented Lagrangian penalty function is defined as follows.

ϕ=f(x)+(∑_(i = 1)^p)(Subscript[v, i]) (((Subscript[h, i])(x))^( ))+μ [(∑_(i = 1)^p)((((Subscript[h, i])(x)))^(2 ))]

It is demonstrated in the following example that if the optimum Lagrange multipliers are known, then the solution of this unconstrained problem corresponds to the solution of the original problem regardless of the value of the penalty parameter.

Example

Consider the following optimization problem

Minimize	f = (1/2) (x^2) + x (y^2)

Subject to	(x +  (y^2) = 10)

Using KT conditions it is easy to compute the optimum solution as follows.

KTSolution[
      ((x^2)/2) + x (y^2), {x + (y^2) - 10 == 0}, {x, y}];

Optimum: ({{x = 10, y = 0, Lagrange multiplier, v = -10}})

In the augmented Lagrangian approach, the unconstrained function is defined by adding the exterior penalty term to the Lagrangian of the original problem. For this example, thus, we have the following unconstrained function.

ϕ =
      ((x^2)/2) + x (y^2) +
        v (x + (y^2) - 10) + μ (((x + (y^2) - 10))^2);

The necessary conditions for minimum of this function give the following equations.

eqns = Thread[Grad[ϕ, {x, y}] == 0]; MatrixForm[eqns]

If v is set to the optimum value of the Lagrange multiplier, we get the following equations.

eqns2 = eqns /. {v -> -10};
    MatrixForm[eqns2]

The second equation can be written as follows.

y (-20  + 2 x  - 40  μ + 4 x  μ +
              4 (y^2) μ) == 0

Thus y = 0 satisfies this equation for any value of μ. Substituting y = 0 in the first equation we get

-10 + x - 20 μ + 2 x μ == 0	or	(-10 + x) (1 + 2 μ) == 0

Thus x = 10 satisfies this equation. Thus the Lagrangian penalty function has the property that the optimum solution of the original problem is recovered, if we know the optimum values of the Lagrange multipliers. Therefore in this sense it is an exact penalty function.

Obviously when we are solving a problem we don't know the optimum Lagrange multipliers. (If they were known we don't need to spend time in developing new algorithms. We can simply use them with the KT conditions to get a solution). However the presence of Lagrange multipliers makes the choice of penalty parameter less critical. In a computational procedure based on using the augmented Lagrangian penalty function we start with arbitrary values of Lagrange multipliers and develop a procedure that moves these multipliers closer to their optimum values. Thus near the optimum, the function ϕ is not as sensitive to values of μ and procedure converges to the true optimum solution.

Augmented Lagrangian Penalty Function — General case

The unconstrained ϕ function can be written as follows.

ϕ(x,u,v,μ)=f(x)+(∑_(i = 1)^p)(Subscript[v, i]) (((Subscript[h, i])(x))^( ))+μ (∑_(i = 1)^p)((((Subscript[h, i])(x)))^(2 ))+μ (∑_(i = 1)^m) (Subscript[r, i]^2)-(∑_(i = 1)^m) ((Subscript[u, i]^2)/(4  μ))

where

(Subscript[r, i])≡Max[(Subscript[g, i])(x)+((Subscript[u, i])/(2  μ)),0],i = 1,…,m

Note that because of the presence of Max, the function (Subscript[r, i]) is not differentiable. However (Subscript[r, i]^2) is continuous and has well defined first derivatives as follows.

({{∂((((Subscript[r, i]^2)))/∂)(Subscript[x, j]) =
                2  (Subscript[r, i]) ∂((Subscript[g, i])/∂)(Subscript[x, j]), j =
                1, …, n}})

The following example numerically illustrates this fact.

Example

For the following function, compute the first derivative as indicated above and graphically compare it with its finite difference approximation.

f(x) = (((Max[0,  Sin[x]]))^2)

The first derivative of the function is as follows.

(df/dx) = 2  Max[0,
            Sin[x]]   Cos[x]

Solution algorithm

(i) Set iteration counter k  =  0. Set multipliers (Subscript[u, i]) = 0, i = 1, …, m and (Subscript[v, i]) = 0, i = 1, …, p. Set multiplier update counter ℓ = 0. Choose a penalty parameter μ and a factor c > 1 to increase penalty parameter value during iterations. Typical values are μ = 10 and c = 2.

(ii) Set up the unconstrained minimization problem.

ϕ(x,u,v,μ)=f(x)+(∑_(i = 1)^p)(Subscript[v, i]) (((Subscript[h, i])(x))^( ))+μ (∑_(i = 1)^p)((((Subscript[h, i])(x)))^(2 ))+μ (∑_(i = 1)^m) (Subscript[r, i]^2)-(∑_(i = 1)^m) ((Subscript[u, i]^2)/(4  μ))

Use a suitable unconstrained minimization method to find the minimum point (x^((k)))of this problem.

(iii) Check for convergence.

A simple convergence criteria is to stop when all constraints are satisfied and the objective function is not changing much between successive iterations. Thus stop if the following two conditions are met otherwise continue to step (iv).

Abs[(f((x^((k))))-f((x^((k - 1))))/f((x^((k))))]<ϵ

V((x^((k))))<ϵ

where V((x^((k))))≡Max[(Subscript[g, i])((x^((k)))),i = 1,…,m,Abs[(Subscript[h, i])((x^((k)))),i = 1,…,p]] is the maximum constraint violation.

(iv) Update the multipliers and the penalty parameter

If V((x^((k))))≤(1/4) V((x^((ℓ)))) then update multipliers

({{(Subscript[u, i]^((ℓ + 1)))=(Subscript[u, i]^((ℓ)))+2μ (Subscript[r, i])((x^((k)))), i =
                1, …, m}})

({{(Subscript[v, i]^((ℓ1)))=(Subscript[v, i]^((ℓ)))+2μ (Subscript[h, i])((x^((k)))), i =
                1, …, p}})

Set ℓ = ℓ +1

Else update the penalty parameter

(μ^((k + 1))) = c (μ^((k)))

(v) Update the iteration counter k = k + 1 and go back to step (ii).

ALPF Function

Needs["OptimizationToolbox`ConstrainedNLP`"];
Needs["OptimizationToolbox`GraphicalSolution`"];

Information["ALPF", LongForm -> False]

OptionsUsage[ALPF]

Example

Minimize	f = (1/2) (Subscript[x, 1]^2) + (Subscript[x, 1]) (Subscript[x, 2]^2)

Subject to	(({{(Subscript[x, 1]) (Subscript[x, 2]^2) = 10}, {-2 (Subscript[x, 1]^2) + (Subscript[x, 2])/3 ≤ 0}, {(Subscript[x, 2]) ≥ -4}}))

vars = {(Subscript[x, 1]), (Subscript[x, 2])}; (pt = {3.0, 1.0};)
f = ((Subscript[x, 1]^2)/2) + (Subscript[x, 1]) (Subscript[x, 2]^2);
cons = {(Subscript[x, 1]) (Subscript[x, 2]^2) == 10, -2 (Subscript[x, 1]^2) + ((Subscript[x, 2])/3) ≤
            0, (Subscript[x, 2]) ≥ -4};

The solution is obtained by using the ALPF method. All intermediate calculations are shown for the first two iterations.

ALPF[f, cons, vars, pt, MaxIterations -> 2, PrintLevel -> 3,
      PenaltyParameters -> {10, 2}];

The procedure is allowed to run until an optimum is found.

{sol, hist} =
        ALPF[f, cons, vars, pt, PrintLevel -> 0,
          PenaltyParameters -> {10, 2}];

TableForm[hist]

The following graph confirms the solution and shows the path taken by the ALPF method.

xhist = Transpose[
          Drop[hist,
            1]][[1]]

GraphicalSolution[f, {(Subscript[x, 1]), 0, 5}, {(Subscript[x, 2]), 0, 5},
        ObjectiveContours -> {5, 10, 20, 50, 70},
        Constraints -> cons, ShadingOffset -> 0.07,
        PlotPoints -> 35,
        Epilog -> {RGBColor[1, 0, 0], Disk[Last[xhist], 0.05],
            Text["Minimum", Last[xhist] + { .1, 0}, {-1, 0}],
            Text["Start", First[xhist] + {0, - .1}],
            Line[xhist]}];

Example

vars = {(Subscript[x, 1]), (Subscript[x, 2])};
f = (E^(-(Subscript[x, 1]) - (Subscript[x, 2]))) + (Subscript[x, 1]) (Subscript[x, 2]) +
          (Subscript[x, 2]^2);
g = {((E^(-(Subscript[x, 1])))/4) + ((Subscript[x, 2]^2)/4) - 1 ≤ 0,
          ((Subscript[x, 1])/5) + ((Subscript[x, 2])/2) - 1 ≤ 0};
pt = {1.0, 1.0};

{sol, hist} =
        ALPF[f, g, vars, pt, PrintLevel -> 0,
          PenaltyParameters -> {10, 2}];

TableForm[hist]

xhist = Transpose[
            Drop[hist,
              1]][[1]];

GraphicalSolution[f, {(Subscript[x, 1]), -3, 12}, {(Subscript[x, 2]), -3, 5},
        ObjectiveContours -> {-20, -10, 2, 15, 30},
        Constraints -> g, ShadingOffset -> 0.05,
        PlotPoints -> 35,
        Epilog -> {RGBColor[1, 0, 0], Disk[Last[xhist], 0.1],
            Text["Minimum", Last[xhist] + {0, - .2}],
            Text["Start", First[xhist] + {0, 0.1}],
            Line[xhist]}];

Building design example

Clear[d, h, w]; 
vars = {d, h, w}; 
f = 1.618 d (w^2); 
cons = {((100 (5.236 h w + 1.618 (w^2)))/225000) ≤
            1, ((1.618 w)/50) ≤
            1, ((0.462286 (d + h) (w^2))/20000) ≥ 1,
          d ≥ 0, h ≥ 0,
          w ≥ 0};
pt = {40, 10, 20};

{sol, hist} =
      ALPF[f, cons, vars, pt, MaxIterations -> 20, PrintLevel -> 0,
        MaxUnconstrainedIterations -> 100,
        PenaltyParameters -> {1000, 2}];

TableForm[hist]

TELOS/Springer-Verlag book

Practical Optimization Methods
With Mathematica Applications

M. Asghar Bhatti, University of Iowa, Iowa City, IA

to be published December 1999

Chapter 1—Problem Formulation

.1 Optimization Problem Formulation

.2 Standard Form of an Optimization Problem

.3 Solution of Optimization Problems

.4 Time value of money

Problems

Chapter 2—Graphical Optimization

2.1 Procedure for Graphical Optimization

2.2 GraphicalSolution Function

2.3 Graphical Optimization Examples

Problems

Chapter 3—Mathematical Preliminaries

3.1 Vectors and Matrices

3.2 Approximation Using Taylor Series

3.3 Solution of Nonlinear Equations

3.4 Quadratic Forms

3.5 Convex Functions and Convex Optimization Problems

Problems

Chapter 4—Optimality Conditions

4.1 Optimality Conditions for Unconstrained Problems

4.2 Additive Property of Constraints

4.3 Karush-Kuhn-Tucker (KT) Conditions

4.4 Geometric Interpretation of KT Conditions

4.5 Sensitivity Analysis

4.6 Optimality Conditions for Convex Problems

4.7 Second Order Sufficient Conditions

4.8 Lagrangian Duality

Problems

Chapter 5—Unconstrained Problems

5.1 Descent Direction

5.2 Line Search Techniques—Step Length Calculations

5.3 Unconstrained Minimization Techniques

Problems

Chapter 6—Linear Programming

6.1 Standard LP Problem

6.2 Solution of Linear System of Equations

6.3 Basic Solutions of an LP Problem

6.4 The Simplex Method

6.5 Unusual Situations Arising During the Simplex Solution

6.6 Post-Optimality Analysis

6.7 Revised Simplex Method

6.8 Sensitivity Analysis Using the Revised Simplex Method

Problems

Chapter 7—Interior Point Methods

7.1 Optimality Conditions for Standard LP

7.2 Primal Affine Scaling Method

7.3 Primal-Dual Interior Point Method

Appendix—Null and Range Spaces

Problems

Chapter 8—Quadratic Programming

8.1 KT Conditions for Standard QP

8.2 Primal Affine Scaling Method for QP

8.3 Primal-Dual Method for QP

8.4 Active Set Method

8.5 Active Set Method for Dual QP Problem

Appendix — Derivation of the Descent Direction Formula for PAS method

Problems

Chapter 9—Nonlinear Constrained Problems

9.1 Normalization

9.2 Penalty Methods

9.3 Linearization of a Nonlinear Problem

9.4 Sequential Linear Programming — SLP

9.5 Basic Sequential Quadratic Programming — SQP

9.6 Refined SQP Methods

Problems

Appendix—Introduction to Mathematica

A.1. Basic manipulations in Mathematica

A.2. Lists and Matrices

A.3. Solving equations

A.4. Plotting in Mathematica

A.5. Programming in Mathematica

A.6. Packages in Mathematica

A.7. On-Line Help

Bibliography

CD containing all OptimizationToolbox functions, additional examples