Home > matpower7.1 > mp-opt-model > lib > qps_cplex.m

qps_cplex

PURPOSE ^

QPS_CPLEX Quadratic Program Solver based on CPLEX.

SYNOPSIS ^

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

DESCRIPTION ^

QPS_CPLEX  Quadratic Program Solver based on CPLEX.
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
       QPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, OPT)
   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_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
       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
           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, 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] = ...
           qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt)

       x = qps_cplex(H, c, A, l, u)
       x = qps_cplex(H, c, A, l, u, xmin, xmax)
       x = qps_cplex(H, c, A, l, u, xmin, xmax, x0)
       x = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt)
       x = qps_cplex(problem), where problem is a struct with fields:
                       H, c, A, l, u, xmin, xmax, x0, opt
                       all fields except 'c', 'A' and 'l' or 'u' are optional
       x = qps_cplex(...)
       [x, f] = qps_cplex(...)
       [x, f, exitflag] = qps_cplex(...)
       [x, f, exitflag, output] = qps_cplex(...)
       [x, f, exitflag, output, lambda] = qps_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] = qps_cplex(H, c, A, l, u, xmin, [], x0, opt);

   See also QPS_MASTER, CPLEXQP, CPLEXLP, CPLEX_OPTIONS.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x, f, eflag, output, lambda] = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt)
