Home > matpower7.1 > mp-opt-model > lib > 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)
   [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/.

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 %   [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

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005