solveCobraMIQP

PURPOSE ^

solveCobraQP Solve constraint-based QP problems

SYNOPSIS ^

function solution = solveCobraMIQP(MIQPproblem,varargin)

DESCRIPTION ^

solveCobraQP Solve constraint-based QP problems

 solution = solveCobraQP(MIQPproblem,solver,verbFlag,solverParams)

 % Solves problems of the type 

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

INPUT
MIQPproblem    Structure containing the following fields describing the QP
               problem to be solved
  A                LHS matrix
  b                RHS vector
  F                F matrix for quadratic objective (see above)
  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 defined in the CBT_MIQP_SOLVER global variable (set using
 changeCobraSolver). Solvers currently available are 'tomlab_cplex'

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)
                       1   Optimal solution found
                       2   Unbounded solution
                       0   Infeasible QP
                      -1   No optimal solution found (time limit etc)
                       3   Solution exists but with problems
  origStat         Original status returned by the specific solver 
  time             Solve time in seconds


 Markus Herrgard 6/8/07
 Tim Harrington  05/18/12 Added support for the Gurobi 5.0 solver

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function solution = solveCobraMIQP(MIQPproblem,varargin)
0002 %solveCobraQP Solve constraint-based QP problems
0003 %
0004 % solution = solveCobraQP(MIQPproblem,solver,verbFlag,solverParams)
0005 %
0006 % % Solves problems of the type
0007 %
0008 %      min   osense * 0.5 x' * F * x + osense * c' * x
0009 %      s/t   lb <= x <= ub
0010 %            A * x  <=/=/>= b
0011 %            xi = integer
0012 %
0013 %INPUT
0014 %MIQPproblem    Structure containing the following fields describing the QP
0015 %               problem to be solved
0016 %  A                LHS matrix
0017 %  b                RHS vector
0018 %  F                F matrix for quadratic objective (see above)
0019 %  c                Objective coeff vector
0020 %  lb               Lower bound vector
0021 %  ub               Upper bound vector
0022 %  osense           Objective sense (-1 max, +1 min)
0023 %  csense           Constraint senses, a string containting the constraint
0024 %                   sense for each row in A ('E', equality, 'G' greater
0025 %                   than, 'L' less than).
0026 %
0027 %OPTIONAL INPUTS
0028 % Optional parameters can be entered using parameters structure or as
0029 % parameter followed by parameter value: i.e. ,'printLevel',3)
0030 %
0031 % parameters    Structure containing optional parameters as fields.
0032 %  printLevel   Print level for solver
0033 %  saveInput    Saves LPproblem to filename specified in field.
0034 %               Setting parameters = 'default' uses default setting set in
0035 %               getCobraSolverParameters.
0036 %
0037 % the solver defined in the CBT_MIQP_SOLVER global variable (set using
0038 % changeCobraSolver). Solvers currently available are 'tomlab_cplex'
0039 %
0040 %OUTPUT
0041 % solution      Structure containing the following fields describing a QP
0042 %               solution
0043 %  full             Full QP solution vector
0044 %  obj              Objective value
0045 %  solver           Solver used to solve QP problem
0046 %  stat             Solver status in standardized form (see below)
0047 %                       1   Optimal solution found
0048 %                       2   Unbounded solution
0049 %                       0   Infeasible QP
0050 %                      -1   No optimal solution found (time limit etc)
0051 %                       3   Solution exists but with problems
0052 %  origStat         Original status returned by the specific solver
0053 %  time             Solve time in seconds
0054 %
0055 %
0056 % Markus Herrgard 6/8/07
0057 % Tim Harrington  05/18/12 Added support for the Gurobi 5.0 solver
0058 
0059 global CBT_MIQP_SOLVER;
0060 solver = CBT_MIQP_SOLVER;
0061 
0062 %optional parameters
0063 optParamNames = {'printLevel', 'saveInput', 'timeLimit'};
0064 parameters = '';
0065 if nargin ~=1
0066     if mod(length(varargin),2)==0
0067         for i=1:2:length(varargin)-1
0068             if ismember(varargin{i},optParamNames)
0069                 parameters.(varargin{i}) = varargin{i+1};
0070             else
0071                 error([varargin{i} ' is not a valid optional parameter']);
0072             end
0073         end
0074     elseif strcmp(varargin{1},'default')
0075         parameters = 'default';
0076     elseif isstruct(varargin{1})
0077         parameters = varargin{1};
0078     else
0079         display('Warning: Invalid number of parameters/values')
0080         solution=[];
0081         return;
0082     end
0083 end
0084 [printLevel, saveInput, timeLimit] = getCobraSolverParams('QP',optParamNames,parameters);
0085 
0086 [A,b,F,c,lb,ub,csense,osense, vartype] = ...
0087     deal(MIQPproblem.A,MIQPproblem.b,MIQPproblem.F,MIQPproblem.c,MIQPproblem.lb,MIQPproblem.ub,...
0088     MIQPproblem.csense,MIQPproblem.osense, MIQPproblem.vartype);
0089 
0090 t_start = clock;
0091 switch solver
0092 %% CPLEX through TOMLAB
0093     case 'tomlab_cplex'
0094         if (~isempty(csense))
0095             b_L(csense == 'E') = b(csense == 'E');
0096             b_U(csense == 'E') = b(csense == 'E');
0097             b_L(csense == 'G') = b(csense == 'G');
0098             b_U(csense == 'G') = inf;
0099             b_L(csense == 'L') = -inf;
0100             b_U(csense == 'L') = b(csense == 'L');
0101         else
0102             b_L = b;
0103             b_U = b;
0104         end
0105         intVars = find((vartype == 'B') | (vartype == 'I'));
0106         %tomlabProblem = qpAssign(osense*F,osense*c,A,b_L,b_U,lb,ub,[],'CobraQP');
0107         tomlabProblem  = miqpAssign(osense*F, osense*c, A, b_L, b_U, lb, ub,[], ...
0108                              intVars, [],[],[],'CobraMIQP');
0109         tomlabProblem.CPLEX.LogFile = 'MIQPproblem.log';
0110         
0111         %optional parameters
0112         PriLvl = printLevel;
0113         
0114         %Save Input if selected
0115         if ~isempty(saveInput)
0116             fileName = parameters.saveInput;
0117             if ~find(regexp(fileName,'.mat'))
0118                 fileName = [fileName '.mat'];
0119             end
0120             display(['Saving MIQPproblem in ' fileName]);
0121             save(fileName,'MIQPproblem')
0122         end
0123         tomlabProblem.MIP.cpxControl.TILIM = timeLimit; % time limit
0124         tomlabProblem.MIP.cpxControl.THREADS = 1; % by default use only one thread
0125         Result = tomRun('cplex', tomlabProblem, PriLvl);
0126         
0127         x = Result.x_k;
0128         f = osense*Result.f_k;
0129         stat = Result.Inform;
0130         if (stat == 1 ||stat == 101 || stat == 102)
0131             solStat = 1; % Optimal
0132         elseif (stat == 3 || stat == 4)
0133             solStat = 0; % Infeasible
0134         elseif (stat == 103)
0135             solStat = 0; % Integer Infeasible
0136         elseif (stat == 2 || stat == 118 || stat == 119)
0137             solStat = 2; % Unbounded
0138         elseif (stat == 106 || stat == 108 || stat == 110 || stat == 112 || stat == 114 || stat == 117)
0139             solStat = -1; % No integer solution exists
0140         elseif (stat >= 10)
0141             solStat = -1; % No optimal solution found (time or other limits reached, other infeasibility problems)
0142         else
0143             solStat = 3; % Solution exists, but either scaling problems or not proven to be optimal
0144         end
0145             %%
0146     case 'gurobi'
0147         % Free academic licenses for the Gurobi solver can be obtained from
0148         % http://www.gurobi.com/html/academic.html
0149         %
0150         % The code below uses Gurobi Mex to interface with Gurobi. It can be downloaded from
0151         % http://www.convexoptimization.com/wikimization/index.php/Gurobi_Mex:_A_MATLAB_interface_for_Gurobi
0152 
0153         clear opts            % Use the default parameter settings
0154         if printLevel == 0
0155            % Version v1.10 of Gurobi Mex has a minor bug. For complete silence
0156            % Remove Line 736 of gurobi_mex.c: mexPrintf("\n");
0157            opts.Display = 0;
0158            opts.DisplayInterval = 0;
0159         else
0160            opts.Display = 1;
0161         end
0162 
0163         
0164         
0165         if (isempty(csense))
0166             clear csense
0167             csense(1:length(b),1) = '=';
0168         else
0169             csense(csense == 'L') = '<';
0170             csense(csense == 'G') = '>';
0171             csense(csense == 'E') = '=';
0172             csense = csense(:);
0173         end
0174         
0175         % Gurobi passes individual terms instead of an F matrix. qrow and
0176         % qcol specify which variables are multipled to get each term,
0177         % while qval specifies the coefficients of each term.
0178 
0179         [qrow,qcol,qval]=find(F);
0180         qrow=qrow'-1;   % -1 because gurobi numbers indices from zero, not one.
0181         qcol=qcol'-1;
0182         qval=0.5*qval';
0183         
0184         opts.QP.qrow = int32(qrow); 
0185         opts.QP.qcol = int32(qcol); 
0186         opts.QP.qval = qval;
0187         opts.Method = 0;    % 0 - primal, 1 - dual
0188         opts.Presolve = -1; % -1 - auto, 0 - no, 1 - conserv, 2 - aggressive
0189         opts.FeasibilityTol = 1e-6;
0190         opts.IntFeasTol = 1e-5;
0191         opts.OptimalityTol = 1e-6;
0192         %opt.Quad=1;
0193         
0194         %gurobi_mex doesn't cast logicals to doubles automatically
0195         c = double(c);
0196         [x,f,origStat,output,y] = gurobi_mex(c,osense,sparse(A),b, ...
0197                                              csense,lb,ub,vartype,opts);
0198         if origStat==2
0199            stat = 1; % Optimal solutuion found
0200         elseif origStat==3
0201            stat = 0; % Infeasible
0202         elseif origStat==5
0203            stat = 2; % Unbounded
0204         elseif origStat==4
0205            stat = 0; % Gurobi reports infeasible *or* unbounded
0206         else
0207            stat = -1; % Solution not optimal or solver problem
0208         end
0209         solStat = stat;
0210     case 'gurobi5'
0211      %% gurobi5
0212      % Free academic licenses for the Gurobi solver can be obtained from
0213         % http://www.gurobi.com/html/academic.html
0214         resultgurobi = struct('x',[],'objval',[],'pi',[]);
0215         clear params            % Use the default parameter settings
0216         if printLevel == 0 
0217             params.OutputFlag = 0;
0218             params.DisplayInterval = 1;
0219         else
0220             params.OutputFlag = 1;
0221             params.DisplayInterval = 5;
0222         end
0223 
0224         params.Method = 0;    %-1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent
0225         params.Presolve = -1; % -1 - auto, 0 - no, 1 - conserv, 2 - aggressive
0226         params.IntFeasTol = 1e-5;
0227         params.FeasibilityTol = 1e-6;
0228         params.OptimalityTol = 1e-6;
0229         
0230         if (isempty(MIQPproblem.csense))
0231             clear MIQPproblem.csense
0232             MIQPproblem.csense(1:length(b),1) = '=';
0233         else
0234             MIQPproblem.csense(MIQPproblem.csense == 'L') = '<';
0235             MIQPproblem.csense(MIQPproblem.csense == 'G') = '>';
0236             MIQPproblem.csense(MIQPproblem.csense == 'E') = '=';
0237             MIQPproblem.csense = MIQPproblem.csense(:);
0238         end
0239     
0240         if MIQPproblem.osense == -1
0241             MIQPproblem.osense = 'max';
0242         else
0243             MIQPproblem.osense = 'min';
0244         end
0245         
0246         MIQPproblem.vtype = vartype;
0247         MIQPproblem.Q = 0.5*sparse(MIQPproblem.F);
0248         MIQPproblem.modelsense = MIQPproblem.osense;
0249         [MIQPproblem.A,MIQPproblem.rhs,MIQPproblem.obj,MIQPproblem.sense] = deal(sparse(MIQPproblem.A),MIQPproblem.b,MIQPproblem.c,MIQPproblem.csense);
0250         resultgurobi = gurobi(MIQPproblem,params);
0251         solStat = resultgurobi.status;
0252         if strcmp(resultgurobi.status,'OPTIMAL')
0253            stat = 1; % Optimal solution found
0254            
0255            if exist('resultgurobi.pi')
0256                [x,f,y] = deal(resultgurobi.x,resultgurobi.objval,resultgurobi.pi);
0257            else
0258                [x,f] = deal(resultgurobi.x,resultgurobi.objval);
0259            end
0260         elseif strcmp(resultgurobi.status,'INFEASIBLE')
0261            stat = 0; % Infeasible
0262         elseif strcmp(resultgurobi.status,'UNBOUNDED')
0263            stat = 2; % Unbounded
0264         elseif strcmp(resultgurobi.status,'INF_OR_UNBD')
0265            stat = 0; % Gurobi reports infeasible *or* unbounded
0266         else
0267            stat = -1; % Solution not optimal or solver problem
0268         end
0269         %%
0270     otherwise
0271         error(['Unknown solver: ' solver]);
0272 end
0273 %%
0274 t = etime(clock, t_start);
0275 
0276 solution.obj = f;
0277 solution.solver = solver;
0278 solution.stat = solStat;
0279 solution.origStat = stat; 
0280 solution.time = t;
0281 solution.full = x;

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