Home > matpower6.0 > 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)
   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
       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 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] = qps_cplex(H, c, A, l, u, xmin, [], x0, opt);

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

Generated on Fri 16-Dec-2016 12:45:37 by m2html © 2005