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