Home > matpower5.1 > miqps_cplex.m

miqps_cplex

PURPOSE ^

MIQPS_CPLEX Mixed Integer Quadratic Program Solver based on CPLEX.

SYNOPSIS ^

function [x, f, eflag, output, lambda] = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)

DESCRIPTION ^

MIQPS_CPLEX  Mixed Integer Quadratic Program Solver based on CPLEX.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       MIQPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
   A wrapper function providing a MATPOWER standardized interface for using
   CPLEXQP or CPLEXLP to solve the following QP (quadratic programming)
   problem:

       min 1/2 X'*H*X + C'*X
        X

   subject to

       L <= A*X <= U       (linear constraints)
       XMIN <= X <= XMAX   (variable bounds)

   Inputs (all optional except H, C, A and L):
       H : matrix (possibly sparse) of quadratic cost coefficients
       C : vector of linear cost coefficients
       A, L, U : define the optional linear constraints. Default
           values for the elements of L and U are -Inf and Inf,
           respectively.
       XMIN, XMAX : optional lower and upper bounds on the
           X variables, defaults are -Inf and Inf, respectively.
       X0 : optional starting value of optimization vector X
       VTYPE : character string of length NX (number of elements in X),
               or 1 (value applies to all variables in x),
               allowed values are 'C' (continuous), 'B' (binary),
               'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer).
       OPT : optional options structure with the following fields,
           all of which are also optional (default values shown in
           parentheses)
           verbose (0) - controls level of progress output displayed
               0 = no progress output
               1 = some progress output
               2 = verbose progress output
           skip_prices (0) - flag that specifies whether or not to
               skip the price computation stage, in which the problem
               is re-solved for only the continuous variables, with all
               others being constrained to their solved values
           cplex_opt - options struct for CPLEX, value in verbose
               overrides these options
       PROBLEM : The inputs can alternatively be supplied in a single
           PROBLEM struct with fields corresponding to the input arguments
           described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt

   Outputs:
       X : solution vector
       F : final objective function value
       EXITFLAG : CPLEXQP/CPLEXLP exit flag
           (see CPLEXQP and CPLEXLP documentation for details)
       OUTPUT : CPLEXQP/CPLEXLP output struct
           (see CPLEXQP and CPLEXLP documentation for details)
       LAMBDA : struct containing the Langrange and Kuhn-Tucker
           multipliers on the constraints, with fields:
           mu_l - lower (left-hand) limit on linear constraints
           mu_u - upper (right-hand) limit on linear constraints
           lower - lower bound on optimization variables
           upper - upper bound on optimization variables

   Note the calling syntax is almost identical to that of QUADPROG
   from MathWorks' Optimization Toolbox. The main difference is that
   the linear constraints are specified with A, L, U instead of
   A, B, Aeq, Beq.

   Calling syntax options:
       [x, f, exitflag, output, lambda] = ...
           miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)

       x = miqps_cplex(H, c, A, l, u)
       x = miqps_cplex(H, c, A, l, u, xmin, xmax)
       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0)
       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype)
       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
       x = miqps_cplex(problem), where problem is a struct with fields:
                       H, c, A, l, u, xmin, xmax, x0, vtype, opt
                       all fields except 'c', 'A' and 'l' or 'u' are optional
       x = miqps_cplex(...)
       [x, f] = miqps_cplex(...)
       [x, f, exitflag] = miqps_cplex(...)
       [x, f, exitflag, output] = miqps_cplex(...)
       [x, f, exitflag, output, lambda] = miqps_cplex(...)


   Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm)
       H = [   1003.1  4.3     6.3     5.9;
               4.3     2.2     2.1     3.9;
               6.3     2.1     3.5     4.8;
               5.9     3.9     4.8     10  ];
       c = zeros(4,1);
       A = [   1       1       1       1;
               0.17    0.11    0.10    0.18    ];
       l = [1; 0.10];
       u = [1; Inf];
       xmin = zeros(4,1);
       x0 = [1; 0; 0; 1];
       opt = struct('verbose', 2);
       [x, f, s, out, lambda] = miqps_cplex(H, c, A, l, u, xmin, [], x0, vtype, opt);

   See also CPLEXMIQP, CPLEXMILP, CPLEXQP, CPLEXLP, CPLEX_OPTIONS.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0002 %MIQPS_CPLEX  Mixed Integer Quadratic Program Solver based on CPLEX.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       MIQPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, VTYPE, OPT)
