Derivative of Function handle (To obtain the derivative of constraints in function fmincon)

Hey, everyone, I'm using fmincon to solve the optimization problem under constraints. To increase the speed of calculation, I would like to provide the derivative matrix of the constraints. However, the constraint conditions are nonliner and complex, e.g as below.

where

the last one implies is also a function of . If I want to calculate the derivative of respect to , it will take a lot of efforts for me to calculate it by myself.

Does the MATLAB provide some tool to obtain the derivative of function handle, especially under using fmincon to solve the optimization problem ? Any suggestions to obtain the derivative of constraint conditions efficiently instead of calculating them by hand?

Answer 1

[Matt J]

You can use the Symbolic Math Toolbox to obtain the expression for the derivative, and then use matlabFunction() to convert from symbolic form to function form.

[sxj]

Hey, Matt, thank you for your message and I tried your suggestions. The code below is a simple example where I use global to define parameters in the symbolic expression and returns the value of constraint derivative at the endogenous variabes to the main function. Am I right with the example?

function [output] = Derivative(v1,v2,v3)

% Input: v1,v2,v3 are the endogenous variables of optimization problem

% Parameters: a b c are parameters defined in the main function

% Output: the vector of derivative of constraint conditions respect to v1,v2,v3

global a b c

a = 3;
b = 4;
c = 5;

syms x y z 

h1(x,y,z) = z^2;
f1(x,y,z)=a*x^2+b*y^2+c*h1(x,y,z)^2;      % the constraint condition in the optimization problem
g1 = matlabFunction(diff(f1,x));
g2 = matlabFunction(diff(f1,y));
g3 = matlabFunction(diff(f1,z));

d11 = g1(v1,v2,v3);
d12 = g2(v1,v2,v3);
d13 = g3(v1,v2,v3);

output = [d11,d12,d13]; 

end

[Torsten]

You shouldn't build g1,g2 and g3 each time the optimizer needs to get the derivatives for certain values v1, v2 and v3.

Build g1, g2 and g3 only once before the start of the optimization, pass them to the optimizer and evaluate them for the v-values given.

[Matt J]

A few more remarks:

First, you should just use gradient or jacobian instead of diff to avoid generating 3 separate functions. Second, because you are interested in speed, I would write the output to a file, which can be more optimally auto-coded. As Torsten said, all of this should be done before you run fmincon.

syms a b c
v=sym('v',[3,1]);

f=[a,b,c]*v.^2; 

matlabFunction(f, gradient(f,v),'File','myobjective.m','Optimize',true)

This example generates a file myobjective.m that looks like this, which is almost entirely the form you need for fmincon. Notice that I made sure that a, b, and c would be arguments to this external file.

function [f,out2] = myobjective(a,b,c,v1,v2,v3)
%MYOBJECTIVE
%    [F,OUT2] = MYOBJECTIVE(A,B,C,V1,V2,V3)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    25-Jun-2019 07:40:24

f = a.*v1.^2+b.*v2.^2+c.*v3.^2;
if nargout > 1
    out2 = [a.*v1.*2.0;b.*v2.*2.0;c.*v3.*2.0];
end

Personally, I would now edit this slightly to accept vector-valued input, which is what fmincon expects,

function [f,out2] = myobjective(v,knowns)

     v=num2cell(v);
[v1,v2,v3]=deal(v{:});

    knowns=num2cell(knowns);
[a,b,c]=deal(knowns{:});

f = a.*v1.^2+b.*v2.^2+c.*v3.^2;
if nargout > 1
    out2 = [a.*v1.*2.0;b.*v2.*2.0;c.*v3.*2.0];
end

Now you are ready to run fmincon as follows

knowns=[3,4,5];
v_opt=fmincon(@(v) myobjective(v,knowns),  v0,____)

[sxj]

Thank you very much for the details, Torsten and Matt, and your suggestions are pretty helpful. I change my code according to your suggestions, but my optimization problem is very complex, so I list the simple version of my problem and relevant code of using fmincon here. Further comments are appreciated.

The optimization problem:

The main code

global C1 C2 a b

%Parameters
C1 = 10;
C2 = 20;
a  = 2;
b  = 4;

%Initial Guess
x_0 = [20;10;0.5]; % I choose the guess values randomly

