Home > matpower5.0 > 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 %   $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

Generated on Mon 26-Jan-2015 15:21:31 by m2html © 2005