MATLAB Answers

Quadratically constrained linear maximisation problem: issues with fmincon

6 views (last 30 days)
CT
CT on 27 Mar 2020
Commented: Matt J on 31 Mar 2020
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 Comments

Matt J
Matt J on 31 Mar 2020
I suspect this should be due to my mistakes/inefficiences in providing gradient, hessian, etc.
I don't think that is the reason. I think the main reason is that Gurobi apparently has specific mechanisms for handling quadratic constraints, as suggested by this segment of code.
model.quadcon(1).Qc=Q;
model.quadcon(1).q=q;
model.quadcon(1).rhs=d;
Conversely, fmincon has no way of distinguishing quadratic constraints from other more general nonlinear constraints and giving them special handling. It handles all nonlinear constraints in the same way.
That said, the following implementation of your functions would be slightly more efficient.
function [y,yeq,grady,gradyeq] = quadconstr(z,Q,d)
Qz=Q*z;
y= z'*Qz - d;
yeq = []; %no quadratic inequalities
if nargout > 2
grady = 2*Qz;
gradyeq = []; %no quadratic inequalities
end
end
function [y,grady] = quadobj(z,f)
grady = -f;
y = grady.'*z;
end
Matt J
Matt J on 31 Mar 2020
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.

Sign in to comment.

Answers (1)

Matt J
Matt J on 30 Mar 2020
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 Comments

Show 7 older comments
Matt J
Matt J on 31 Mar 2020
How fast does Gurobi solve the problem when the quadratic inequality is replaced with equality? And how fast when the quadratic constraint is removed altogether?
CT
CT on 31 Mar 2020
1) With quadratic equality: I need to investigate because it is non-convex
2) Without quadratic constraint: 0.16 sec.
Matt J
Matt J on 31 Mar 2020
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?

Sign in to comment.

Sign in to answer this question.