Home > matpower5.1 > qps_bpmpd.m

qps_bpmpd

PURPOSE ^

QPS_BPMPD Quadratic Program Solver based on BPMPD_MEX.

SYNOPSIS ^

function [x, f, eflag, output, lambda] = qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt)

DESCRIPTION ^

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/.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005