solveCobraQP

PURPOSE ^

solveCobraQP Solve constraint-based QP problems

SYNOPSIS ^

function solution = solveCobraQP(QPproblem,varargin)

DESCRIPTION ^

solveCobraQP Solve constraint-based QP problems

 solution = solveCobraQP(QPproblem,parameters)

 % Solves problems of the type 

      min   0.5 x' * F * x + osense * c' * x
      s/t   lb <= x <= ub
            A * x  <=/=/>= b

INPUT
 QPproblem Structure containing the following fields describing the QP
 problem to be solved
  A      LHS matrix
  b      RHS vector
  F      F matrix for quadratic objective (must be positive definite)
  c      Objective coeff vector
  lb     Lower bound vector
  ub     Upper bound vector
  osense Objective sense (-1 max, +1 min)
  csense Constraint senses, a string containting the constraint sense for
         each row in A ('E', equality, 'G' greater than, 'L' less than).

OPTIONAL INPUTS
 Optional parameters can be entered using parameters structure or as
 parameter followed by parameter value: i.e. ,'printLevel',3)

 parameters    Structure containing optional parameters as fields.
  printLevel   Print level for solver
  saveInput    Saves LPproblem to filename specified in field. 
               Setting parameters = 'default' uses default setting set in
               getCobraSolverParameters.

 The solver is defined in the CBT_MILP_SOLVER global variable
 (set using changeCobraSolver). Solvers currently available are
 'tomlab_cplex', 'mosek' and 'qpng' (limited support for small problems)

