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

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005