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