Use Supported Operations When Possible

Generally, create your objective or nonlinear constraint functions by using supported operations on optimization variables and expressions. Doing so has these advantages:

In general, supported operations include all elementary mathematical operations: addition, subtraction, multiplication, division, powers, and elementary functions such as exponential and trigonometric functions and their inverses. Nonsmooth operations such as max, abs, if, and case are not supported. For the complete description, see Supported Operations on Optimization Variables and Expressions.

For example, suppose that your objective function is

where $$r$$ is a parameter that you supply, and the problem is to minimize $$f$$ over $$x$$ and $$y$$ . This objective function is a sum of squares, and takes the minimal value of 0 at the point $$x=r$$, $$y=r^2$$.

The objective function is a polynomial, so you can write it in terms of elementary operations on optimization variables.

r = 2;
x = optimvar('x');
y = optimvar('y');
f = 100*(y - x^2)^2 + (r - x)^2;
prob = optimproblem("Objective",f);
x0.x = -1;
x0.y = 2;
[sol,fval] = solve(prob,x0) 

Solving problem using lsqnonlin.

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

sol = struct with fields:
    x: 2.0000
    y: 4.0000

fval = 1.4644e-20

To solve the same problem by converting the objective function using fcn2optimexpr (not recommended), first write the objective as an anonymous function.

fun = @(x,y)100*(y - x^2)^2 + (r - x)^2;
prob.Objective = fcn2optimexpr(fun,x,y);
[sol2,fval2] = solve(prob,x0)

Solving problem using fminunc.

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

sol2 = struct with fields:
    x: 2.0000
    y: 3.9998

fval2 = 1.7143e-09

Notice that solve uses fminunc this time instead of the more efficient lsqnonlin, and the reported solution for y is slightly different than the correct solution 4. Furthermore, the reported fval is about 1e-9 instead of 1e-20 (the actual solution value is exactly 0). These slight inaccuracies are due to solve not using the more efficient solver.

The remainder of this example shows how to convert a function to an optimization expression using fcn2optimexpr.

Function File

To use a function file in the problem-based approach, you must convert the file to an expression using fcn2optimexpr.

For example, the expfn3.m file contains the following code:

type expfn3.m

function [f,g,mineval] = expfn3(u,v)
mineval = min(eig(u));
f = v'*u*v;
f = -exp(-f);
t = u*v;
g = t'*t + sum(t) - 3;

This function is not entirely composed of supported operations because of min(eig(u)). Therefore, to use expfn3(u,v) as an optimization expression, you must first convert it using fcn2optimexpr.

To use expfn3 as an optimization expression, first create optimization variables of the appropriate sizes.

u = optimvar('u',3,3,'LowerBound',-1,'UpperBound',1); % 3-by-3 variable
v = optimvar('v',3,'LowerBound',-2,'UpperBound',2); % 3-by-1 variable

Convert the function file to an optimization expressions using fcn2optimexpr.

[f,g,mineval] = fcn2optimexpr(@expfn3,u,v);

Because all returned expressions are scalar, you can save computing time by specifying the expression sizes using the 'OutputSize' name-value pair argument. Also, because expfn3 computes all of the outputs, you can save more computing time by using the ReuseEvaluation name-value pair.

[f,g,mineval] = fcn2optimexpr(@expfn3,u,v,'OutputSize',[1,1],'ReuseEvaluation',true)

f = 
  Nonlinear OptimizationExpression

    [argout,~,~] = expfn3(u, v)

g = 
  Nonlinear OptimizationExpression

    [~,argout,~] = expfn3(u, v)

mineval = 
  Nonlinear OptimizationExpression

    [~,~,argout] = expfn3(u, v)

Anonymous Function

To use a general nonlinear function handle in the problem-based approach, convert the handle to an optimization expression using fcn2optimexpr. For example, write a function handle equivalent to mineval and convert it.

fun = @(u)min(eig(u));
funexpr = fcn2optimexpr(fun,u,'OutputSize',[1,1])

funexpr = 
  Nonlinear OptimizationExpression

    anonymousFunction2(u)

  where:

    anonymousFunction2 = @(u)min(eig(u));

Create Objective

To use the objective expression as an objective function, create an optimization problem.

prob = optimproblem;
prob.Objective = f;

Define Constraints

Define the constraint g <= 0 in the optimization problem.

prob.Constraints.nlcons1 = g <= 0;

Also define the constraints that u is symmetric and that $$mineval\ge -1/2$$.

prob.Constraints.sym = u == u.';
prob.Constraints.mineval = mineval >= -1/2;

View the problem.

show(prob)

  OptimizationProblem : 

	Solve for:
       u, v

	minimize :
       [argout,~,~] = expfn3(u, v)


	subject to nlcons1:
       arg_LHS <= 0

       where:

         [~,arg_LHS,~] = expfn3(u, v);

	subject to sym:
       u(2, 1) - u(1, 2) == 0
       u(3, 1) - u(1, 3) == 0
       -u(2, 1) + u(1, 2) == 0
       u(3, 2) - u(2, 3) == 0
       -u(3, 1) + u(1, 3) == 0
       -u(3, 2) + u(2, 3) == 0

	subject to mineval:
       arg_LHS >= (-0.5)

       where:

         [~,~,arg_LHS] = expfn3(u, v);

	variable bounds:
       -1 <= u(1, 1) <= 1
       -1 <= u(2, 1) <= 1
       -1 <= u(3, 1) <= 1
       -1 <= u(1, 2) <= 1
       -1 <= u(2, 2) <= 1
       -1 <= u(3, 2) <= 1
       -1 <= u(1, 3) <= 1
       -1 <= u(2, 3) <= 1
       -1 <= u(3, 3) <= 1

       -2 <= v(1) <= 2
       -2 <= v(2) <= 2
       -2 <= v(3) <= 2

Solve Problem

To solve the problem, call solve. Set an initial point x0.

rng default % For reproducibility
x0.u = 0.25*randn(3);
x0.u = x0.u + x0.u.';
x0.v = 2*randn(3,1);
[sol,fval,exitflag,output] = solve(prob,x0)

Solving problem using fmincon.

Feasible point with lower objective function value found.


Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

sol = struct with fields:
    u: [3x3 double]
    v: [3x1 double]

fval = -403.4288

exitflag = 
    OptimalSolution

output = struct with fields:
         iterations: 74
          funcCount: 1156
    constrviolation: 0
           stepsize: 1.6235e-06
          algorithm: 'interior-point'
      firstorderopt: 1.1203e-04
       cgiterations: 122
            message: '...'
       bestfeasible: [1x1 struct]
             solver: 'fmincon'

View the solution.

disp(sol.u)

    0.5486    0.2067   -0.8420
    0.2067    0.1909    0.4842
   -0.8420    0.4842    0.8262

disp(sol.v)

    2.0000
   -2.0000
    2.0000

The solution matrix u is symmetric. All values of v are at the bounds.

Copyright 2018–2020 The MathWorks, Inc.