Home > matpower5.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 %   $Id: qps_cplex.m 2340 2014-06-27 19:15:10Z ray $
0095 %   by Ray Zimmerman, PSERC Cornell
0096 %   Copyright (c) 2010-2013 by Power System Engineering Research Center (PSERC)
0097 %
0098 %   This file is part of MATPOWER.
0099 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0100 %
0101 %   MATPOWER is free software: you can redistribute it and/or modify
0102 %   it under the terms of the GNU General Public License as published
0103 %   by the Free Software Foundation, either version 3 of the License,
0104 %   or (at your option) any later version.
0105 %
0106 %   MATPOWER is distributed in the hope that it will be useful,
0107 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0108 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0109 %   GNU General Public License for more details.
0110 %
0111 %   You should have received a copy of the GNU General Public License
0112 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0113 %
0114 %   Additional permission under GNU GPL version 3 section 7
0115 %
0116 %   If you modify MATPOWER, or any covered work, to interface with
0117 %   other modules (such as MATLAB code and MEX-files) available in a
0118 %   MATLAB(R) or comparable environment containing parts covered
0119 %   under other licensing terms, the licensors of MATPOWER grant
0120 %   you additional permission to convey the resulting work.
0121 
0122 %% check for CPLEX
0123 % if ~have_fcn('cplexqp')
0124 %     error('qps_cplex: requires the Matlab interface for CPLEX');
0125 % end
0126 
0127 %%----- input argument handling  -----
0128 %% gather inputs
0129 if nargin == 1 && isstruct(H)       %% problem struct
0130     p = H;
0131     if isfield(p, 'opt'),   opt = p.opt;    else,   opt = [];   end
0132     if isfield(p, 'x0'),    x0 = p.x0;      else,   x0 = [];    end
0133     if isfield(p, 'xmax'),  xmax = p.xmax;  else,   xmax = [];  end
0134     if isfield(p, 'xmin'),  xmin = p.xmin;  else,   xmin = [];  end
0135     if isfield(p, 'u'),     u = p.u;        else,   u = [];     end
0136     if isfield(p, 'l'),     l = p.l;        else,   l = [];     end
0137     if isfield(p, 'A'),     A = p.A;        else,   A = [];     end
0138     if isfield(p, 'c'),     c = p.c;        else,   c = [];     end
0139     if isfield(p, 'H'),     H = p.H;        else,   H = [];     end
0140 else                                %% individual args
0141     if nargin < 9
0142         opt = [];
0143         if nargin < 8
0144             x0 = [];
0145             if nargin < 7
0146                 xmax = [];
0147                 if nargin < 6
0148                     xmin = [];
0149                 end
0150             end
0151         end
0152     end
0153 end
0154 
0155 %% define nx, set default values for missing optional inputs
0156 if isempty(H) || ~any(any(H))
0157     if isempty(A) && isempty(xmin) && isempty(xmax)
0158         error('qps_cplex: LP problem must include constraints or variable bounds');
0159     else
0160         if ~isempty(A)
0161             nx = size(A, 2);
0162         elseif ~isempty(xmin)
0163             nx = length(xmin);
0164         else    % if ~isempty(xmax)
0165             nx = length(xmax);
0166         end
0167     end
0168 else
0169     nx = size(H, 1);
0170 end
0171 if isempty(c)
0172     c = zeros(nx, 1);
0173 end
0174 if  ~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
0175                    (isempty(u) || all(u == Inf))
0176     A = sparse(0,nx);           %% no limits => no linear constraints
0177 end
0178 nA = size(A, 1);                %% number of original linear constraints
0179 if isempty(u)                   %% By default, linear inequalities are ...
0180     u = Inf * ones(nA, 1);      %% ... unbounded above and ...
0181 end
0182 if isempty(l)
0183     l = -Inf * ones(nA, 1);     %% ... unbounded below.
0184 end
0185 if isempty(xmin)                %% By default, optimization variables are ...
0186     xmin = -Inf * ones(nx, 1);  %% ... unbounded below and ...
0187 end
0188 if isempty(xmax)
0189     xmax = Inf * ones(nx, 1);   %% ... unbounded above.
0190 end
0191 if isempty(x0)
0192     x0 = zeros(nx, 1);
0193 end
0194 
0195 %% default options
0196 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
0197     verbose = opt.verbose;
0198 else
0199     verbose = 0;
0200 end
0201 
0202 %% split up linear constraints
0203 ieq = find( abs(u-l) <= eps );          %% equality
0204 igt = find( u >=  1e10 & l > -1e10 );   %% greater than, unbounded above
0205 ilt = find( l <= -1e10 & u <  1e10 );   %% less than, unbounded below
0206 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
0207 Ae = A(ieq, :);
0208 be = u(ieq);
0209 Ai  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
0210 bi  = [ u(ilt);    -l(igt);    u(ibx);    -l(ibx)];
0211 
0212 %% grab some dimensions
0213 nlt = length(ilt);      %% number of upper bounded linear inequalities
0214 ngt = length(igt);      %% number of lower bounded linear inequalities
0215 nbx = length(ibx);      %% number of doubly bounded linear inequalities
0216 
0217 %% set up options struct for CPLEX
0218 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt)
0219     cplex_opt = cplex_options(opt.cplex_opt);
0220 else
0221     cplex_opt = cplex_options;
0222 end
0223 
0224 cplex = Cplex('null');
0225 vstr = cplex.getVersion;
0226 [s,e,tE,m,t] = regexp(vstr, '(\d+\.\d+)\.');
0227 vnum = str2num(t{1}{1});
0228 vrb = max([0 verbose-1]);
0229 cplex_opt.barrier.display   = vrb;
0230 cplex_opt.conflict.display  = vrb;
0231 cplex_opt.mip.display       = vrb;
0232 cplex_opt.sifting.display   = vrb;
0233 cplex_opt.simplex.display   = vrb;
0234 cplex_opt.tune.display      = vrb;
0235 if vrb && vnum > 12.2
0236     cplex_opt.diagnostics   = 'on';
0237 end
0238 if verbose > 2
0239     cplex_opt.Display = 'iter';
0240 elseif verbose > 1
0241     cplex_opt.Display = 'on';
0242 elseif verbose > 0
0243     cplex_opt.Display = 'off';
0244 end
0245 
0246 if isempty(Ai) && isempty(Ae)
0247     unconstrained = 1;
0248     Ae = sparse(1, nx);
0249     be = 0;
0250 else
0251     unconstrained = 0;
0252 end
0253 
0254 %% call the solver
0255 if verbose
0256     methods = {
0257         'default',
0258         'primal simplex',
0259         'dual simplex',
0260         'network simplex',
0261         'barrier',
0262         'sifting',
0263         'concurrent'
0264     };
0265 end
0266 if isempty(H) || ~any(any(H))
0267     if verbose
0268         fprintf('CPLEX Version %s -- %s LP solver\n', ...
0269             vstr, methods{cplex_opt.lpmethod+1});
0270     end
0271     [x, f, eflag, output, lam] = ...
0272         cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
0273 else
0274     if verbose
0275         fprintf('CPLEX Version %s --  %s QP solver\n', ...
0276             vstr, methods{cplex_opt.qpmethod+1});
0277     end
0278     %% ensure H is numerically symmetric
0279     if ~isequal(H, H')
0280         H = (H + H')/2;
0281     end
0282     [x, f, eflag, output, lam] = ...
0283         cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt);
0284 end
0285 
0286 %% workaround for eflag == 5 (which we have seen return infeasible results)
0287 %%          cplexstatus: 6
0288 %%    cplexstatusstring: 'non-optimal'
0289 %%              message: 'Solution with numerical issues'
0290 if eflag > 1
0291     warning('qps_cplex: Undocumented ''exitflag'' value (%d)\n          cplexstatus: %d\n    cplexstatusstring: ''%s''\n              message: ''%s''', eflag, output.cplexstatus, output.cplexstatusstring, output.message);
0292     eflag = -100 - eflag;
0293 end
0294 
0295 %% check for empty results (in case optimization failed)
0296 if isempty(x)
0297     x = NaN(nx, 1);
0298 end
0299 if isempty(f)
0300     f = NaN;
0301 end
0302 if isempty(lam)
0303     lam.ineqlin = NaN(length(bi), 1);
0304     lam.eqlin   = NaN(length(be), 1);
0305     lam.lower   = NaN(nx, 1);
0306     lam.upper   = NaN(nx, 1);
0307     mu_l        = NaN(nA, 1);
0308     mu_u        = NaN(nA, 1);
0309 else
0310     mu_l        = zeros(nA, 1);
0311     mu_u        = zeros(nA, 1);
0312 end
0313 if unconstrained
0314     lam.eqlin = [];
0315 end
0316 
0317 %% negate prices depending on version
0318 if vnum < 12.3
0319     lam.eqlin   = -lam.eqlin;
0320     lam.ineqlin = -lam.ineqlin;
0321 end
0322 
0323 %% repackage lambdas
0324 kl = find(lam.eqlin < 0);   %% lower bound binding
0325 ku = find(lam.eqlin > 0);   %% upper bound binding
0326 
0327 mu_l(ieq(kl)) = -lam.eqlin(kl);
0328 mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
0329 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
0330 
0331 mu_u(ieq(ku)) = lam.eqlin(ku);
0332 mu_u(ilt) = lam.ineqlin(1:nlt);
0333 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
0334 
0335 lambda = struct( ...
0336     'mu_l', mu_l, ...
0337     'mu_u', mu_u, ...
0338     'lower', lam.lower, ...
0339     'upper', lam.upper ...
0340 );

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005