% Using fmincon to solve the optimization problem 
tol = 1.0E-13;
options = optimset( ...
    'Display',         'off', ...
    'GradObj',         'on', ...
    'GradConstr',      'on', ...
    'DerivativeCheck', 'off', ...
    'FinDiffType',     'central', ...
    'TolFun',          tol, ...
    'TolX',            tol, ...
    'TolCon',          tol, ...
    'algorithm',       'active-set', ...
    'MaxFunEvals',     inf, ...
    'MaxIter',         5000);

lb = -10.*x_0;
ub =  10.*x_0;

opt = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)constraint(x), options);

The objective function

function [obj,grad] = objective(x)
obj = -x(1,1);  
     if nargout>1
       grad=zeros(size(x)); 
        grad(1,1)=-1;
     end
end

The constraint function

function [c, ceq, GC, GCeq] = constraint(x)

global C1 C2 a b

Y  = x(1,1);
x1 = x(1,2);
x2 = x(1,3);

% Constraint conditions 
COND1 = C1 - a*x1^(2)-b*x2^(2);
COND2 = C2 - a*x1^(3)-2*b*x2;
COND3 = x1/x2;

CSTRT1 = [COND1; COND2];
CSTRT2 = Y - COND3;

ceq  = [CSTRT1; CSTRT2];

%The derivative of constraint conditions
GCeq = derivative_cons(Y,x1,x2);

%The other outputs
c    = [];
GC   = [];

end

The derivative function of constraint conditions

function [output] = derivative_cons(v1,v2,v3)

global C1 C2 a b 

syms Y x1 x2
z = [Y,x1,x2];
f1(z) = C1 - a*x1^(2)-b*x2^(2);
f2(z) = C2 - a*x1^(3)-2*b*x2;
f3(z) = Y  - x1/x2;

grad1 = matlabFunction(gradient(f1,z));
grad2 = matlabFunction(gradient(f2,z));
grad3 = matlabFunction(gradient(f3,z));

output = [grad1(v1,v2,v3);grad2(v1,v2,v3);grad3(v1,v2,v3)];

end

[Matt J]

No. Again, you should not be calling matlabFunction inside derivative_cons or inside any other function that is called iteratively during the optimization. matlabFunction is a code writing tool. You should be using it to write derivative_cons and you would do this before running the actual optimization. Also, you should never have to call matlabFunction more than once. Use the jacobian command as in the code below to do all derivatives in one step.

syms a b c
syms Y x1 x2
z = [Y,x1,x2];
f1(z) = C1 - a*x1^(2)-b*x2^(2);
f2(z) =  C2 - a*x1^(3)-2*b*x2;
f3(z) =  Y  - x1/x2;
output = matlabFunction(jacobian([f1;f2;f3],z).','File','derivative_cons');

Finally, you should try to avoid using global variables. They are slow and dangerous. Instead, use one of the other methods discussed here,

https://www.mathworks.com/help/optim/ug/passing-extra-parameters.html

[sxj]

Hey, Matt, thank you very much for these further comments. I adjust my code for the simple example as below, including using matlabFunction out of loops and avoid using global variables. But when I run the code, the error message is 'Error using derivative_cons

Too many output arguments.' I couldn't figure out how to fix the output issue of constraint and its gradient condition. Your comments are appreciated.

The main code

% Parameters
 C1 = 12;
 C2 = 20;
 a  = 2;
 b  = 4;

syms Y x1 x2
z = [Y,x1,x2];
f1(z) =  C1 - a*x1^(2)-b*x2^(2);
f2(z) =  C2- a*x1^(3)-2*b*x2;
f3(z) =  Y  - x1/x2;
cons = matlabFunction([f1;f2;f3],jacobian([f1;f2;f3],z).','File','derivative_cons');

%Use 'fmincon' to solve out the optimization problem
x_0 = [20;10;0.5]; % I choose the guess values randomly
tol = 1.0E-13;
options = optimset( ...
    'Display',         'off', ...
    'GradObj',         'on', ...
    'GradConstr',      'on', ...
    'DerivativeCheck', 'off', ...
    'FinDiffType',     'central', ...
    'TolFun',          tol, ...
    'TolX',            tol, ...
    'TolCon',          tol, ...
    'algorithm',       'active-set', ...
    'MaxFunEvals',     inf, ...
    'MaxIter',         5000);
lb = -10.*x_0;
ub =  10.*x_0;
opt = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)cons(x), options);

