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 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 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 https://v8doc.sas.com/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-2016, Power Systems Engineering Research Center (PSERC) 0099 % by Ray Zimmerman, PSERC Cornell 0100 % 0101 % This file is part of MATPOWER. 0102 % Covered by the 3-clause BSD License (see LICENSE file for details). 0103 % See https://matpower.org for more info. 0104 0105 %% check for BPMPD_MEX 0106 % if ~have_fcn('bpmpd') 0107 % error('qps_bpmpd: requires BPMPD_MEX (http://www.pserc.cornell.edu/bpmpd/)'); 0108 % end 0109 0110 %%----- input argument handling ----- 0111 %% gather inputs 0112 if nargin == 1 && isstruct(H) %% problem struct 0113 p = H; 0114 if isfield(p, 'opt'), opt = p.opt; else, opt = []; end 0115 if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end 0116 if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end 0117 if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end 0118 if isfield(p, 'u'), u = p.u; else, u = []; end 0119 if isfield(p, 'l'), l = p.l; else, l = []; end 0120 if isfield(p, 'A'), A = p.A; else, A = []; end 0121 if isfield(p, 'c'), c = p.c; else, c = []; end 0122 if isfield(p, 'H'), H = p.H; else, H = []; end 0123 else %% individual args 0124 if nargin < 9 0125 opt = []; 0126 if nargin < 8 0127 x0 = []; 0128 if nargin < 7 0129 xmax = []; 0130 if nargin < 6 0131 xmin = []; 0132 end 0133 end 0134 end 0135 end 0136 end 0137 0138 %% define nx, set default values for missing optional inputs 0139 if isempty(H) || ~any(any(H)) 0140 if isempty(A) && isempty(xmin) && isempty(xmax) 0141 error('qps_bpmpd: LP problem must include constraints or variable bounds'); 0142 else 0143 if ~isempty(A) 0144 nx = size(A, 2); 0145 elseif ~isempty(xmin) 0146 nx = length(xmin); 0147 else % if ~isempty(xmax) 0148 nx = length(xmax); 0149 end 0150 end 0151 else 0152 nx = size(H, 1); 0153 end 0154 if isempty(c) 0155 c = zeros(nx, 1); 0156 end 0157 if isempty(A) || (~isempty(A) && (isempty(l) || all(l == -Inf)) && ... 0158 (isempty(u) || all(u == Inf))) 0159 A = sparse(0,nx); %% no limits => no linear constraints 0160 end 0161 nA = size(A, 1); %% number of original linear constraints 0162 if isempty(u) %% By default, linear inequalities are ... 0163 u = Inf(nA, 1); %% ... unbounded above and ... 0164 end 0165 if isempty(l) 0166 l = -Inf(nA, 1); %% ... unbounded below. 0167 end 0168 if isempty(xmin) %% By default, optimization variables are ... 0169 xmin = -Inf(nx, 1); %% ... unbounded below and ... 0170 end 0171 if isempty(xmax) 0172 xmax = Inf(nx, 1); %% ... unbounded above. 0173 end 0174 if isempty(x0) 0175 x0 = zeros(nx, 1); 0176 end 0177 0178 %% default options 0179 if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose) 0180 verbose = opt.verbose; 0181 else 0182 verbose = 0; 0183 end 0184 0185 %% make sure args are sparse/full as expected by BPMPD 0186 if ~isempty(H) 0187 if ~issparse(H) 0188 H = sparse(H); 0189 end 0190 end 0191 if ~issparse(A) 0192 A = sparse(A); 0193 end 0194 if issparse(c) 0195 c = full(c); %% avoid a crash 0196 end 0197 0198 %% split up linear constraints 0199 ieq = find( abs(u-l) <= eps ); %% equality 0200 igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above 0201 ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below 0202 ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) ); 0203 Ae = A(ieq, :); 0204 be = u(ieq); 0205 Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ]; 0206 bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)]; 0207 0208 %% grab some dimensions 0209 neq = length(ieq); %% number of equality constraints 0210 niq = length(bi); %% number of inequality constraints 0211 nlt = length(ilt); %% number of upper bounded linear inequalities 0212 ngt = length(igt); %% number of lower bounded linear inequalities 0213 nbx = length(ibx); %% number of doubly bounded linear inequalities 0214 0215 %% set up linear constraints 0216 if neq || niq 0217 AA = [Ae; Ai]; 0218 bb = [be; bi]; 0219 ee = [zeros(neq, 1); -ones(niq, 1)]; 0220 else 0221 AA = []; bb = []; ee = []; 0222 end 0223 0224 %% set up variable bounds and initial value 0225 if ~isempty(xmin) 0226 llist = find(xmin > -1e15); % interpret limits <= -1e15 as unbounded 0227 if isempty(llist) 0228 llist = []; 0229 lval = []; 0230 else 0231 lval = xmin(llist); 0232 end 0233 else 0234 llist = []; 0235 lval = []; 0236 end 0237 if ~isempty(xmax) 0238 ulist = find(xmax < 1e15); % interpret limits >= 1e15 as unbounded 0239 if isempty(ulist) 0240 ulist = []; 0241 uval = []; 0242 else 0243 uval = xmax(ulist); 0244 end 0245 else 0246 ulist = []; 0247 uval = []; 0248 end 0249 0250 %% set up options 0251 if ~isempty(opt) && isfield(opt, 'bp_opt') && ~isempty(opt.bp_opt) 0252 bp_opt = opt.bp_opt; 0253 else 0254 bp_opt = bpopt; %% use default options 0255 % bp_opt(14)= 1e-3; % TPIV1 first relative pivot tolerance (desired) 0256 % bp_opt(20)= 1e-8; % TOPT1 stop if feasible and rel. dual gap less than this 0257 % bp_opt(22)= 1e-7; % TFEAS1 relative primal feasibility tolerance 0258 % bp_opt(23)= 1e-7; % TFEAS2 relative dual feasibility tolerance 0259 % bp_opt(29)= 1e-9; % TRESX acceptable primal residual 0260 % bp_opt(30)= 1e-9; % TRESY acceptable dual residual 0261 % bp_opt(38)= 2; % SMETHOD1 prescaling method 0262 end 0263 if verbose > 1 0264 prnlev = 1; 0265 else 0266 prnlev = 0; 0267 end 0268 if strcmp(computer, 'PCWIN') 0269 if prnlev 0270 fprintf('Windows version of BPMPD_MEX cannot print to screen.\n'); 0271 end 0272 prnlev = 0; % The Windows incarnation of bp was born mute and deaf, 0273 end % probably because of acute shock after realizing its fate. 0274 % Can't be allowed to try to speak or its universe crumbles. 0275 0276 %% call the solver 0277 [x, y, s, w, output.message] = bp(H, AA, bb, c, ee, llist, lval, ... 0278 ulist, uval, bp_opt, prnlev); 0279 0280 %% compute final objective 0281 if nargout > 1 0282 f = 0; 0283 if ~isempty(c) 0284 f = f + c' * x; 0285 end 0286 if ~isempty(H) 0287 f = f + 0.5 * x' * H * x; 0288 end 0289 end 0290 0291 %% set exit flag 0292 if strcmp(output.message, 'optimal solution') 0293 eflag = 1; 0294 elseif strcmp(output.message, 'suboptimal solution') 0295 eflag = -1; 0296 elseif strcmp(output.message, 'infeasible primal') 0297 eflag = -2; 0298 elseif strcmp(output.message, 'infeasible dual') 0299 eflag = -3; 0300 elseif strcmp(output.message, 'not enough memory') 0301 eflag = -10; 0302 else 0303 eflag = 0; 0304 end 0305 0306 %% zero out lambdas smaller than a certain tolerance 0307 y(abs(y) < 1e-9) = 0; 0308 w(abs(w) < 1e-9) = 0; 0309 0310 %% necessary for proper operation of constr.m 0311 if eflag == -2 %% infeasible primal 0312 y = zeros(size(y)); 0313 w = zeros(size(w)); 0314 end 0315 0316 %% repackage lambdas 0317 lam.eqlin = -y(1:neq); 0318 lam.ineqlin = -y(neq+(1:niq)); 0319 kl = find(lam.eqlin < 0); %% lower bound binding 0320 ku = find(lam.eqlin > 0); %% upper bound binding 0321 0322 mu_l = zeros(nA, 1); 0323 mu_l(ieq(kl)) = -lam.eqlin(kl); 0324 mu_l(igt) = lam.ineqlin(nlt+(1:ngt)); 0325 mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx)); 0326 0327 mu_u = zeros(nA, 1); 0328 mu_u(ieq(ku)) = lam.eqlin(ku); 0329 mu_u(ilt) = lam.ineqlin(1:nlt); 0330 mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx)); 0331 0332 lam.lower = zeros(nx, 1); 0333 lam.upper = zeros(nx, 1); 0334 kl = find(w > 0); %% lower bound binding 0335 ku = find(w < 0); %% upper bound binding 0336 lam.lower(kl) = w(kl); 0337 lam.upper(ku) = -w(ku); 0338 0339 lambda = struct( ... 0340 'mu_l', mu_l, ... 0341 'mu_u', mu_u, ... 0342 'lower', lam.lower, ... 0343 'upper', lam.upper ... 0344 ); 0345 0346 %% Note: BPMPD_MEX has a bug which causes it to return an incorrect 0347 %% (infeasible) solution for some problems. 0348 %% So we need to double-check for feasibility 0349 if eflag > 0 0350 bpmpd_bug_fatal = 0; 0351 err_tol = 5e-4; 0352 if ~isempty(xmin) 0353 lb_violation = xmin - x; 0354 if verbose > 1 0355 fprintf('max variable lower bound violatation: %g\n', max(lb_violation)); 0356 end 0357 else 0358 lb_violation = zeros(nx, 1); 0359 end 0360 if ~isempty(xmax) 0361 ub_violation = x - xmax; 0362 if verbose > 1 0363 fprintf('max variable upper bound violation: %g\n', max(ub_violation)); 0364 end 0365 else 0366 ub_violation = zeros(nx, 1); 0367 end 0368 if neq > 0 0369 eq_violation = abs( Ae * x - be ); 0370 if verbose > 1 0371 fprintf('max equality constraint violation: %g\n', max(eq_violation)); 0372 end 0373 else 0374 eq_violation = zeros(neq, 1); 0375 end 0376 if niq 0377 ineq_violation = Ai * x - bi; 0378 if verbose > 1 0379 fprintf('max inequality constraint violation: %g\n', max(ineq_violation)); 0380 end 0381 else 0382 ineq_violation = zeros(niq, 1); 0383 end 0384 if any( [ lb_violation; 0385 ub_violation; 0386 eq_violation; 0387 ineq_violation ] > err_tol) 0388 err_cnt = 0; 0389 if any( lb_violation > err_tol ) 0390 err_cnt = err_cnt + 1; 0391 errs{err_cnt} = ... 0392 sprintf('variable lower bound violated by %g', ... 0393 max(lb_violation)); 0394 end 0395 if any( ub_violation > err_tol ) 0396 err_cnt = err_cnt + 1; 0397 errs{err_cnt} = ... 0398 sprintf('variable upper bound violated by %g', ... 0399 max(ub_violation)); 0400 end 0401 if any( eq_violation > err_tol ) 0402 err_cnt = err_cnt + 1; 0403 errs{err_cnt} = ... 0404 sprintf('equality constraint violated by %g', ... 0405 max(eq_violation)); 0406 end 0407 if any( ineq_violation > err_tol ) 0408 err_cnt = err_cnt + 1; 0409 errs{err_cnt} = ... 0410 sprintf('inequality constraint violated by %g', ... 0411 max(ineq_violation)); 0412 end 0413 if verbose > 0 0414 fprintf('\nWARNING: This version of BPMPD_MEX has a bug which caused it to return\n'); 0415 fprintf( ' an incorrect (infeasible) solution for this particular problem.\n'); 0416 end 0417 for err_idx = 1:err_cnt 0418 fprintf(' %s\n', errs{err_idx}); 0419 end 0420 if bpmpd_bug_fatal 0421 error('%s\n%s', ... 0422 'To avoid this bug in BPMPD_MEX you will need', ... 0423 'to use a different QP solver for this problem.'); 0424 end 0425 eflag = -99; 0426 output.message = [output.message '\nBPMPD bug: returned infeasible solution']; 0427 end 0428 end