fmincon pause after 2 iterations when using gradients

Without gradients my code runs 'fine'. Meaning that it is slow to converge, but at least it is continuously iterating towards a better solution.

Adding gradients results in a freeze/pauze, usually after the second iteration. After a while it picks up again and it's better/faster than without gradients. Consecutive runs of the same script don't have this freeze/pause, which leads me to believe it's some sort of caching related issue.

'checkGradients' tells me that the gradients are good. I get my gradients from a different script by taking the symbolic jacobian of ceq, and using matlabFunction to generate a sparse optimized matlab function file.

If anyone could explain why this happens and maybe how I could resolve the slowdown I would be very grateful!

(the slowdown is not because of the 'checkGradients' as it also happens without the check. I only generate the function once, not every time I try to optimise.)

CheckGradients Information

Objective function derivatives:

Maximum relative difference between supplied

and finite-difference derivatives = 2.0994e-07.

Nonlinear equality constraint derivatives:

Maximum relative difference between supplied

and finite-difference derivatives = 4.20079e-08.

CheckGradients successfully passed.

____________________________________________________________

First-order Norm of

Iter F-count f(x) Feasibility optimality step

0 1 1.960200e+01 9.900e-01 5.822e-01

1 2 2.011323e+01 9.888e-01 1.048e+00 7.881e-01

2 3 2.174207e+01 9.884e-01 1.435e+00 4.568e-01

[Matt J]

We'll need to see code, preferably in a form that we can run ourselves.

[Hendrik van Gils]

My script is dependent on tens of files, so sadly sending it in a way that you can run yourself isn't the most viable option. But here are some relevant pieces of code.

Where and how I call fmincon:

%% fmincon stuff

options = optimoptions('fmincon',...

'Algorithm','interior-point',...

'Display','iter',...

'MaxFunEvals',2e5,'MaxIterations',1e3,...

'SpecifyConstraintGradient',true,...

'SpecifyObjectiveGradient',true,...

'CheckGradients',false);

% 'OutputFcn', @(z,optimValues,state) outputFunction(z,optimValues,state,@unpackVariables,pack,figures),...

problem.x0 = zGuess;

problem.lb = zLowerBound;

problem.ub = zUpperBound;

problem.Aineq = Aineq; problem.bineq = bineq;

problem.Aeq = Aeq; problem.beq = beq;

problem.options = options;

problem.solver = 'fmincon';

problem.objective = @(z) objectiveFunction(z,pack,objective); % Kelly uses weights for trapezoidal integration, why?

problem.nonlcon = @(z) constraintFunction(@unpackVariables,z,pack,dynamics);

[zSolution, cost, exitFlag, output] = fmincon(problem);

The constraint function:

function [c,ceq,gradc,gradceq] = constraintFunction(unpackFunction,z,pack,dynamics)

% This constraint function is meant for the NONLINEAR constraints ONLY.

% Therefore only the second half of the state vector is used, as it is

% nonlinear. The first half will be used in the main script to formulate

% efficient linear constraints.

[t,x,u] = unpackFunction(z,pack); % Unpack variables

xHalf1 = 1:size(x,1)/2; % Index upper half state vector

xHalf2 = xHalf1+size(x,1)/2; % Index lower half state vector

dt = t(:,2:end) - t(:,1:end-1);

dx = x(xHalf2,2:end) - x(xHalf2,1:end-1);

f = dynamics(t,x,u);

% Specify inequality constraints

c = [];

% Specify equality constraints

defects = dx - 0.5*dt.*(f(xHalf2,1:end-1) + f(xHalf2,2:end)); % Trapazoid rule, should equal 0

ceq = defects(:);

%% If Gradient is required

if nargout > 2

gradc = gradientInequalityConstraints(z);

gradceq = gradientEqualityConstraints(z);

end

end

How the gradient function is made:

%% Constraint function gradient

xHalf1 = 1:size(x,1)/2; % Index upper half state vector

xHalf2 = xHalf1+size(x,1)/2; % Index lower half state vector

dt = (t(end) - t(1))/(nCollocationPoints - 1);

dx = x(xHalf2,2:end) - x(xHalf2,1:end-1);

f = dynamics(t,x,u);

% Specify inequality constraints

c = [];

% Specify equality constraints

defects = dx - 0.5*dt.*(f(xHalf2,1:end-1) + f(xHalf2,2:end)); % Trapazoid rule, should equal 0

ceq = defects(:);

% Calculate gradient

gradientInequalityConstraints = jacobian(c,z).';

gradientEqualityConstraints = jacobian(ceq,z).';

if generateConstraintGradients

matlabFunction(gradientInequalityConstraints,'File','gradientInequalityConstraints','Vars',{z},'Sparse',true)

matlabFunction(gradientEqualityConstraints,'File','gradientEqualityConstraints','Vars',{z},'Sparse',true)

end

I hope this helps!

[Matt J]

Have you tried running with "Pause on Nan or Inf"

Answer 1

[Alan Weiss]

I don't know for sure, but MATLAB has a just-in-time compiler that does indeed make a sort of cache for code that is repeatedly called. On first calling a function, the compiler makes a compiled version that is then reused on subsequent calls. Your matlabFunction call probably creates a long, complex function that is slow to interpret the first time, and is fast to run thereafter.

Or I could be wrong, and what is happening is your matlabFunction call creates a call to an obscure portion of the code base that takes a while to load. Sorry, I don't really know which is more likely, or how you would check it. But perhaps running your code for an iteration or two before calling your optimization will smooth things out either way.

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 2

[Matt J]

I get my gradients from a different script by taking the symbolic jacobian of ceq, and using matlabFunction to generate a sparse optimized matlab function file.

If you're repeating this process every time your constraint function is called, it's probably a bad idea. You should generate your matlabFunction once, before the optimization is launched.

[Hendrik van Gils]

Sorry if my formulation was vague, I generate the function once, before the optimization, and then just reuse it.