OUTPUT
 solution  Structure containing the following fields describing a QP
           solution
  full     Full QP solution vector
  obj      Objective value
  solver   Solver used to solve QP problem
  stat     Solver status in standardized form (see below)
  origStat Original status returned by the specific solver 
  time     Solve time in seconds

  stat     Solver status in standardized form
           1   Optimal solution
           2   Unbounded solution
           0   Infeasible
           -1  No solution reported (timelimit, numerical problem etc)

 Markus Herrgard        6/8/07
 Ronan Fleming         12/07/09  Added support for mosek
 Ronan Fleming         18 Jan 10 Added support for pdco
 Josh Lerman           04/17/10 changed def. parameters, THREADS, QPMETHOD
 Tim Harrington        05/18/12 Added support for the Gurobi 5.0 solver

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function solution = solveCobraQP(QPproblem,varargin)
0002 %solveCobraQP Solve constraint-based QP problems
0003 %
0004 % solution = solveCobraQP(QPproblem,parameters)
0005 %
0006 % % Solves problems of the type
0007 %
0008 %      min   0.5 x' * F * x + osense * c' * x
0009 %      s/t   lb <= x <= ub
0010 %            A * x  <=/=/>= b
0011 %
0012 %INPUT
0013 % QPproblem Structure containing the following fields describing the QP
0014 % problem to be solved
0015 %  A      LHS matrix
0016 %  b      RHS vector
0017 %  F      F matrix for quadratic objective (must be positive definite)
0018 %  c      Objective coeff vector
0019 %  lb     Lower bound vector
0020 %  ub     Upper bound vector
0021 %  osense Objective sense (-1 max, +1 min)
0022 %  csense Constraint senses, a string containting the constraint sense for
0023 %         each row in A ('E', equality, 'G' greater than, 'L' less than).
0024 %
0025 %OPTIONAL INPUTS
0026 % Optional parameters can be entered using parameters structure or as
0027 % parameter followed by parameter value: i.e. ,'printLevel',3)
0028 %
0029 % parameters    Structure containing optional parameters as fields.
0030 %  printLevel   Print level for solver
0031 %  saveInput    Saves LPproblem to filename specified in field.
0032 %               Setting parameters = 'default' uses default setting set in
0033 %               getCobraSolverParameters.
0034 %
0035 % The solver is defined in the CBT_MILP_SOLVER global variable
0036 % (set using changeCobraSolver). Solvers currently available are
0037 % 'tomlab_cplex', 'mosek' and 'qpng' (limited support for small problems)
0038 %
0039 %OUTPUT
0040 % solution  Structure containing the following fields describing a QP
0041 %           solution
0042 %  full     Full QP solution vector
0043 %  obj      Objective value
0044 %  solver   Solver used to solve QP problem
0045 %  stat     Solver status in standardized form (see below)
0046 %  origStat Original status returned by the specific solver
0047 %  time     Solve time in seconds
0048 %
0049 %  stat     Solver status in standardized form
0050 %           1   Optimal solution
0051 %           2   Unbounded solution
0052 %           0   Infeasible
0053 %           -1  No solution reported (timelimit, numerical problem etc)
0054 %
0055 % Markus Herrgard        6/8/07
0056 % Ronan Fleming         12/07/09  Added support for mosek
0057 % Ronan Fleming         18 Jan 10 Added support for pdco
0058 % Josh Lerman           04/17/10 changed def. parameters, THREADS, QPMETHOD
0059 % Tim Harrington        05/18/12 Added support for the Gurobi 5.0 solver
0060 
0061 global CBT_QP_SOLVER;
0062 
0063 if (~isempty(CBT_QP_SOLVER))
0064     solver = CBT_QP_SOLVER;
0065 else
0066     error('No solver found');
0067 end
0068 
0069 optParamNames = {'printLevel','saveInput'};
0070 parameters = '';
0071 if nargin ~=1
0072     if mod(length(varargin),2)==0
0073         for i=1:2:length(varargin)-1
0074             if ismember(varargin{i},optParamNames)
0075                 parameters.(varargin{i}) = varargin{i+1};
0076             else
0077                 error([varargin{i} ' is not a valid optional parameter']);
0078             end
0079         end
0080     elseif strcmp(varargin{1},'default')
0081         parameters = 'default';
0082     elseif isstruct(varargin{1})
0083         parameters = varargin{1};
0084     else
0085         display('Warning: Invalid number of parameters/values')
0086         solution=[];
0087         return;
0088     end
0089 end
0090 
0091 % Defaults in case the solver does not return anything
0092 x = [];
0093 y = [];
0094 w = [];
0095 f = [];
0096 xInt = [];
0097 xCont = [];
0098 stat = -99;
0099 solStat = -99;
0100 
0101 %parameters
0102 [printLevel, saveInput] = getCobraSolverParams('QP',optParamNames,parameters);
0103 
0104 [A,b,F,c,lb,ub,csense,osense] = ...
0105     deal(QPproblem.A,QPproblem.b,QPproblem.F,QPproblem.c,QPproblem.lb,QPproblem.ub,...
0106     QPproblem.csense,QPproblem.osense);
0107 
0108 t_start = clock;
0109 switch solver
0110 %%
0111     case 'tomlab_cplex'
0112         if (~isempty(csense))
0113             b_L(csense == 'E') = b(csense == 'E');
0114             b_U(csense == 'E') = b(csense == 'E');
0115             b_L(csense == 'G') = b(csense == 'G');
0116             b_U(csense == 'G') = inf;
0117             b_L(csense == 'L') = -inf;
0118             b_U(csense == 'L') = b(csense == 'L');
0119         else
0120             b_L = b;
0121             b_U = b;
0122         end
0123         tomlabProblem = qpAssign(F,osense*c,A,b_L,b_U,lb,ub,[],'CobraQP');
0124         
0125         %optional parameters
0126         tomlabProblem.PriLvl=printLevel;
0127         tomlabProblem.MIP.cpxControl.QPMETHOD = 1;
0128         tomlabProblem.MIP.cpxControl.THREADS = 1;
0129                 
0130         %Save Input if selected
0131         if ~isempty(saveInput)
0132             fileName = saveInput;
0133             if ~find(regexp(fileName,'.mat'))
0134                 fileName = [fileName '.mat'];
0135             end
0136             display(['Saving QPproblem in ' fileName]);
0137             save(fileName,'QPproblem')
0138         end
0139         
0140         Result = tomRun('cplex', tomlabProblem);
0141         x = Result.x_k;
0142         f = osense*Result.f_k;
0143         origStat = Result.Inform;
0144         if (origStat == 1)
0145             stat = 1; % Optimal
0146         elseif (origStat == 3 || origStat == 4)
0147             stat = 0; % Infeasible
0148         elseif (origStat == 2)
0149             stat = 2; % Unbounded
0150         elseif (origStat >= 10)
0151             stat = -1; % No optimal solution found (time or other limits reached, other infeasibility problems)
0152         else
0153             stat = 3; % Solution exists, but either scaling problems or not proven to be optimal
0154         end
0155         %%
0156     case 'cplex_direct'
0157         %% Tomlab cplex.m direct
0158         %Used with the current script, only some of the control affoarded with
0159         %this interface is provided. Primarily, this is to change the print
0160         %level and whether to minimise the Euclidean Norm of the internal
0161         %fluxes or not.
0162         %See solveCobraLPCPLEX.m for more refined control of cplex
0163         %Ronan Fleming 11/12/2008
0164                
0165         solution=solveCobraLPCPLEX(QPproblem,printLevel,[],[],[],minNorm);
0166         %%
0167    case 'qpng'
0168         % qpng.m This file is part of GLPKMEX.
0169         % Copyright 2006-2007 Nicolo Giorgetti.
0170         %
0171         % Solves the general quadratic program
0172         %      min 0.5 x'*H*x + q'*x
0173         %       x
0174         % subject to
0175         %      A*x [ "=" | "<=" | ">=" ] b
0176         %      lb <= x <= ub
0177         ctype=csense;
0178         ctype(('G'==csense))='L';
0179         ctype(('E'==csense))='E';
0180         ctype(('L'==csense))='U';
0181         
0182         x0=ones(size(QPproblem.A,2),1);
0183         %equality constraint matrix must be full row rank
0184         [x, f, y, info] = qpng (QPproblem.F, QPproblem.c*QPproblem.osense, full(QPproblem.A), QPproblem.b, ctype, QPproblem.lb, QPproblem.ub, x0);
0185         
0186         f = 0.5*x'*QPproblem.F*x + c'*x;
0187         
0188         w=[];
0189         
0190         if (info.status == 0)
0191             stat = 1;
0192         elseif (info.status == 1)
0193             stat = 0;
0194         else
0195             stat = -1;
0196         end
0197         origStat=info.status;
0198         %%
0199     case 'mosek'
0200         if (~isempty(csense))
0201             b_L(csense == 'E') = b(csense == 'E');
0202             b_U(csense == 'E') = b(csense == 'E');
0203             b_L(csense == 'G') = b(csense == 'G');
0204             b_U(csense == 'G') = inf;
0205             b_L(csense == 'L') = -inf;
0206             b_U(csense == 'L') = b(csense == 'L');
0207         else
0208             b_L = b;
0209             b_U = b;
0210         end
0211                 
0212         if printLevel>0
0213             cmd='minimize';
0214         else
0215             cmd='minimize echo(0)';
0216         end
0217         
0218         % Optimize the problem.
0219         % min 0.5*x'*F*x + osense*c'*x
0220         % st. blc <= A*x <= buc
0221         %     bux <= x   <= bux
0222         [res] = mskqpopt(F,osense*c,A,b_L,b_U,lb,ub,[],cmd);
0223         
0224         if isempty(res)
0225             stat=3;
0226         else
0227             if isfield(res,'sol')
0228                 origStat=res.sol.itr.solsta;
0229                 if strcmp(res.sol.itr.prosta,'PRIMAL_AND_DUAL_FEASIBLE') &&  (strcmp(res.sol.itr.solsta,'OPTIMAL') || strcmp(res.sol.itr.solsta,'NEAR_OPTIMAL'))
0230                     stat=1;
0231                     % x solution.
0232                     x = res.sol.itr.xx;
0233                     f = 0.5*x'*F*x + c'*x;
0234                     
0235                     %dual to equality
0236                     y=res.sol.itr.y;
0237                     
0238                     %dual to lower and upper bounds
0239                     w=res.sol.itr.slx - res.sol.itr.sux;
0240                 else
0241                     stat=3;
0242                 end
0243             else
0244                 stat=3;
0245                 origStat=res.rmsg;
0246             end
0247         end
0248         % stat   Solver status
0249         %           1   Optimal solution found
0250         %           2   Unbounded solution
0251         %           0   Infeasible QP
0252         %           3   Other problem (time limit etc)
0253         %%
0254       case 'pdco'
0255         %-----------------------------------------------------------------------
0256         % pdco.m: Primal-Dual Barrier Method for Convex Objectives (16 Dec 2008)
0257         %-----------------------------------------------------------------------
0258         % AUTHOR:
0259         %    Michael Saunders, Systems Optimization Laboratory (SOL),
0260         %    Stanford University, Stanford, California, USA.
0261         %Interfaced with Cobra toolbox by Ronan Fleming, 18 Jan 2010
0262         [nMet,nRxn]=size(A);
0263         d1=ones(nRxn,1)*1e-4;
0264         %dont minimise the norm of reactions in linear objective
0265         d1(c~=0)=0;
0266         d2=1e-5;
0267         options = pdcoSet;
0268         
0269         x0 = ones(nRxn,1);
0270         y0 = ones(nMet,1);
0271         z0 = ones(nRxn,1);
0272         xsize = 1000;
0273         zsize = 1000;
0274         options.Method=2; %QR
0275         options.MaxIter=100;
0276         options.Print=printLevel;
0277         %get handle to helper function for objective
0278         pdObjHandle = @(x) QPObj(x);
0279         %solve the QP
0280         [x,y,w,inform,PDitns,CGitns,time] = ...
0281             pdco(pdObjHandle,A,b,lb,ub,d1,d2,options,x0,y0,z0,xsize,zsize);
0282         f= c'*x + 0.5*x'*F*x;
0283         % inform = 0 if a solution is found;
0284 %        = 1 if too many iterations were required;
0285 %        = 2 if the linesearch failed too often;
0286 %        = 3 if the step lengths became too small;
0287 %        = 4 if Cholesky said ADDA was not positive definite.
0288         if (inform == 0)
0289             stat = 1;
0290         elseif (inform == 1 || inform == 2 || inform == 3)
0291             stat = 0;
0292         else
0293             stat = -1;
0294         end
0295         origStat=inform;
0296     %%
0297     case 'gurobi'
0298         % Free academic licenses for the Gurobi solver can be obtained from
0299         % http://www.gurobi.com/html/academic.html
0300         %
0301         % The code below uses Gurobi Mex to interface with Gurobi. It can be downloaded from
0302         % http://www.convexoptimization.com/wikimization/index.php/Gurobi_Mex:_A_MATLAB_interface_for_Gurobi
0303 
0304         clear opts            % Use the default parameter settings
0305         if printLevel == 0
0306            % Version v1.10 of Gurobi Mex has a minor bug. For complete silence
0307            % Remove Line 736 of gurobi_mex.c: mexPrintf("\n");
0308            opts.Display = 0;
0309            opts.DisplayInterval = 0;
0310         else
0311            opts.Display = 1;
0312         end
0313 
0314         
0315         
0316         if (isempty(csense))
0317             clear csense
0318             csense(1:length(b),1) = '=';
0319         else
0320             csense(csense == 'L') = '<';
0321             csense(csense == 'G') = '>';
0322             csense(csense == 'E') = '=';
0323             csense = csense(:);
0324         end
0325         
0326         % Gurobi passes individual terms instead of an F matrix. qrow and
0327         % qcol specify which variables are multipled to get each term,
0328         % while qval specifies the coefficients of each term.
0329 
0330         [qrow,qcol,qval]=find(F);
0331         qrow=qrow'-1;   % -1 because gurobi numbers indices from zero, not one.
0332         qcol=qcol'-1;
0333         qval=0.5*qval';
0334         
0335         opts.QP.qrow = int32(qrow); 
0336         opts.QP.qcol = int32(qcol); 
0337         opts.QP.qval = qval;
0338         opts.Method = 0;    % 0 - primal, 1 - dual
0339         opts.Presolve = -1; % -1 - auto, 0 - no, 1 - conserv, 2 - aggressive
0340         opts.FeasibilityTol = 1e-6;
0341         opts.IntFeasTol = 1e-5;
0342         opts.OptimalityTol = 1e-6;
0343         %opt.Quad=1;
0344         
0345         %gurobi_mex doesn't cast logicals to doubles automatically
0346         c = double(c);
0347         [x,f,origStat,output,y] = gurobi_mex(c,osense,sparse(A),b, ...
0348                                              csense,lb,ub,[],opts);
0349         if origStat==2
0350            stat = 1; % Optimal solutuion found
0351         elseif origStat==3
0352            stat = 0; % Infeasible
0353         elseif origStat==5
0354            stat = 2; % Unbounded
0355         elseif origStat==4
0356            stat = 0; % Gurobi reports infeasible *or* unbounded
0357         else
0358            stat = -1; % Solution not optimal or solver problem
0359         end
0360         
0361     case 'gurobi5'
0362      %% gurobi 5
0363      % Free academic licenses for the Gurobi solver can be obtained from
0364         % http://www.gurobi.com/html/academic.html
0365         resultgurobi = struct('x',[],'objval',[],'pi',[]);
0366         clear params            % Use the default parameter settings
0367         if printLevel == 0 
0368             params.OutputFlag = 0;
0369             params.DisplayInterval = 1;
0370         else
0371             params.OutputFlag = 1;
0372             params.DisplayInterval = 5;
0373         end
0374 
0375         params.Method = 0;    %-1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent
0376         params.Presolve = -1; % -1 - auto, 0 - no, 1 - conserv, 2 - aggressive
0377         params.IntFeasTol = 1e-5;
0378         params.FeasibilityTol = 1e-6;
0379         params.OptimalityTol = 1e-6;
0380         %params.Quad = 1;
0381         
0382         if (isempty(QPproblem.csense))
0383             clear QPproblem.csense
0384             QPproblem.csense(1:length(b),1) = '=';
0385         else
0386             QPproblem.csense(QPproblem.csense == 'L') = '<';
0387             QPproblem.csense(QPproblem.csense == 'G') = '>';
0388             QPproblem.csense(QPproblem.csense == 'E') = '=';
0389             QPproblem.csense = QPproblem.csense(:);
0390         end
0391     
0392         if QPproblem.osense == -1
0393             QPproblem.osense = 'max';
0394         else
0395             QPproblem.osense = 'min';
0396         end
0397         
0398         QPproblem.Q = 0.5*sparse(QPproblem.F);
0399         QPproblem.modelsense = QPproblem.osense;
0400         [QPproblem.A,QPproblem.rhs,QPproblem.obj,QPproblem.sense] = deal(sparse(QPproblem.A),QPproblem.b,double(QPproblem.c),QPproblem.csense);
0401         resultgurobi = gurobi(QPproblem,params);
0402         origStat = resultgurobi.status;
0403         if strcmp(resultgurobi.status,'OPTIMAL')
0404            stat = 1; % Optimal solution found
0405            [x,f,y] = deal(resultgurobi.x,resultgurobi.objval,resultgurobi.pi);
0406         elseif strcmp(resultgurobi.status,'INFEASIBLE')
0407            stat = 0; % Infeasible
0408         elseif strcmp(resultgurobi.status,'UNBOUNDED')
0409            stat = 2; % Unbounded
0410         elseif strcmp(resultgurobi.status,'INF_OR_UNBD')
0411            stat = 0; % Gurobi reports infeasible *or* unbounded
0412         else
0413            stat = -1; % Solution not optimal or solver problem
0414         end
0415         
0416     %%
0417     otherwise
0418         error(['Unknown solver: ' solver]);
0419 end
0420 %%
0421 t = etime(clock, t_start);
0422 solution.obj = f;
0423 solution.solver = solver;
0424 solution.stat = stat;
0425 solution.origStat = origStat; 
0426 solution.time = t;
0427 solution.full = x;
0428 solution.dual = y;
0429 solution.rcost = w;
0430 
0431 %Helper function for pdco
0432 %%
0433     function [obj,grad,hess] = QPObj(x)
0434         obj  = c'*x + 0.5*x'*F*x;
0435         grad = c + F*x;
0436         hess = F;
0437     end
0438 end

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003