Home > matpower4.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
           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, 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 %           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, lambda] = 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.11 2011/09/09 15:26:08 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         if verbose > 0
0445             fprintf('\nWARNING: This version of BPMPD_MEX has a bug which caused it to return\n');
0446             fprintf(  '         an incorrect (infeasible) solution for this particular problem.\n');
0447         end
0448         for err_idx = 1:err_cnt
0449             fprintf('         %s\n', errs{err_idx});
0450         end
0451         if bpmpd_bug_fatal
0452             error('%s\n%s', ...
0453                 'To avoid this bug in BPMPD_MEX you will need', ...
0454                 'to use a different QP solver for this problem.');
0455         end
0456         eflag = -99;
0457         output.message = [output.message '\nBPMPD bug: returned infeasible solution'];
0458     end
0459 end

Generated on Mon 26-Jan-2015 15:00:13 by m2html © 2005