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.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 0034 % verbose 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,v 1.8 2011/12/12 16:07:39 cvs Exp $ 0095 % by Ray Zimmerman, PSERC Cornell 0096 % Copyright (c) 2010, 2011 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 % if ~isempty(opt) && isfield(opt, 'max_it') && ~isempty(opt.max_it) 0202 % max_it = opt.max_it; 0203 % else 0204 % max_it = 0; 0205 % end 0206 0207 %% split up linear constraints 0208 ieq = find( abs(u-l) <= eps ); %% equality 0209 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0210 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0211 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0212 Ae = A(ieq, :); 0213 be = u(ieq); 0214 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0215 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0216 0217 %% grab some dimensions 0218 nlt = length(ilt); %% number of upper bounded linear inequalities 0219 ngt = length(igt); %% number of lower bounded linear inequalities 0220 nbx = length(ibx); %% number of doubly bounded linear inequalities 0221 0222 %% set up options struct for CPLEX 0223 if ~isempty(opt) && isfield(opt, 'cplex_opt') && ~isempty(opt.cplex_opt) 0224 cplex_opt = cplex_options(opt.cplex_opt); 0225 else 0226 cplex_opt = cplex_options; 0227 end 0228 0229 cplex = Cplex('null'); 0230 vstr = cplex.getVersion; 0231 [s,e,tE,m,t] = regexp(vstr, '(\d+\.\d+)\.'); 0232 vnum = str2num(t{1}{1}); 0233 vrb = max([0 verbose-1]); 0234 cplex_opt.barrier.display = vrb; 0235 cplex_opt.conflict.display = vrb; 0236 cplex_opt.mip.display = vrb; 0237 cplex_opt.sifting.display = vrb; 0238 cplex_opt.simplex.display = vrb; 0239 cplex_opt.tune.display = vrb; 0240 if vrb && vnum > 12.2 0241 cplex_opt.diagnostics = 'on'; 0242 end 0243 % if max_it 0244 % cplex_opt. %% not sure what to set here 0245 % end 0246 0247 if isempty(Ai) && isempty(Ae) 0248 unconstrained = 1; 0249 Ae = sparse(1, nx); 0250 be = 0; 0251 else 0252 unconstrained = 0; 0253 end 0254 0255 %% call the solver 0256 if verbose 0257 methods = { 0258 'default', 0259 'primal simplex', 0260 'dual simplex', 0261 'network simplex', 0262 'barrier', 0263 'sifting', 0264 'concurrent' 0265 }; 0266 end 0267 if isempty(H) || ~any(any(H)) 0268 if verbose 0269 fprintf('CPLEX Version %s -- %s LP solver\n', ... 0270 vstr, methods{cplex_opt.lpmethod+1}); 0271 end 0272 [x, f, eflag, output, lam] = ... 0273 cplexlp(c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0274 else 0275 if verbose 0276 fprintf('CPLEX Version %s -- %s QP solver\n', ... 0277 vstr, methods{cplex_opt.qpmethod+1}); 0278 end 0279 %% ensure H is numerically symmetric 0280 if ~isequal(H, H') 0281 H = (H + H')/2; 0282 end 0283 [x, f, eflag, output, lam] = ... 0284 cplexqp(H, c, Ai, bi, Ae, be, xmin, xmax, x0, cplex_opt); 0285 end 0286 0287 %% check for empty results (in case optimization failed) 0288 if isempty(x) 0289 x = NaN(nx, 1); 0290 end 0291 if isempty(f) 0292 f = NaN; 0293 end 0294 if isempty(lam) 0295 lam.ineqlin = NaN(length(bi), 1); 0296 lam.eqlin = NaN(length(be), 1); 0297 lam.lower = NaN(nx, 1); 0298 lam.upper = NaN(nx, 1); 0299 mu_l = NaN(nA, 1); 0300 mu_u = NaN(nA, 1); 0301 else 0302 mu_l = zeros(nA, 1); 0303 mu_u = zeros(nA, 1); 0304 end 0305 if unconstrained 0306 lam.eqlin = []; 0307 end 0308 0309 %% negate prices depending on version 0310 if vnum < 12.3 0311 lam.eqlin = -lam.eqlin; 0312 lam.ineqlin = -lam.ineqlin; 0313 end 0314 0315 %% repackage lambdas 0316 kl = find(lam.eqlin < 0); %% lower bound binding 0317 ku = find(lam.eqlin > 0); %% upper bound binding 0318 0319 mu_l(ieq(kl)) = -lam.eqlin(kl); 0320 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0321 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0322 0323 mu_u(ieq(ku)) = lam.eqlin(ku); 0324 mu_u(ilt) = lam.ineqlin(1:nlt); 0325 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0326 0327 lambda = struct( ... 0328 'mu_l', mu_l, ... 0329 'mu_u', mu_u, ... 0330 'lower', lam.lower, ... 0331 'upper', lam.upper ... 0332 );