0005 %   A wrapper function providing a MATPOWER standardized interface for using
0006 %   CPLEXQP or CPLEXLP to solve the following QP (quadratic programming)
0007 %   problem:
0008 %
0009 %       min 1/2 X'*H*X + C'*X
0010 %        X
0011 %
0012 %   subject to
0013 %
0014 %       L <= A*X <= U       (linear constraints)
0015 %       XMIN <= X <= XMAX   (variable bounds)
0016 %
0017 %   Inputs (all optional except H, C, A and L):
0018 %       H : matrix (possibly sparse) of quadratic cost coefficients
0019 %       C : vector of linear cost coefficients
0020 %       A, L, U : define the optional linear constraints. Default
0021 %           values for the elements of L and U are -Inf and Inf,
0022 %           respectively.
0023 %       XMIN, XMAX : optional lower and upper bounds on the
0024 %           X variables, defaults are -Inf and Inf, respectively.
0025 %       X0 : optional starting value of optimization vector X
0026 %       VTYPE : character string of length NX (number of elements in X),
0027 %               or 1 (value applies to all variables in x),
0028 %               allowed values are 'C' (continuous), 'B' (binary),
0029 %               'I' (integer), 'S' (semi-continuous), or 'N' (semi-integer).
0030 %       OPT : optional options structure with the following fields,
0031 %           all of which are also optional (default values shown in
0032 %           parentheses)
0033 %           verbose (0) - controls level of progress output displayed
0034 %               0 = no progress output
0035 %               1 = some progress output
0036 %               2 = verbose progress output
0037 %           skip_prices (0) - flag that specifies whether or not to
0038 %               skip the price computation stage, in which the problem
0039 %               is re-solved for only the continuous variables, with all
0040 %               others being constrained to their solved values
0041 %           cplex_opt - options struct for CPLEX, value in verbose
0042 %               overrides these options
0043 %       PROBLEM : The inputs can alternatively be supplied in a single
0044 %           PROBLEM struct with fields corresponding to the input arguments
0045 %           described above: H, c, A, l, u, xmin, xmax, x0, vtype, opt
0046 %
0047 %   Outputs:
0048 %       X : solution vector
0049 %       F : final objective function value
0050 %       EXITFLAG : CPLEXQP/CPLEXLP exit flag
0051 %           (see CPLEXQP and CPLEXLP documentation for details)
0052 %       OUTPUT : CPLEXQP/CPLEXLP output struct
0053 %           (see CPLEXQP and CPLEXLP documentation for details)
0054 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0055 %           multipliers on the constraints, with fields:
0056 %           mu_l - lower (left-hand) limit on linear constraints
0057 %           mu_u - upper (right-hand) limit on linear constraints
0058 %           lower - lower bound on optimization variables
0059 %           upper - upper bound on optimization variables
0060 %
0061 %   Note the calling syntax is almost identical to that of QUADPROG
0062 %   from MathWorks' Optimization Toolbox. The main difference is that
0063 %   the linear constraints are specified with A, L, U instead of
0064 %   A, B, Aeq, Beq.
0065 %
0066 %   Calling syntax options:
0067 %       [x, f, exitflag, output, lambda] = ...
0068 %           miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0069 %
0070 %       x = miqps_cplex(H, c, A, l, u)
0071 %       x = miqps_cplex(H, c, A, l, u, xmin, xmax)
0072 %       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0)
0073 %       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype)
0074 %       x = miqps_cplex(H, c, A, l, u, xmin, xmax, x0, vtype, opt)
0075 %       x = miqps_cplex(problem), where problem is a struct with fields:
0076 %                       H, c, A, l, u, xmin, xmax, x0, vtype, opt
0077 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0078 %       x = miqps_cplex(...)
0079 %       [x, f] = miqps_cplex(...)
0080 %       [x, f, exitflag] = miqps_cplex(...)
0081 %       [x, f, exitflag, output] = miqps_cplex(...)
0082 %       [x, f, exitflag, output, lambda] = miqps_cplex(...)
0083 %
0084 %
0085 %   Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm)
0086 %       H = [   1003.1  4.3     6.3     5.9;
0087 %               4.3     2.2     2.1     3.9;
0088 %               6.3     2.1     3.5     4.8;
0089 %               5.9     3.9     4.8     10  ];
0090 %       c = zeros(4,1);
0091 %       A = [   1       1       1       1;
0092 %               0.17    0.11    0.10    0.18    ];
0093 %       l = [1; 0.10];
0094 %       u = [1; Inf];
0095 %       xmin = zeros(4,1);
0096 %       x0 = [1; 0; 0; 1];
0097 %       opt = struct('verbose', 2);
0098 %       [x, f, s, out, lambda] = miqps_cplex(H, c, A, l, u, xmin, [], x0, vtype, opt);
0099 %
0100 %   See also CPLEXMIQP, CPLEXMILP, CPLEXQP, CPLEXLP, CPLEX_OPTIONS.
0101 
0102 %   MATPOWER
0103 %   Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC)
0104 %   by Ray Zimmerman, PSERC Cornell
0105 %
0106 %   $Id: miqps_cplex.m 2661 2015-03-20 17:02:46Z ray $
0107 %
0108 %   This file is part of MATPOWER.
0109 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0110 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0111 
0112 %% check for CPLEX
0113 % if ~have_fcn('cplexqp')
0114 %     error('miqps_cplex: requires the Matlab interface for CPLEX');
0115 % end
0116 
0117 %%----- input argument handling  -----
0118 %% gather inputs
0119 if nargin == 1 && isstruct(H)       %% problem struct
0120     p = H;
0121     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0122     if isfield(p, 'vtype'), vtype = p.vtype;else,   vtype = []; end
0123     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0124     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0125     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0126     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0127     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0128     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0129     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0130     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0131 else                                %% individual args
0132     if nargin < 10
0133         opt = [];
0134         if nargin < 9
0135             vtype = [];
0136             if nargin < 8
0137                 x0 = [];
0138                 if nargin < 7
0139                     xmax = [];
0140                     if nargin < 6
0141                         xmin = [];
0142                     end
0143                 end
0144             end
0145         end
0146     end
0147 end
0148 
0149 %% define nx, set default values for missing optional inputs
0150 if isempty(H) || ~any(any(H))
0151     if isempty(A) && isempty(xmin) && isempty(xmax)
0152         error('miqps_cplex: LP problem must include constraints or variable bounds');
0153     else
0154         if ~isempty(A)
0155             nx = size(A, 2);
0156         elseif ~isempty(xmin)
0157             nx = length(xmin);
0158         else    % if ~isempty(xmax)
0159             nx = length(xmax);
0160         end
0161     end
0162 else
0163     nx = size(H, 1);
0164 end
0165 if isempty(c)
0166     c = zeros(nx, 1);
0167 end
0168 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0169                                  (isempty(u) || all(u == Inf)))
0170     A = sparse(0,nx);           %% no limits => no linear constraints
0171 end
0172 nA = size(A, 1);                %% number of original linear constraints
0173 if isempty(u)                   %% By default, linear inequalities are ...
0174     u = Inf(nA, 1);             %% ... unbounded above and ...
0175 end
0176 if isempty(l)
0177     l = -Inf(nA, 1);            %% ... unbounded below.
0178 end
0179 if isempty(xmin)                %% By default, optimization variables are ...
0180     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0181 end
0182 if isempty(xmax)
0183     xmax = Inf(nx, 1);          %% ... unbounded above.
0184 end
0185 if isempty(x0)
0186     x0 = zeros(nx, 1);
0187 end
0188 
0189 %% default options
0190 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0191     verbose = opt.verbose;
0192 else
0193     verbose = 0;
0194 end
0195 
0196 %% split up linear constraints
0197 ieq = find( abs(u-l) <= eps );          %% equality
0198 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0199 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0200 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0201 Ae = A(ieq, :);
0202 be = u(ieq);
0203 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0204 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0205 
0206 %% grab some dimensions
0207 nlt = length(ilt);      %% number of upper bounded linear inequalities
0208 ngt = length(igt);      %% number of lower bounded linear inequalities
0209 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0210 
0211 %% set up options struct for CPLEX
0212 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt)
0213     cplex_opt = cplex_options(opt.cplex_opt);
0214 else
0215     cplex_opt = cplex_options;
0216 end
0217 
0218 vstr = have_fcn('cplex', 'vstr');
0219 vnum = have_fcn('cplex', 'vnum');
0220 vrb = max([0 verbose-1]);
0221 cplex_opt.barrier.display   = vrb;
0222 cplex_opt.conflict.display  = vrb;
0223 cplex_opt.mip.display       = vrb;
0224 cplex_opt.sifting.display   = vrb;
0225 cplex_opt.simplex.display   = vrb;
0226 cplex_opt.tune.display      = vrb;
0227 if vrb && vnum > 12.002
0228     cplex_opt.diagnostics   = 'on';
0229 end
0230 if verbose > 2
0231     cplex_opt.Display = 'iter';
0232 elseif verbose > 1
0233     cplex_opt.Display = 'on';
0234 elseif verbose > 0
0235     cplex_opt.Display = 'off';
0236 end
0237 
0238 if isempty(Ai) && isempty(Ae)
0239     unconstrained = 1;
0240     Ae = sparse(1, nx);
0241     be = 0;
0242 else
0243     unconstrained = 0;
0244 end
0245 
0246 %% call the solver
0247 if verbose
0248     methods = {
0249         'default',
0250         'primal simplex',
0251         'dual simplex',
0252         'network simplex',
0253         'barrier',
0254         'sifting',
0255         'concurrent'
0256     };
0257 end
0258 if isempty(vtype) || isempty(find(vtype == 'B' | vtype == 'I' | ...
0259         vtype == 'S' | vtype == 'N'))
0260     mi = 0;
0261 else
0262     mi = 1;
0263     %% expand vtype to nx elements if necessary
0264     if length(vtype) == 1 && nx > 1
0265         vtype = char(vtype * ones(1,nx));
0266     end
0267 end
0268 
0269 if mi
0270     if isempty(H) || ~any(any(H))
0271         if verbose
0272             fprintf('CPLEX Version %s -- %s MILP solver\n', ...
0273                 vstr, methods{cplex_opt.lpmethod+1});
0274         end
0275         [x, f, eflag, output] = ...
0276             cplexmilp(c, Ai, bi, Ae, be, [], [], [], xmin, xmax, vtype, x0, cplex_opt);
0277         lam = [];
0278     else
0279         if verbose
0280             fprintf('CPLEX Version %s --  %s MIQP solver\n', ...
0281                 vstr, methods{cplex_opt.qpmethod+1});
0282         end
0283         %% ensure H is numerically symmetric
0284         if ~isequal(H, H')
0285             H = (H + H')/2;
0286         end
0287         [x, f, eflag, output] = ...
0288             cplexmiqp(H, c, Ai, bi, Ae, be, [], [], [], xmin, xmax, vtype, x0, cplex_opt);
0289     end
0290     lam = [];
0291 else
0292     if isempty(H) || ~any(any(H))
0293         if verbose
0294             fprintf('CPLEX Version %s -- %s LP solver\n', ...
0295                 vstr, methods{cplex_opt.lpmethod+1});
0296         end
0297         [x, f, eflag, output, lam] = ...
0298             cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
0299     else
0300         if verbose
0301             fprintf('CPLEX Version %s --  %s QP solver\n', ...
0302                 vstr, methods{cplex_opt.qpmethod+1});
0303         end
0304         %% ensure H is numerically symmetric
0305         if ~isequal(H, H')
0306             H = (H + H')/2;
0307         end
0308         [x, f, eflag, output, lam] = ...
0309             cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
0310     end
0311 end
0312 
0313 %% workaround for eflag == 5 (which we have seen return infeasible results)
0314 %%          cplexstatus: 6
0315 %%    cplexstatusstring: 'non-optimal'
0316 %%              message: 'Solution with numerical issues'
0317 if eflag > 1
0318     warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n          cplexstatus: %d\n    cplexstatusstring: ''%s''\n              message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message);
0319     if eflag == 5 && mi
0320         eflag = 1;      %% give it a try for the MI phase
0321     else
0322         eflag = -100 - eflag;
0323     end
0324 end
0325 
0326 %% check for empty results (in case optimization failed)
0327 if isempty(x)
0328     x = NaN(nx, 1);
0329 end
0330 if isempty(f)
0331     f = NaN;
0332 end
0333 if isempty(lam)
0334     lam.ineqlin = NaN(length(bi), 1);
0335     lam.eqlin   = NaN(length(be), 1);
0336     lam.lower   = NaN(nx, 1);
0337     lam.upper   = NaN(nx, 1);
0338     mu_l        = NaN(nA, 1);
0339     mu_u        = NaN(nA, 1);
0340 else
0341     mu_l        = zeros(nA, 1);
0342     mu_u        = zeros(nA, 1);
0343 end
0344 if unconstrained
0345     lam.eqlin = [];
0346 end
0347 
0348 %% negate prices depending on version
0349 if vnum < 12.003
0350     lam.eqlin   = -lam.eqlin;
0351     lam.ineqlin = -lam.ineqlin;
0352 end
0353 
0354 %% repackage lambdas
0355 kl = find(lam.eqlin < 0);   %% lower bound binding
0356 ku = find(lam.eqlin > 0);   %% upper bound binding
0357 
0358 mu_l(ieq(kl)) = -lam.eqlin(kl);
0359 mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0360 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0361 
0362 mu_u(ieq(ku)) = lam.eqlin(ku);
0363 mu_u(ilt) = lam.ineqlin(1:nlt);
0364 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0365 
0366 lambda = struct( ...
0367     'mu_l', mu_l, ...
0368     'mu_u', mu_u, ...
0369     'lower', lam.lower, ...
0370     'upper', lam.upper ...
0371 );
0372 
0373 if mi && eflag == 1 && (~isfield(opt, 'skip_prices') || ~opt.skip_prices)
0374     if verbose
0375         fprintf('--- Integer stage complete, starting price computation stage ---\n');
0376     end
0377     tol = 1e-7;
0378     k = find(vtype == 'I' | vtype == 'B' | vtype == 'N' | ...
0379             (vtype == 'S' & x' == 0));
0380     x(k) = round(x(k));
0381     xmin(k) = x(k);
0382     xmax(k) = x(k);
0383     x0 = x;
0384     opt.cplex_opt.lpmethod = 1;     %% primal simplex
0385     opt.cplex_opt.qpmethod = 1;     %% primal simplex
0386     
0387     [x_, f_, eflag_, output_, lambda] = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt);
0388     if eflag ~= eflag_
0389         error('miqps_cplex: EXITFLAG from price computation stage = %d', eflag_);
0390     end
0391     if abs(f - f_)/max(abs(f), 1) > tol
0392         warning('miqps_cplex: relative mismatch in objective function value from price computation stage = %g', abs(f - f_)/max(abs(f), 1));
0393     end
0394     xn = x;
0395     xn(abs(xn)<1) = 1;
0396     [mx, k] = max(abs(x - x_) ./ xn);
0397     if mx > tol
0398         warning('miqps_cplex: max relative mismatch in x from price computation stage = %g (%g)', mx, x(k));
0399     end
0400     output.price_stage = output_;
0401 end

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005