Home > matpower5.1 > 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-2015 by Power System Engineering Research Center (PSERC)
0095 %   by Ray Zimmerman, PSERC Cornell
0096 %
0097 %   $Id: qps_cplex.m 2661 2015-03-20 17:02:46Z ray $
0098 %
0099 %   This file is part of MATPOWER.
0100 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0101 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0102 
0103 %% check for CPLEX
0104 % if ~have_fcn('cplexqp')
0105 %     error('qps_cplex: requires the Matlab interface for CPLEX');
0106 % end
0107 
0108 %%----- input argument handling  -----
0109 %% gather inputs
0110 if nargin == 1 && isstruct(H)       %% problem struct
0111     p = H;
0112     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0113     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0114     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0115     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0116     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0117     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0118     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0119     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0120     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0121 else                                %% individual args
0122     if nargin < 9
0123         opt = [];
0124         if nargin < 8
0125             x0 = [];
0126             if nargin < 7
0127                 xmax = [];
0128                 if nargin < 6
0129                     xmin = [];
0130                 end
0131             end
0132         end
0133     end
0134 end
0135 
0136 %% define nx, set default values for missing optional inputs
0137 if isempty(H) || ~any(any(H))
0138     if isempty(A) && isempty(xmin) && isempty(xmax)
0139         error('qps_cplex: LP problem must include constraints or variable bounds');
0140     else
0141         if ~isempty(A)
0142             nx = size(A, 2);
0143         elseif ~isempty(xmin)
0144             nx = length(xmin);
0145         else    % if ~isempty(xmax)
0146             nx = length(xmax);
0147         end
0148     end
0149 else
0150     nx = size(H, 1);
0151 end
0152 if isempty(c)
0153     c = zeros(nx, 1);
0154 end
0155 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0156                                  (isempty(u) || all(u == Inf)))
0157     A = sparse(0,nx);           %% no limits => no linear constraints
0158 end
0159 nA = size(A, 1);                %% number of original linear constraints
0160 if isempty(u)                   %% By default, linear inequalities are ...
0161     u = Inf(nA, 1);             %% ... unbounded above and ...
0162 end
0163 if isempty(l)
0164     l = -Inf(nA, 1);            %% ... unbounded below.
0165 end
0166 if isempty(xmin)                %% By default, optimization variables are ...
0167     xmin = -Inf(nx, 1);         %% ... unbounded below and ...
0168 end
0169 if isempty(xmax)
0170     xmax = Inf(nx, 1);          %% ... unbounded above.
0171 end
0172 if isempty(x0)
0173     x0 = zeros(nx, 1);
0174 end
0175 
0176 %% default options
0177 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0178     verbose = opt.verbose;
0179 else
0180     verbose = 0;
0181 end
0182 
0183 %% split up linear constraints
0184 ieq = find( abs(u-l) <= eps );          %% equality
0185 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0186 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0187 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0188 Ae = A(ieq, :);
0189 be = u(ieq);
0190 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0191 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0192 
0193 %% grab some dimensions
0194 nlt = length(ilt);      %% number of upper bounded linear inequalities
0195 ngt = length(igt);      %% number of lower bounded linear inequalities
0196 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0197 
0198 %% set up options struct for CPLEX
0199 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt)
0200     cplex_opt = cplex_options(opt.cplex_opt);
0201 else
0202     cplex_opt = cplex_options;
0203 end
0204 
0205 vstr = have_fcn('cplex', 'vstr');
0206 vnum = have_fcn('cplex', 'vnum');
0207 vrb = max([0 verbose-1]);
0208 cplex_opt.barrier.display   = vrb;
0209 cplex_opt.conflict.display  = vrb;
0210 cplex_opt.mip.display       = vrb;
0211 cplex_opt.sifting.display   = vrb;
0212 cplex_opt.simplex.display   = vrb;
0213 cplex_opt.tune.display      = vrb;
0214 if vrb && vnum > 12.002
0215     cplex_opt.diagnostics   = 'on';
0216 end
0217 if verbose > 2
0218     cplex_opt.Display = 'iter';
0219 elseif verbose > 1
0220     cplex_opt.Display = 'on';
0221 elseif verbose > 0
0222     cplex_opt.Display = 'off';
0223 end
0224 
0225 if isempty(Ai) && isempty(Ae)
0226     unconstrained = 1;
0227     Ae = sparse(1, nx);
0228     be = 0;
0229 else
0230     unconstrained = 0;
0231 end
0232 
0233 %% call the solver
0234 if verbose
0235     methods = {
0236         'default',
0237         'primal simplex',
0238         'dual simplex',
0239         'network simplex',
0240         'barrier',
0241         'sifting',
0242         'concurrent'
0243     };
0244 end
0245 if isempty(H) || ~any(any(H))
0246     if verbose
0247         fprintf('CPLEX Version %s -- %s LP solver\n', ...
0248             vstr, methods{cplex_opt.lpmethod+1});
0249     end
0250     [x, f, eflag, output, lam] = ...
0251         cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
0252 else
0253     if verbose
0254         fprintf('CPLEX Version %s --  %s QP solver\n', ...
0255             vstr, methods{cplex_opt.qpmethod+1});
0256     end
0257     %% ensure H is numerically symmetric
0258     if ~isequal(H, H')
0259         H = (H + H')/2;
0260     end
0261     [x, f, eflag, output, lam] = ...
0262         cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
0263 end
0264 
0265 %% workaround for eflag == 5 (which we have seen return infeasible results)
0266 %%          cplexstatus: 6
0267 %%    cplexstatusstring: 'non-optimal'
0268 %%              message: 'Solution with numerical issues'
0269 if eflag > 1
0270     warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n          cplexstatus: %d\n    cplexstatusstring: ''%s''\n              message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message);
0271     eflag = -100 - eflag;
0272 end
0273 
0274 %% check for empty results (in case optimization failed)
0275 if isempty(x)
0276     x = NaN(nx, 1);
0277 end
0278 if isempty(f)
0279     f = NaN;
0280 end
0281 if isempty(lam)
0282     lam.ineqlin = NaN(length(bi), 1);
0283     lam.eqlin   = NaN(length(be), 1);
0284     lam.lower   = NaN(nx, 1);
0285     lam.upper   = NaN(nx, 1);
0286     mu_l        = NaN(nA, 1);
0287     mu_u        = NaN(nA, 1);
0288 else
0289     mu_l        = zeros(nA, 1);
0290     mu_u        = zeros(nA, 1);
0291 end
0292 if unconstrained
0293     lam.eqlin = [];
0294 end
0295 
0296 %% negate prices depending on version
0297 if vnum < 12.003
0298     lam.eqlin   = -lam.eqlin;
0299     lam.ineqlin = -lam.ineqlin;
0300 end
0301 
0302 %% repackage lambdas
0303 kl = find(lam.eqlin < 0);   %% lower bound binding
0304 ku = find(lam.eqlin > 0);   %% upper bound binding
0305 
0306 mu_l(ieq(kl)) = -lam.eqlin(kl);
0307 mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0308 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0309 
0310 mu_u(ieq(ku)) = lam.eqlin(ku);
0311 mu_u(ilt) = lam.ineqlin(1:nlt);
0312 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0313 
0314 lambda = struct( ...
0315     'mu_l', mu_l, ...
0316     'mu_u', mu_u, ...
0317     'lower', lam.lower, ...
0318     'upper', lam.upper ...
0319 );

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