Main Content

Constrained Minimization Using Pattern Search, Problem-Based

This example shows how to minimize an objective function, subject to nonlinear inequality constraints and bounds, using pattern search in the problem-based approach. For a solver-based version of this problem, see Constrained Minimization Using Pattern Search, Solver-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

Figure contains an axes object. The axes object contains 3 objects of type surface, line.

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 terminated: 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

View the solution point.

disp(sol.x)
    0.8122

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       148        92267            0        1e-05   Increase penalty
    3       374      91333.2            0        1e-07   Increase penalty
    4       639        91324            0        1e-09   Increase penalty
Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.

Figure Pattern Search contains 2 axes objects. Axes object 1 with title Best Function Value: 91324 contains an object of type line. Axes object 2 with title Max Constraint Violation: 0 contains an object of type line.

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 xy10 you have the constraint I1(x)+I1(y)10, where I1(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 I1(x)+I1(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    1.659e-08        0.001   Increase penalty
    2        78      962.841            0        1e-05   Increase penalty
    3       394      960.942            0        1e-07   Increase penalty
    4       531       960.94            0    8.511e-15   Update multipliers
Optimization terminated: mesh size less than options.MeshTolerance
 and constraint violation is less than options.ConstraintTolerance.

Figure Pattern Search contains 2 axes objects. Axes object 1 with title Best Function Value: 960.94 contains an object of type line. Axes object 2 with title Max Constraint Violation: 0 contains an object of type line.

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.

See Also

|

Related Topics