QPS_BPMPD Quadratic Program Solver based on BPMPD_MEX.
[X, F, EXITFLAG, OUTPUT, LAMBDA] = ...
QPS_BPMPD(H, C, A, L, U, XMIN, XMAX, X0, OPT)
A wrapper function providing a MATPOWER standardized interface for using
BPMPD_MEX (http://www.pserc.cornell.edu/bpmpd/) 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
bp_opt - options vector for BP, 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 : exit flag,
1 = optimal solution
-1 = suboptimal solution
-2 = infeasible primal
-3 = infeasible dual
-10 = not enough memory
-99 = BPMPD bug: returned infeasible solution
OUTPUT : output struct with the following fields:
message - exit message
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_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_bpmpd(H, c, A, l, u)
x = qps_bpmpd(H, c, A, l, u, xmin, xmax)
x = qps_bpmpd(H, c, A, l, u, xmin, xmax, x0)
x = qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt)
x = qps_bpmpd(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_bpmpd(...)
[x, f] = qps_bpmpd(...)
[x, f, exitflag] = qps_bpmpd(...)
[x, f, exitflag, output] = qps_bpmpd(...)
[x, f, exitflag, output, lambda] = qps_bpmpd(...)
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_bpmpd(H, c, A, l, u, xmin, [], x0, opt);
See also BPMPD_MEX, http://www.pserc.cornell.edu/bpmpd/.0001 function [x, f, eflag, output, lambda] = qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt) 0002 %QPS_BPMPD Quadratic Program Solver based on BPMPD_MEX. 0003 % [X, F, EXITFLAG, OUTPUT, LAMBDA] = ... 0004 % QPS_BPMPD(H, C, A, L, U, XMIN, XMAX, X0, OPT) 0005 % A wrapper function providing a MATPOWER standardized interface for using 0006 % BPMPD_MEX (http://www.pserc.cornell.edu/bpmpd/) to solve the 0007 % following QP (quadratic programming) 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 % bp_opt - options vector for BP, 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 : exit flag, 0043 % 1 = optimal solution 0044 % -1 = suboptimal solution 0045 % -2 = infeasible primal 0046 % -3 = infeasible dual 0047 % -10 = not enough memory 0048 % -99 = BPMPD bug: returned infeasible solution 0049 % OUTPUT : output struct with the following fields: 0050 % message - exit message 0051 % LAMBDA : struct containing the Langrange and Kuhn-Tucker 0052 % multipliers on the constraints, with fields: 0053 % mu_l - lower (left-hand) limit on linear constraints 0054 % mu_u - upper (right-hand) limit on linear constraints 0055 % lower - lower bound on optimization variables 0056 % upper - upper bound on optimization variables 0057 % 0058 % Note the calling syntax is almost identical to that of QUADPROG 0059 % from MathWorks' Optimization Toolbox. The main difference is that 0060 % the linear constraints are specified with A, L, U instead of 0061 % A, B, Aeq, Beq. 0062 % 0063 % Calling syntax options: 0064 % [x, f, exitflag, output, lambda] = ... 0065 % qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt) 0066 % 0067 % x = qps_bpmpd(H, c, A, l, u) 0068 % x = qps_bpmpd(H, c, A, l, u, xmin, xmax) 0069 % x = qps_bpmpd(H, c, A, l, u, xmin, xmax, x0) 0070 % x = qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt) 0071 % x = qps_bpmpd(problem), where problem is a struct with fields: 0072 % H, c, A, l, u, xmin, xmax, x0, opt 0073 % all fields except 'c', 'A' and 'l' or 'u' are optional 0074 % x = qps_bpmpd(...) 0075 % [x, f] = qps_bpmpd(...) 0076 % [x, f, exitflag] = qps_bpmpd(...) 0077 % [x, f, exitflag, output] = qps_bpmpd(...) 0078 % [x, f, exitflag, output, lambda] = qps_bpmpd(...) 0079 % 0080 % Example: (problem from from http://www.jmu.edu/docs/sasdoc/sashtml/iml/chap8/sect12.htm) 0081 % H = [ 1003.1 4.3 6.3 5.9; 0082 % 4.3 2.2 2.1 3.9; 0083 % 6.3 2.1 3.5 4.8; 0084 % 5.9 3.9 4.8 10 ]; 0085 % c = zeros(4,1); 0086 % A = [ 1 1 1 1; 0087 % 0.17 0.11 0.10 0.18 ]; 0088 % l = [1; 0.10]; 0089 % u = [1; Inf]; 0090 % xmin = zeros(4,1); 0091 % x0 = [1; 0; 0; 1]; 0092 % opt = struct('verbose', 2); 0093 % [x, f, s, out, lambda] = qps_bpmpd(H, c, A, l, u, xmin, [], x0, opt); 0094 % 0095 % See also BPMPD_MEX, http://www.pserc.cornell.edu/bpmpd/. 0096 0097 % MATPOWER 0098 % Copyright (c) 2010-2015 by Power System Engineering Research Center (PSERC) 0099 % by Ray Zimmerman, PSERC Cornell 0100 % 0101 % $Id: qps_bpmpd.m 2661 2015-03-20 17:02:46Z ray $ 0102 % 0103 % This file is part of MATPOWER. 0104 % Covered by the 3-clause BSD License (see LICENSE file for details). 0105 % See http://www.pserc.cornell.edu/matpower/ for more info. 0106 0107 %% check for BPMPD_MEX 0108 % if ~have_fcn('bpmpd') 0109 % error('qps_bpmpd: requires BPMPD_MEX (http://www.pserc.cornell.edu/bpmpd/)'); 0110 % end 0111 0112 %%----- input argument handling ----- 0113 %% gather inputs 0114 if nargin == 1 && isstruct(H) %% problem struct 0115 p = H; 0116 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0117 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0118 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0119 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0120 if isfield(p, 'u'), u = p.u; else, u = []; end 0121 if isfield(p, 'l'), l = p.l; else, l = []; end 0122 if isfield(p, 'A'), A = p.A; else, A = []; end 0123 if isfield(p, 'c'), c = p.c; else, c = []; end 0124 if isfield(p, 'H'), H = p.H; else, H = []; end 0125 else %% individual args 0126 if nargin < 9 0127 opt = []; 0128 if nargin < 8 0129 x0 = []; 0130 if nargin < 7 0131 xmax = []; 0132 if nargin < 6 0133 xmin = []; 0134 end 0135 end 0136 end 0137 end 0138 end 0139 0140 %% define nx, set default values for missing optional inputs 0141 if isempty(H) || ~any(any(H)) 0142 if isempty(A) && isempty(xmin) && isempty(xmax) 0143 error('qps_bpmpd: LP problem must include constraints or variable bounds'); 0144 else 0145 if ~isempty(A) 0146 nx = size(A, 2); 0147 elseif ~isempty(xmin) 0148 nx = length(xmin); 0149 else % if ~isempty(xmax) 0150 nx = length(xmax); 0151 end 0152 end 0153 else 0154 nx = size(H, 1); 0155 end 0156 if isempty(c) 0157 c = zeros(nx, 1); 0158 end 0159 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0160 (isempty(u) || all(u == Inf))) 0161 A = sparse(0,nx); %% no limits => no linear constraints 0162 end 0163 nA = size(A, 1); %% number of original linear constraints 0164 if isempty(u) %% By default, linear inequalities are ... 0165 u = Inf(nA, 1); %% ... unbounded above and ... 0166 end 0167 if isempty(l) 0168 l = -Inf(nA, 1); %% ... unbounded below. 0169 end 0170 if isempty(xmin) %% By default, optimization variables are ... 0171 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0172 end 0173 if isempty(xmax) 0174 xmax = Inf(nx, 1); %% ... unbounded above. 0175 end 0176 if isempty(x0) 0177 x0 = zeros(nx, 1); 0178 end 0179 0180 %% default options 0181 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0182 verbose = opt.verbose; 0183 else 0184 verbose = 0; 0185 end 0186 0187 %% make sure args are sparse/full as expected by BPMPD 0188 if ~isempty(H) 0189 if ~issparse(H) 0190 H = sparse(H); 0191 end 0192 end 0193 if ~issparse(A) 0194 A = sparse(A); 0195 end 0196 if issparse(c) 0197 c = full(c); %% avoid a crash 0198 end 0199 0200 %% split up linear constraints 0201 ieq = find( abs(u-l) <= eps ); %% equality 0202 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0203 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0204 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0205 Ae = A(ieq, :); 0206 be = u(ieq); 0207 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0208 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0209 0210 %% grab some dimensions 0211 neq = length(ieq); %% number of equality constraints 0212 niq = length(bi); %% number of inequality constraints 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 linear constraints 0218 if neq || niq 0219 AA = [Ae; Ai]; 0220 bb = [be; bi]; 0221 ee = [zeros(neq, 1); -ones(niq, 1)]; 0222 else 0223 AA = []; bb = []; ee = []; 0224 end 0225 0226 %% set up variable bounds and initial value 0227 if ~isempty(xmin) 0228 llist = find(xmin > -1e15); % interpret limits <= -1e15 as unbounded 0229 if isempty(llist) 0230 llist = []; 0231 lval = []; 0232 else 0233 lval = xmin(llist); 0234 end 0235 else 0236 llist = []; 0237 lval = []; 0238 end 0239 if ~isempty(xmax) 0240 ulist = find(xmax < 1e15); % interpret limits >= 1e15 as unbounded 0241 if isempty(ulist) 0242 ulist = []; 0243 uval = []; 0244 else 0245 uval = xmax(ulist); 0246 end 0247 else 0248 ulist = []; 0249 uval = []; 0250 end 0251 0252 %% set up options 0253 if ~isempty(opt) && isfield(opt, 'bp_opt') && ~isempty(opt.bp_opt) 0254 bp_opt = opt.bp_opt; 0255 else 0256 bp_opt = bpopt; %% use default options 0257 % bp_opt(14)= 1e-3; % TPIV1 first relative pivot tolerance (desired) 0258 % bp_opt(20)= 1e-8; % TOPT1 stop if feasible and rel. dual gap less than this 0259 % bp_opt(22)= 1e-7; % TFEAS1 relative primal feasibility tolerance 0260 % bp_opt(23)= 1e-7; % TFEAS2 relative dual feasibility tolerance 0261 % bp_opt(29)= 1e-9; % TRESX acceptable primal residual 0262 % bp_opt(30)= 1e-9; % TRESY acceptable dual residual 0263 % bp_opt(38)= 2; % SMETHOD1 prescaling method 0264 end 0265 if verbose > 1 0266 prnlev = 1; 0267 else 0268 prnlev = 0; 0269 end 0270 if strcmp(computer, 'PCWIN') 0271 if prnlev 0272 fprintf('Windows version of BPMPD_MEX cannot print to screen.\n'); 0273 end 0274 prnlev = 0; % The Windows incarnation of bp was born mute and deaf, 0275 end % probably because of acute shock after realizing its fate. 0276 % Can't be allowed to try to speak or its universe crumbles. 0277 0278 %% call the solver 0279 [x, y, s, w, output.message] = bp(H, AA, bb, c, ee, llist, lval, ... 0280 ulist, uval, bp_opt, prnlev); 0281 0282 %% compute final objective 0283 if nargout > 1 0284 f = 0; 0285 if ~isempty(c) 0286 f = f + c' * x; 0287 end 0288 if ~isempty(H) 0289 f = f + 0.5 * x' * H * x; 0290 end 0291 end 0292 0293 %% set exit flag 0294 if strcmp(output.message, 'optimal solution') 0295 eflag = 1; 0296 elseif strcmp(output.message, 'suboptimal solution') 0297 eflag = -1; 0298 elseif strcmp(output.message, 'infeasible primal') 0299 eflag = -2; 0300 elseif strcmp(output.message, 'infeasible dual') 0301 eflag = -3; 0302 elseif strcmp(output.message, 'not enough memory') 0303 eflag = -10; 0304 else 0305 eflag = 0; 0306 end 0307 0308 %% zero out lambdas smaller than a certain tolerance 0309 y(abs(y) < 1e-9) = 0; 0310 w(abs(w) < 1e-9) = 0; 0311 0312 %% necessary for proper operation of constr.m 0313 if eflag == -2 %% infeasible primal 0314 y = zeros(size(y)); 0315 w = zeros(size(w)); 0316 end 0317 0318 %% repackage lambdas 0319 lam.eqlin = -y(1:neq); 0320 lam.ineqlin = -y(neq+(1:niq)); 0321 kl = find(lam.eqlin < 0); %% lower bound binding 0322 ku = find(lam.eqlin > 0); %% upper bound binding 0323 0324 mu_l = zeros(nA, 1); 0325 mu_l(ieq(kl)) = -lam.eqlin(kl); 0326 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0327 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0328 0329 mu_u = zeros(nA, 1); 0330 mu_u(ieq(ku)) = lam.eqlin(ku); 0331 mu_u(ilt) = lam.ineqlin(1:nlt); 0332 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0333 0334 lam.lower = zeros(nx, 1); 0335 lam.upper = zeros(nx, 1); 0336 kl = find(w > 0); %% lower bound binding 0337 ku = find(w < 0); %% upper bound binding 0338 lam.lower(kl) = w(kl); 0339 lam.upper(ku) = -w(ku); 0340 0341 lambda = struct( ... 0342 'mu_l', mu_l, ... 0343 'mu_u', mu_u, ... 0344 'lower', lam.lower, ... 0345 'upper', lam.upper ... 0346 ); 0347 0348 %% Note: BPMPD_MEX has a bug which causes it to return an incorrect 0349 %% (infeasible) solution for some problems. 0350 %% So we need to double-check for feasibility 0351 if eflag > 0 0352 bpmpd_bug_fatal = 0; 0353 err_tol = 5e-4; 0354 if ~isempty(xmin) 0355 lb_violation = xmin - x; 0356 if verbose > 1 0357 fprintf('max variable lower bound violatation: %g\n', max(lb_violation)); 0358 end 0359 else 0360 lb_violation = zeros(nx, 1); 0361 end 0362 if ~isempty(xmax) 0363 ub_violation = x - xmax; 0364 if verbose > 1 0365 fprintf('max variable upper bound violation: %g\n', max(ub_violation)); 0366 end 0367 else 0368 ub_violation = zeros(nx, 1); 0369 end 0370 if neq > 0 0371 eq_violation = abs( Ae * x - be ); 0372 if verbose > 1 0373 fprintf('max equality constraint violation: %g\n', max(eq_violation)); 0374 end 0375 else 0376 eq_violation = zeros(neq, 1); 0377 end 0378 if niq 0379 ineq_violation = Ai * x - bi; 0380 if verbose > 1 0381 fprintf('max inequality constraint violation: %g\n', max(ineq_violation)); 0382 end 0383 else 0384 ineq_violation = zeros(niq, 1); 0385 end 0386 if any( [ lb_violation; 0387 ub_violation; 0388 eq_violation; 0389 ineq_violation ] > err_tol) 0390 err_cnt = 0; 0391 if any( lb_violation > err_tol ) 0392 err_cnt = err_cnt + 1; 0393 errs{err_cnt} = ... 0394 sprintf('variable lower bound violated by %g', ... 0395 max(lb_violation)); 0396 end 0397 if any( ub_violation > err_tol ) 0398 err_cnt = err_cnt + 1; 0399 errs{err_cnt} = ... 0400 sprintf('variable upper bound violated by %g', ... 0401 max(ub_violation)); 0402 end 0403 if any( eq_violation > err_tol ) 0404 err_cnt = err_cnt + 1; 0405 errs{err_cnt} = ... 0406 sprintf('equality constraint violated by %g', ... 0407 max(eq_violation)); 0408 end 0409 if any( ineq_violation > err_tol ) 0410 err_cnt = err_cnt + 1; 0411 errs{err_cnt} = ... 0412 sprintf('inequality constraint violated by %g', ... 0413 max(ineq_violation)); 0414 end 0415 if verbose > 0 0416 fprintf('\nWARNING: This version of BPMPD_MEX has a bug which caused it to return\n'); 0417 fprintf( ' an incorrect (infeasible) solution for this particular problem.\n'); 0418 end 0419 for err_idx = 1:err_cnt 0420 fprintf(' %s\n', errs{err_idx}); 0421 end 0422 if bpmpd_bug_fatal 0423 error('%s\n%s', ... 0424 'To avoid this bug in BPMPD_MEX you will need', ... 0425 'to use a different QP solver for this problem.'); 0426 end 0427 eflag = -99; 0428 output.message = [output.message '\nBPMPD bug: returned infeasible solution']; 0429 end 0430 end