The objective function

function [obj,grad] = objective(x)
obj = -x(1,1);  
     if nargout>1
       grad=zeros(size(x)); 
        grad(1,1)=-1;
     end
end

[Walter Roberson]

With regards to vector valued inputs: if you use the 'vars' option of matlabFunction and pass it a cell array that contains a vector of variable names, then it will bundle those all into a vector for calling purposes.

[Walter Roberson]

cons = matlabFunction([f1;f2;f3],jacobian([f1;f2;f3],z).','File','derivative_cons');

Should not have the Jacobian second parameter.

[sxj]

Hey, Walter, thank you for the quick reply. I transfer the code to be

cons = matlabFunction([f1;f2;f3],jacobian([f1;f2;f3]).','File','derivative_cons');

as you said, but MATLAB still provide the same error message. Do I misunderstand what you mean?

According to my understanding, the 'cons' contains two function: '[f1;f2;f3]' represents the constraint conditions, while 'jacobian([f1;f2;f3],z).'' represents the gradient of constraints. These two function in 'cons' are consistent with the requirement of 'fmincon'.

[Walter Roberson]

I had not realized that you could output constraints on the gradients. Looking at the documentation, when you use that option, you must return four values:

If the gradients of the constraints can also be computed and the SpecifyConstraintGradient option is true, as set by

options = optimoptions('fmincon','SpecifyConstraintGradient',true)

then nonlcon must also return, in the third and fourth output arguments, GC, the gradient of c(x), and GCeq, the gradient of ceq(x). GC and GCeq can be sparse or dense. If GC or GCeq is large, with relatively few nonzero entries, save running time and memory in the interior-point algorithm by representing them as sparse matrices. For more information, see Nonlinear Constraints.

Your matlabFunction call would have to have four arguments before the options (and you will need to use the 'vars' option the way I described earlier.)

[Matt J]

@sxj, could you show the full error message and also the full contents of derivative_cons.m?

[sxj]

To Walter and Matt, according to your suggestions, I change the relevant codes as below

cons = matlabFunction([],[f1;f2;f3],[],jacobian([f1;f2;f3]).','File','derivative_cons');
opt = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)cons(x), options);

where I use [] to assign empty matrix to relevant condition about 'c' and 'GC'. In this way, the error 'Too many output argument' disappear but the following error message appears:

I tried the following code to fix it but the error still exists

opt = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(Y,x1,x2)cons(Y,x1,x2), options);

how many inputs do I need here? the full message of derivative_cons.m is as below:

[OUT1,OUT2,OUT3,OUT4] = derivative_cons(Y,x1,x2)

which comes from the following code

syms Y x1 x2
z = [Y,x1,x2];
f1(z) =  C1 - a*x1^(2)-b*x2^(2);
f2(z) =  C2- a*x1^(3)-2*b*x2;
f3(z) =  Y  - x1/x2;
cons = matlabFunction([],[f1;f2;f3],[],jacobian([f1;f2;f3]).','File','derivative_cons');

[Matt J]

Remember, fmincon will pass inputs to (and expect outputs from) your constraint function in vector form. However, derivative_cons expects each x(i) as a separate argument. You must therefore do things like this,

@(x)cons(x(1),x(2),x(3))

[sxj]

thank you for the quick reply, Matt. The error of "not enough input argument' disappear now. But the code

opt = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)cons(Y,x1,x2), options);

makes the new error message:

I think 'matlabFunction' has already transfer the constraint condition and its derivatives to be the function handle, but why this error comes out here if I provide inputs here.

(Sorry, I use 'matlabFunction' and 'fmincon' for the first time so I'm not familiar with their characteristics)

[Matt J]

It's hard to say, because you haven't shown the complete error message. Incidentally, I encourage you not to post screen shots of error text. Copy/pasing the text will consume much less disk space in the forum.

In any case, why not take the message you have shown at face value? Test the constraint and objective functions separately and see if they are returning type double or not.

[Walter Roberson]

@(x)cons(x(1),x(2),x(3))

You do not need to do this if you use the vars option to matlabFunction like I have been recommending.

[Walter Roberson]

@(x)cons(Y,x1,x2)

Why are you ignoring the current state of minimization from the x input and calculating the constraints based upon the symbolic variables Y x1 x2?

[sxj]

Thanks for your help, Matt and Walter, it's my mistake and I have fix it now. The code works well without errors.

But for your suggestion using vars option in matlabFunction, I tried the following ways, but all of them still report the error message 'not enough input arguments'.

cons = matlabFunction([],[f1;f2;f3],[],jacobian([f1;f2;f3]).','File','derivative_cons1','Vars',[Y,x1,x2]);
opt = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)cons(x), options);

