Home > matpower7.0 > 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)
   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/.

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

Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005