fmincon pause after 2 iterations when using gradients

4 visualizzazioni (ultimi 30 giorni)
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
  3 Commenti
Hendrik van Gils
Hendrik van Gils il 8 Mar 2022
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
Matt J il 8 Mar 2022
Have you tried running with "Pause on Nan or Inf"

Accedi per commentare.

Risposta accettata

Alan Weiss
Alan Weiss il 9 Mar 2022
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

Più risposte (1)

Matt J
Matt J il 8 Mar 2022
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.
  1 Commento
Hendrik van Gils
Hendrik van Gils il 8 Mar 2022
Sorry if my formulation was vague, I generate the function once, before the optimization, and then just reuse it.

Accedi per commentare.

Prodotti


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by