or

w = {'Y','x1','x2'};
cons = matlabFunction([],[f1;f2;f3],[],jacobian([f1;f2;f3]).','File','derivative_cons1','Vars',w);
opt = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)cons(x), options);

or

v = {Y,x1,x2};
cons = matlabFunction([],[f1;f2;f3],[],jacobian([f1;f2;f3]).','File','derivative_cons1','Vars',v);
opt = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)cons(x), options);

Still need to use @(x)cons(x(1),x(2),x(3)) in these situations. Do I misunderstand your words?

[Walter Roberson]

As I wrote above,

"if you use the 'vars' option of matlabFunction and pass it a cell array that contains a vector of variable names, then it will bundle those all into a vector for calling purposes."

You tried passing a vector of values. You tried passing a cell array of values. You did not try passing a cell array containing a vector.

cons = matlabFunction([],[f1;f2;f3],[],jacobian([f1;f2;f3]).','File','derivative_cons1','Vars',{[Y,x1,x2]});

[sxj]

Hey, Walter, I tried your suggestions

cons = matlabFunction([],[f1;f2;f3],[],jacobian([f1;f2;f3]).','File','derivative_cons1','Vars',{[Y,x1,x2]});
opt  = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)cons(x), options);

But the following error message appears:

Index in position 2 exceeds array bounds (must not exceed 1).
Error in derivative_cons1 (line 9)
x1 = in1(:,2);
Error in @(x)cons(x)
Error in fmincon (line 656)
      [ctmp,ceqtmp,initVals.gnc,initVals.gnceq] = feval(confcn{3},X,varargin{:});
Caused by:
    Failure in initial nonlinear constraint function evaluation. FMINCON cannot continue.

where the derivative_cons1 from matlabFunction is

function [out1,out2,out3,out4] = derivative_cons1(in1)
%DERIVATIVE_CONS1
%    [OUT1,OUT2,OUT3,OUT4] = DERIVATIVE_CONS1(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.3.
%    28-Jun-2019 15:59:37

Y = in1(:,1);
x1 = in1(:,2);
x2 = in1(:,3);
out1 = zeros(0,0);
if nargout > 1
    t2 = x1.^2;
    t3 = x2.*8.0;
    t4 = 1.0./x2;
    t5 = -t3;
    out2 = [t2.*-2.0-1.0./t4.^2.*4.0+1.2e+1;t5-x1.^3.*2.0+2.0e+1;Y-t4.*x1];
end
if nargout > 2
    out3 = zeros(0,0);
end
if nargout > 3
    out4 = reshape([0.0,x1.*-4.0,t5,0.0,t2.*-6.0,-8.0,1.0,-t4,t4.^2.*x1],[3,3]);
end

Do I misunderstand your words?

[Matt J]

Change the first few lines in derivative_cons1,

Y = in1(1);
x1 = in1(2);
x2 = in1(3);

or make x_0 a row vector,

x_0 = [20,10,0.5]; 

[sxj]

Thanks for your quick reply, Matt. I follow your suggestion and it works well. But I think it is inconvenient to change derivative_cons1 during running the whole code. So make x_0 a row vector will be better.

I also tried the following way

cons = matlabFunction([],[f1;f2;f3],[],jacobian([f1;f2;f3]).','File','derivative_cons1','Vars',{[Y;x1;x2]});
x_0 = [20;10;0.5];
opt  = fmincon(@(x)objective(x), x_0,[], [], [],[],lb,ub, @(x)cons(x), options);

where I change the vector in the option 'Vars' from 1X3 to 3X1. In this way, it also works.

So far I have solved the simple version of my optimization problem and will turn to my original one. Thank you very much for all the discussions above, especially to Matt and Walter. I hope the discussions will also help others.