Quadratically constrained linear maximisation problem: issues with fmincon

7 次查看(过去 30 天)
I would like to solve the following quadratically constrained linear programming problem.
I have written a Matlab code (R2091b) that solve the problem using Gurobi. Now, I would like to rewrite the code using fmincon instead of Gurobi. This is because the optimisation problem will have to be solve thousands of times and Gurobi's academic license does not allow to parallelise via array jobs in a cluster. However, I'm encountering a huge problem: Gurobi takes 0.23 second to give a solution, fmincon takes 13 sec. I suspect this should be due to my mistakes/inefficiences in providing gradient, hessian, etc. Could you kindly help me to improve below?
Also, Gurobi gives me 0.2 as solution, fmincon gives 0.089. Can the accuracy of fmincon be improved without trying other starting points?
This is my code with Gurobi
clear
rng default
load matrices
%1) GUROBI
model.A=[Aineq; Aeq];
model.obj=f;
model.modelsense='max';
model.sense=[repmat('<', size(Aineq,1),1); repmat('=', size(Aeq,1),1)];
model.rhs=[bineq; beq];
model.ub=ub;
model.lb=lb;
model.quadcon(1).Qc=Q;
model.quadcon(1).q=q;
model.quadcon(1).rhs=d;
params.outputflag = 0;
result=gurobi(model, params);
max_problem_Gurobi=result.objval;
This is my code with fmincon
%2) FMINCON
options = optimoptions(@fmincon,'Algorithm','interior-point',...
'SpecifyObjectiveGradient',true,...
'SpecifyConstraintGradient',true,...
'HessianFcn',@(z,lambda)quadhess(z,lambda,Q));
fun = @(z)quadobj(z,f.');
nonlconstr = @(z)quadconstr(z,Q,d);
[~,fval] = fmincon(fun,z0,Aineq,bineq,Aeq,beq,lb,ub,nonlconstr,options);
max_problem_fmincon=-fval;
function [y,yeq,grady,gradyeq] = quadconstr(z,Q,d)
y= z'*Q*z -d;
yeq = []; %no quadratic inequalities
if nargout > 2
grady = 2*Q*z;
end
gradyeq = []; %no quadratic inequalities
end
function hess = quadhess(z,lambda,Q) %#ok<INUSL>
hess = 2*lambda.ineqnonlin*Q;
end
function [y,grady] = quadobj(z,f)
y = -f'*z;
if nargout > 1
grady = -f;
end
Further, if I run the code with fmincon with as starting point the optimal point given by Gurobi, I still get the solution 0.089 (instead of 0.2 as in Gurobi). Why?
%3) FMINCON
z0=result_Gurobi.x;
options = optimoptions(@fmincon,'Algorithm','interior-point',...
'SpecifyObjectiveGradient',true,...
'SpecifyConstraintGradient',true,...
'HessianFcn',@(z,lambda)quadhess(z,lambda,Q));
fun = @(z)quadobj(z,f.');
nonlconstr = @(z)quadconstr(z,Q,d);
[~,fval] = fmincon(fun,z0,Aineq,bineq,Aeq,beq,lb,ub,nonlconstr,options);
max_problem_fmincon=-fval;
function [y,yeq,grady,gradyeq] = quadconstr(z,Q,d)
y= z'*Q*z -d;
yeq = []; %no quadratic inequalities
if nargout > 2
grady = 2*Q*z;
end
gradyeq = []; %no quadratic inequalities
end
function hess = quadhess(z,lambda,Q) %#ok<INUSL>
hess = 2*lambda.ineqnonlin*Q;
end
function [y,grady] = quadobj(z,f)
y = -f'*z;
if nargout > 1
grady = -f;
end
end
  18 个评论
Matt J
Matt J 2020-3-31
编辑:Matt J 2020-3-31
There was no expectation of great gains, but I think it has to be marginally more efficient, maybe a reduction from 36.05 sec to 36 sec. You can plainly see that fewer vector arithmetic operations are done in this version.

请先登录,再进行评论。

回答(1 个)

Matt J
Matt J 2020-3-30
Well, it would be interesting to know what algorithm Gurobi uses, but the issue of the objective function difference appears to be a matter of the tolerances
options = optimoptions(@fmincon,'Algorithm','interior-point',...
'SpecifyObjectiveGradient',true,...
'SpecifyConstraintGradient',true,...
'HessianFcn',@(z,lambda)quadhess(z,lambda,Q),...
'StepTolerance',1e-30,'OptimalityTolerance',1e-10);
fun = @(z)quadobj(z,f);
nonlconstr = @(z)quadconstr(z,Q,d);
tic;
[~,fval] = fmincon(fun,z0(:),Aineq,bineq,Aeq,beq,lb,ub,nonlconstr,options);
toc
max_problem_fmincon=-fval
max_problem_fmincon =
0.2000
  10 个评论
CT
CT 2020-3-31
1) With quadratic equality: I need to investigate because it is non-convex
2) Without quadratic constraint: 0.16 sec.
Matt J
Matt J 2020-3-31
And does the problem data from the thousands of problem instances that you are trying to solve change in a continuous incremental way? If you had the optimal solution for one instance of the problem, would it serve as a good initial estimate for the next instance?

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Quadratic Programming and Cone Programming 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by