0002 %QPS_CPLEX  Quadratic Program Solver based on CPLEX.
0003 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
0004 %       QPS_CPLEX(H, C, A, L, U, XMIN, XMAX, X0, OPT)
0005 %   [X, F, EXITFLAG, OUTPUT, LAMBDA] = QPS_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 %       OPT : optional options structure with the following fields,
0028 %           all of which are also optional (default values shown in
0029 %           parentheses)
0030 %           verbose (0) - controls level of progress output displayed
0031 %               0 = no progress output
0032 %               1 = some progress output
0033 %               2 = verbose progress output
0034 %           cplex_opt - options struct for CPLEX, value in verbose
0035 %                   overrides these options
0036 %       PROBLEM : The inputs can alternatively be supplied in a single
0037 %           PROBLEM struct with fields corresponding to the input arguments
0038 %           described above: H, c, A, l, u, xmin, xmax, x0, opt
0039 %
0040 %   Outputs:
0041 %       X : solution vector
0042 %       F : final objective function value
0043 %       EXITFLAG : CPLEXQP/CPLEXLP exit flag
0044 %           (see CPLEXQP and CPLEXLP documentation for details)
0045 %       OUTPUT : CPLEXQP/CPLEXLP output struct
0046 %           (see CPLEXQP and CPLEXLP documentation for details)
0047 %       LAMBDA : struct containing the Langrange and Kuhn-Tucker
0048 %           multipliers on the constraints, with fields:
0049 %           mu_l - lower (left-hand) limit on linear constraints
0050 %           mu_u - upper (right-hand) limit on linear constraints
0051 %           lower - lower bound on optimization variables
0052 %           upper - upper bound on optimization variables
0053 %
0054 %   Note the calling syntax is almost identical to that of QUADPROG
0055 %   from MathWorks' Optimization Toolbox. The main difference is that
0056 %   the linear constraints are specified with A, L, U instead of
0057 %   A, B, Aeq, Beq.
0058 %
0059 %   Calling syntax options:
0060 %       [x, f, exitflag, output, lambda] = ...
0061 %           qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt)
0062 %
0063 %       x = qps_cplex(H, c, A, l, u)
0064 %       x = qps_cplex(H, c, A, l, u, xmin, xmax)
0065 %       x = qps_cplex(H, c, A, l, u, xmin, xmax, x0)
0066 %       x = qps_cplex(H, c, A, l, u, xmin, xmax, x0, opt)
0067 %       x = qps_cplex(problem), where problem is a struct with fields:
0068 %                       H, c, A, l, u, xmin, xmax, x0, opt
0069 %                       all fields except 'c', 'A' and 'l' or 'u' are optional
0070 %       x = qps_cplex(...)
0071 %       [x, f] = qps_cplex(...)
0072 %       [x, f, exitflag] = qps_cplex(...)
0073 %       [x, f, exitflag, output] = qps_cplex(...)
0074 %       [x, f, exitflag, output, lambda] = qps_cplex(...)
0075 %
0076 %
0077 %   Example: (problem from from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm)
0078 %       H = [   1003.1  4.3     6.3     5.9;
0079 %               4.3     2.2     2.1     3.9;
0080 %               6.3     2.1     3.5     4.8;
0081 %               5.9     3.9     4.8     10  ];
0082 %       c = zeros(4,1);
0083 %       A = [   1       1       1       1;
0084 %               0.17    0.11    0.10    0.18    ];
0085 %       l = [1; 0.10];
0086 %       u = [1; Inf];
0087 %       xmin = zeros(4,1);
0088 %       x0 = [1; 0; 0; 1];
0089 %       opt = struct('verbose', 2);
0090 %       [x, f, s, out, lambda] = qps_cplex(H, c, A, l, u, xmin, [], x0, opt);
0091 %
0092 %   See also QPS_MASTER, CPLEXQP, CPLEXLP, CPLEX_OPTIONS.
0093 
0094 %   MP-Opt-Model
0095 %   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
0096 %   by Ray Zimmerman, PSERC Cornell
0097 %
0098 %   This file is part of MP-Opt-Model.
0099 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0100 %   See https://github.com/MATPOWER/mp-opt-model for more info.
0101 
0102 %% check for CPLEX
0103 % if ~have_feature('cplexqp')
0104 %     error('qps_cplex: requires the MATLAB interface for CPLEX');
0105 % end
0106 
0107 %%----- input argument handling  -----
0108 %% gather inputs
0109 if nargin == 1 && isstruct(H)       %% problem struct
0110     p = H;
0111     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0112     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0113     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0114     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0115     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0116     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0117     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0118     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0119     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0120 else                                %% individual args
0121     if nargin < 9
0122         opt = [];
0123         if nargin < 8
0124             x0 = [];
0125             if nargin < 7
0126                 xmax = [];
0127                 if nargin < 6
0128                     xmin = [];
0129                 end
0130             end
0131         end
0132     end
0133 end
0134 
0135 %% define nx, set default values for missing optional inputs
0136 if isempty(H) || ~any(any(H))
0137     if isempty(A) && isempty(xmin) && isempty(xmax)
0138         error('qps_cplex: LP problem must include constraints or variable bounds');
0139     else
0140         if ~isempty(A)
0141             nx = size(A, 2);
0142         elseif ~isempty(xmin)
0143             nx = length(xmin);
0144         else    % if ~isempty(xmax)
0145             nx = length(xmax);
0146         end
0147     end
0148 else
0149     nx = size(H, 1);
0150 end
0151 if isempty(c)
0152     c = zeros(nx, 1);
0153 end
0154 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0155                                  (isempty(u) || all(u == Inf)))
0156     A = sparse(0,nx);           %% no limits => no linear constraints
0157 end
0158 nA = size(A, 1);                %% number of original linear constraints
0159 if isempty(u)                   %% By default, linear inequalities are ...
0160     u = Inf(nA, 1);             %% ... unbounded above and ...
0161 end
0162 if isempty(l)
0163     l = -Inf(nA, 1);            %% ... unbounded below.
0164 end
0165 if isempty(xmin)                %% By default, optimization variables are ...
0166     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0167 end
0168 if isempty(xmax)
0169     xmax = Inf(nx, 1);          %% ... unbounded above.
0170 end
0171 if isempty(x0)
0172     x0 = zeros(nx, 1);
0173 end
0174 
0175 %% default options
0176 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0177     verbose = opt.verbose;
0178 else
0179     verbose = 0;
0180 end
0181 
0182 %% split up linear constraints
0183 ieq = find( abs(u-l) <= eps );          %% equality
0184 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0185 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0186 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0187 Ae = A(ieq, :);
0188 be = u(ieq);
0189 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0190 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0191 
0192 %% grab some dimensions
0193 nlt = length(ilt);      %% number of upper bounded linear inequalities
0194 ngt = length(igt);      %% number of lower bounded linear inequalities
0195 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0196 
0197 %% set up options struct for CPLEX
0198 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt)
0199     cplex_opt = cplex_options(opt.cplex_opt);
0200 else
0201     cplex_opt = cplex_options;
0202 end
0203 
0204 vstr = have_feature('cplex', 'vstr');
0205 vnum = have_feature('cplex', 'vnum');
0206 vrb = max([0 verbose-1]);
0207 if vrb && vnum > 12.002 && vnum < 12.007
0208     cplex_opt.diagnostics   = 'on';
0209 end
0210 if verbose > 2
0211     cplex_opt.display = 'iter';
0212 elseif verbose > 1
0213     cplex_opt.display = 'on';
0214 elseif verbose > 0
0215     cplex_opt.display = 'off';
0216 end
0217 
0218 if isempty(Ai) && isempty(Ae)
0219     unconstrained = 1;
0220     Ae = sparse(1, nx);
0221     be = 0;
0222 else
0223     unconstrained = 0;
0224 end
0225 
0226 %% call the solver
0227 if verbose
0228     alg_names = {
0229         'default',
0230         'primal simplex',
0231         'dual simplex',
0232         'network simplex',
0233         'barrier',
0234         'sifting',
0235         'concurrent'
0236     };
0237 end
0238 if isempty(H) || ~any(any(H))
0239     if verbose
0240         fprintf('CPLEX Version %s -- %s LP solver\n', ...
0241             vstr, alg_names{cplex_opt.lpmethod+1});
0242     end
0243     [x, f, eflag, output, lam] = ...
0244         cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
0245 else
0246     if verbose
0247         fprintf('CPLEX Version %s --  %s QP solver\n', ...
0248             vstr, alg_names{cplex_opt.qpmethod+1});
0249     end
0250     %% ensure H is numerically symmetric
0251     if ~isequal(H, H')
0252         H = (H + H')/2;
0253     end
0254     [x, f, eflag, output, lam] = ...
0255         cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
0256 end
0257 
0258 %% workaround for eflag == 5 (which we have seen return infeasible results)
0259 %%          cplexstatus: 6
0260 %%    cplexstatusstring: 'non-optimal'
0261 %%              message: 'Solution with numerical issues'
0262 if eflag > 1
0263     warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n          cplexstatus: %d\n    cplexstatusstring: ''%s''\n              message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message);
0264     eflag = -100 - eflag;
0265 end
0266 
0267 %% check for empty results (in case optimization failed)
0268 if isempty(x)
0269     x = NaN(nx, 1);
0270 end
0271 if isempty(f)
0272     f = NaN;
0273 end
0274 if isempty(lam)
0275     lam.ineqlin = NaN(length(bi), 1);
0276     lam.eqlin   = NaN(length(be), 1);
0277     lam.lower   = NaN(nx, 1);
0278     lam.upper   = NaN(nx, 1);
0279     mu_l        = NaN(nA, 1);
0280     mu_u        = NaN(nA, 1);
0281 else
0282     mu_l        = zeros(nA, 1);
0283     mu_u        = zeros(nA, 1);
0284 end
0285 if unconstrained
0286     lam.eqlin = [];
0287 end
0288 
0289 %% negate prices depending on version
0290 if vnum < 12.003
0291     lam.eqlin   = -lam.eqlin;
0292     lam.ineqlin = -lam.ineqlin;
0293 end
0294 
0295 %% repackage lambdas
0296 kl = find(lam.eqlin < 0);   %% lower bound binding
0297 ku = find(lam.eqlin > 0);   %% upper bound binding
0298 
0299 mu_l(ieq(kl)) = -lam.eqlin(kl);
0300 mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0301 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0302 
0303 mu_u(ieq(ku)) = lam.eqlin(ku);
0304 mu_u(ilt) = lam.ineqlin(1:nlt);
0305 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0306 
0307 lambda = struct( ...
0308     'mu_l', mu_l, ...
0309     'mu_u', mu_u, ...
0310     'lower', lam.lower, ...
0311     'upper', lam.upper ...
0312 );

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