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