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

Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005