Constrained Minimization Using Pattern Search, Problem-Based

Constrained Minimization Problem

For this problem, the objective function to minimize is a simple function of 2-D variables X and Y:

camxy = @(X,Y)(4 - 2.1.*X.^2 + X.^4./3).*X.^2 + X.*Y + (-4 + 4.*Y.^2).*Y.^2;

This function is known as "cam," as described in L.C.W. Dixon and G.P. Szego [1].

Additionally, the problem has nonlinear constraints and bounds.

   x(1)*x(2) + x(1) - x(2) + 1.5 <= 0  (nonlinear constraint)
   10 - x(1)*x(2) <= 0                 (nonlinear constraint)
   0 <= x(1) <= 1                      (bound)
   0 <= x(2) <= 13                     (bound)

Plot the nonlinear constraint region on a surface plot of the objective function. The constraints limit the solution to the small region above both red curves.

x1 = linspace(0,1);
y1 = (-x1 - 1.5)./(x1 - 1);
y2 = 10./x1;
[X,Y] = meshgrid(x1,linspace(0,13));
Z = camxy(X,Y);
surf(X,Y,Z,"LineStyle","none")
hold on
z1 = camxy(x1,y1);
z2 = camxy(x1,y2);
plot3(x1,y1,z1,'r-',x1,y2,z2,'r-')
xlim([0 1])
ylim([0 13])
zlim([0,max(Z,[],"all")])
hold off

Create Optimization Variables, Problem, and Constraints

To set up this problem, create optimization variables x and y. Set the bounds as you create the variables.

x = optimvar("x","LowerBound",0,"UpperBound",1);
y = optimvar("y","LowerBound",0,"UpperBound",13);

Create the objective as an optimization expression.

cam = camxy(x,y);

Create an optimization problem with this objective function.

prob = optimproblem("Objective",cam);

Create the two nonlinear inequality constraints, and include them in the problem.

prob.Constraints.cons1 = x*y + x - y + 1.5 <= 0;
prob.Constraints.cons2 = 10 - x*y <= 0;

Review the problem.

show(prob)

  OptimizationProblem : 

	Solve for:
       x, y

	minimize :
       (((((4 - (2.1 .* x.^2)) + (x.^4 ./ 3)) .* x.^2) + (x .* y)) + (((-4) + (4 .* y.^2)) .* y.^2))


	subject to cons1:
       ((((x .* y) + x) - y) + 1.5) <= 0

	subject to cons2:
       (10 - (x .* y)) <= 0

	variable bounds:
       0 <= x <= 1

       0 <= y <= 13

Set Initial Point and Solve

Set the initial point as a structure with field x equal to 0.5 and y equal to 0.5.

x0.x = 0.5;
x0.y = 0.5;

Solve the problem specifying the patternsearch solver.

[sol,fval] = solve(prob,x0,"Solver","patternsearch")

Solving problem using patternsearch.
Optimization finished: mesh size less than options.MeshTolerance 
and constraint violation is less than options.ConstraintTolerance.

sol = struct with fields:
    x: 0.8122
    y: 12.3122

fval = 
9.1324e+04

patternsearch finds the solution point x = 0.8122, y = 12.3122 with objective function value 9.1324e4.

Add Visualization

To observe the solver's progress, specify options that select two plot functions. The plot function psplotbestf plots the best objective function value at every iteration, and the plot function psplotmaxconstr plots the maximum constraint violation at every iteration. Set these two plot functions in a cell array. Also, display information about the solver's progress in the Command Window by setting the Display option to 'iter'.

options = optimoptions(@patternsearch,...
    "PlotFcn",{@psplotbestf,@psplotmaxconstr},...
    "Display","iter");

Run the solver, including the options argument.

[sol,fval] = solve(prob,x0,"Solver","patternsearch","Options",options)

Solving problem using patternsearch.

                                      Max
  Iter   Func-count       f(x)      Constraint   MeshSize      Method
    0         1     0.373958         9.75       0.9086    
    1        18       113581    1.617e-10        0.001   Increase penalty
    2       147        92267            0        1e-05   Increase penalty
    3       373      91333.2            0        1e-07   Increase penalty
    4       638        91324            0        1e-09   Increase penalty
Optimization finished: mesh size less than options.MeshTolerance 
and constraint violation is less than options.ConstraintTolerance.

sol = struct with fields:
    x: 0.8122
    y: 12.3122

fval = 
9.1324e+04

Nonlinear constraints cause patternsearch to solve many subproblems at each iteration. As shown in both the plots and the iterative display, the solution process has few iterations. However, the Func-count column in the iterative display shows many function evaluations per iteration. Both the plots and the iterative display show that the initial point is infeasible, and that the objective function is low at the initial point. During the solution process, the objective function value initially increases, then decreases to its final value.

Unsupported Functions

If your objective or nonlinear constraint functions are not Supported Operations for Optimization Variables and Expressions, use fcn2optimexpr to convert them to a form suitable for the problem-based approach. For example, suppose that instead of the constraint $$xy\ge 10$$ you have the constraint $$I_1(x)+I_1(y)\ge 10$$, where $$I_1(x)$$ is the modified Bessel function besseli(1,x). (The Bessel functions are not supported functions.) Create this constraint using fcn2optimexpr as follows. First create an optimization expression for $$I_1(x)+I_1(y)$$.

bfun = fcn2optimexpr(@(t,u)besseli(1,t) + besseli(1,u),x,y);

Next, replace the constraint cons2 with the constraint bfun >= 10.

prob.Constraints.cons2 = bfun >= 10;

Solve the problem. The solution differs because the constraint region is different.

[sol2,fval2] = solve(prob,x0,"Solver","patternsearch","Options",options)

Solving problem using patternsearch.

                                      Max
  Iter   Func-count       f(x)      Constraint   MeshSize      Method
    0         1     0.373958        9.484       0.9307    
    1        18       113581            0        0.001   Increase penalty
    2       183      962.841            0        1e-05   Increase penalty
    3       499      960.942            0        1e-07   Increase penalty
    4       636       960.94            0    8.511e-15   Update multipliers
Optimization finished: mesh size less than options.MeshTolerance 
and constraint violation is less than options.ConstraintTolerance.

sol2 = struct with fields:
    x: 0.4998
    y: 3.9981

fval2 = 
960.9401

References

[1] Dixon, L. C. W., and G .P. Szego (eds.). Towards Global Optimisation 2. North-Holland: Elsevier Science Ltd., Amsterdam, 1978.