www.gusucode.com > 基于lingo求所以解，对潮流计算求出所有解 > matpower4.1/qps_bpmpd.m
function [x, f, eflag, output, lambda] = qps_bpmpd(H, c, A, l, u, xmin, xmax, x0, opt)
%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/.
% MATPOWER
% $Id: qps_bpmpd.m,v 1.11 2011/09/09 15:26:08 cvs Exp $
% by Ray Zimmerman, PSERC Cornell
% Copyright (c) 2010 by Power System Engineering Research Center (PSERC)
%
% This file is part of MATPOWER.
% See http://www.pserc.cornell.edu/matpower/ for more info.
%
% MATPOWER is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published
% by the Free Software Foundation, either version 3 of the License,
% or (at your option) any later version.
%
% MATPOWER is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
%
% Additional permission under GNU GPL version 3 section 7
%
% If you modify MATPOWER, or any covered work, to interface with
% other modules (such as MATLAB code and MEX-files) available in a
% MATLAB(R) or comparable environment containing parts covered
% under other licensing terms, the licensors of MATPOWER grant
% you additional permission to convey the resulting work.
%% check for BPMPD_MEX
% if ~have_fcn('bpmpd')
% error('qps_bpmpd: requires BPMPD_MEX (http://www.pserc.cornell.edu/bpmpd/)');
% end
%%----- input argument handling -----
%% gather inputs
if nargin == 1 && isstruct(H) %% problem struct
p = H;
if isfield(p, 'opt'), opt = p.opt; else, opt = []; end
if isfield(p, 'x0'), x0 = p.x0; else, x0 = []; end
if isfield(p, 'xmax'), xmax = p.xmax; else, xmax = []; end
if isfield(p, 'xmin'), xmin = p.xmin; else, xmin = []; end
if isfield(p, 'u'), u = p.u; else, u = []; end
if isfield(p, 'l'), l = p.l; else, l = []; end
if isfield(p, 'A'), A = p.A; else, A = []; end
if isfield(p, 'c'), c = p.c; else, c = []; end
if isfield(p, 'H'), H = p.H; else, H = []; end
else %% individual args
if nargin < 9
opt = [];
if nargin < 8
x0 = [];
if nargin < 7
xmax = [];
if nargin < 6
xmin = [];
end
end
end
end
end
%% define nx, set default values for missing optional inputs
if isempty(H) || ~any(any(H))
if isempty(A) && isempty(xmin) && isempty(xmax)
error('qps_bpmpd: LP problem must include constraints or variable bounds');
else
if ~isempty(A)
nx = size(A, 2);
elseif ~isempty(xmin)
nx = length(xmin);
else % if ~isempty(xmax)
nx = length(xmax);
end
end
else
nx = size(H, 1);
end
if isempty(c)
c = zeros(nx, 1);
end
if ~isempty(A) && (isempty(l) || all(l == -Inf)) && ...
(isempty(u) || all(u == Inf))
A = sparse(0,nx); %% no limits => no linear constraints
end
nA = size(A, 1); %% number of original linear constraints
if isempty(u) %% By default, linear inequalities are ...
u = Inf * ones(nA, 1); %% ... unbounded above and ...
end
if isempty(l)
l = -Inf * ones(nA, 1); %% ... unbounded below.
end
if isempty(xmin) %% By default, optimization variables are ...
xmin = -Inf * ones(nx, 1); %% ... unbounded below and ...
end
if isempty(xmax)
xmax = Inf * ones(nx, 1); %% ... unbounded above.
end
if isempty(x0)
x0 = zeros(nx, 1);
end
%% default options
if ~isempty(opt) && isfield(opt, 'verbose') && ~isempty(opt.verbose)
verbose = opt.verbose;
else
verbose = 0;
end
if ~isempty(opt) && isfield(opt, 'max_it') && ~isempty(opt.max_it)
max_it = opt.max_it;
else
max_it = 0;
end
%% make sure args are sparse/full as expected by BPMPD
if ~isempty(H)
if ~issparse(H)
H = sparse(H);
end
end
if ~issparse(A)
A = sparse(A);
end
if issparse(c)
c = full(c); %% avoid a crash
end
%% split up linear constraints
ieq = find( abs(u-l) <= eps ); %% equality
igt = find( u >= 1e10 & l > -1e10 ); %% greater than, unbounded above
ilt = find( l <= -1e10 & u < 1e10 ); %% less than, unbounded below
ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10) );
Ae = A(ieq, :);
be = u(ieq);
Ai = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
bi = [ u(ilt); -l(igt); u(ibx); -l(ibx)];
%% grab some dimensions
neq = length(ieq); %% number of equality constraints
niq = length(bi); %% number of inequality constraints
nlt = length(ilt); %% number of upper bounded linear inequalities
ngt = length(igt); %% number of lower bounded linear inequalities
nbx = length(ibx); %% number of doubly bounded linear inequalities
%% set up linear constraints
if neq || niq
AA = [Ae; Ai];
bb = [be; bi];
ee = [zeros(neq, 1); -ones(niq, 1)];
else
AA = []; bb = []; ee = [];
end
%% set up variable bounds and initial value
if ~isempty(xmin)
llist = find(xmin > -1e15); % interpret limits <= -1e15 as unbounded
if isempty(llist)
llist = [];
lval = [];
else
lval = xmin(llist);
end
else
llist = [];
lval = [];
end
if ~isempty(xmax)
ulist = find(xmax < 1e15); % interpret limits >= 1e15 as unbounded
if isempty(ulist)
ulist = [];
uval = [];
else
uval = xmax(ulist);
end
else
ulist = [];
uval = [];
end
%% set up options
if ~isempty(opt) && isfield(opt, 'bp_opt') && ~isempty(opt.bp_opt)
bp_opt = opt.bp_opt;
else
bp_opt = bpopt; %% use default options
% bp_opt(14)= 1e-3; % TPIV1 first relative pivot tolerance (desired)
% bp_opt(20)= 1e-8; % TOPT1 stop if feasible and rel. dual gap less than this
% bp_opt(22)= 1e-7; % TFEAS1 relative primal feasibility tolerance
% bp_opt(23)= 1e-7; % TFEAS2 relative dual feasibility tolerance
% bp_opt(29)= 1e-9; % TRESX acceptable primal residual
% bp_opt(30)= 1e-9; % TRESY acceptable dual residual
% bp_opt(38)= 2; % SMETHOD1 prescaling method
end
if max_it
bp_opt(26) = max_it; %% MAXITER
end
if verbose > 1
prnlev = 1;
else
prnlev = 0;
end
if strcmp(computer, 'PCWIN')
if prnlev
fprintf('Windows version of BPMPD_MEX cannot print to screen.\n');
end
prnlev = 0; % The Windows incarnation of bp was born mute and deaf,
end % probably because of acute shock after realizing its fate.
% Can't be allowed to try to speak or its universe crumbles.
%% call the solver
[x, y, s, w, output.message] = bp(H, AA, bb, c, ee, llist, lval, ...
ulist, uval, bp_opt, prnlev);
%% compute final objective
if nargout > 1
f = 0;
if ~isempty(c)
f = f + c' * x;
end
if ~isempty(H)
f = f + 0.5 * x' * H * x;
end
end
%% set exit flag
if strcmp(output.message, 'optimal solution')
eflag = 1;
elseif strcmp(output.message, 'suboptimal solution')
eflag = -1;
elseif strcmp(output.message, 'infeasible primal')
eflag = -2;
elseif strcmp(output.message, 'infeasible dual')
eflag = -3;
elseif strcmp(output.message, 'not enough memory')
eflag = -10;
else
eflag = 0;
end
%% zero out lambdas smaller than a certain tolerance
y(abs(y) < 1e-9) = 0;
w(abs(w) < 1e-9) = 0;
%% necessary for proper operation of constr.m
if eflag == -2 %% infeasible primal
y = zeros(size(y));
w = zeros(size(w));
end
%% repackage lambdas
lam.eqlin = -y(1:neq);
lam.ineqlin = -y(neq+(1:niq));
kl = find(lam.eqlin < 0); %% lower bound binding
ku = find(lam.eqlin > 0); %% upper bound binding
mu_l = zeros(nA, 1);
mu_l(ieq(kl)) = -lam.eqlin(kl);
mu_l(igt) = lam.ineqlin(nlt+(1:ngt));
mu_l(ibx) = lam.ineqlin(nlt+ngt+nbx+(1:nbx));
mu_u = zeros(nA, 1);
mu_u(ieq(ku)) = lam.eqlin(ku);
mu_u(ilt) = lam.ineqlin(1:nlt);
mu_u(ibx) = lam.ineqlin(nlt+ngt+(1:nbx));
lam.lower = zeros(nx, 1);
lam.upper = zeros(nx, 1);
kl = find(w > 0); %% lower bound binding
ku = find(w < 0); %% upper bound binding
lam.lower(kl) = w(kl);
lam.upper(ku) = -w(ku);
lambda = struct( ...
'mu_l', mu_l, ...
'mu_u', mu_u, ...
'lower', lam.lower, ...
'upper', lam.upper ...
);
%% Note: BPMPD_MEX has a bug which causes it to return an incorrect
%% (infeasible) solution for some problems.
%% So we need to double-check for feasibility
if eflag > 0
bpmpd_bug_fatal = 0;
err_tol = 5e-4;
if ~isempty(xmin)
lb_violation = xmin - x;
if verbose > 1
fprintf('max variable lower bound violatation: %g\n', max(lb_violation));
end
else
lb_violation = zeros(nx, 1);
end
if ~isempty(xmax)
ub_violation = x - xmax;
if verbose > 1
fprintf('max variable upper bound violation: %g\n', max(ub_violation));
end
else
ub_violation = zeros(nx, 1);
end
if neq > 0
eq_violation = abs( Ae * x - be );
if verbose > 1
fprintf('max equality constraint violation: %g\n', max(eq_violation));
end
else
eq_violation = zeros(neq, 1);
end
if niq
ineq_violation = Ai * x - bi;
if verbose > 1
fprintf('max inequality constraint violation: %g\n', max(ineq_violation));
end
else
ineq_violation = zeros(niq, 1);
end
if any( [ lb_violation;
ub_violation;
eq_violation;
ineq_violation ] > err_tol)
err_cnt = 0;
if any( lb_violation > err_tol )
err_cnt = err_cnt + 1;
errs{err_cnt} = ...
sprintf('variable lower bound violated by %g', ...
max(lb_violation));
end
if any( ub_violation > err_tol )
err_cnt = err_cnt + 1;
errs{err_cnt} = ...
sprintf('variable upper bound violated by %g', ...
max(ub_violation));
end
if any( eq_violation > err_tol )
err_cnt = err_cnt + 1;
errs{err_cnt} = ...
sprintf('equality constraint violated by %g', ...
max(eq_violation));
end
if any( ineq_violation > err_tol )
err_cnt = err_cnt + 1;
errs{err_cnt} = ...
sprintf('inequality constraint violated by %g', ...
max(ineq_violation));
end
if verbose > 0
fprintf('\nWARNING: This version of BPMPD_MEX has a bug which caused it to return\n');
fprintf( ' an incorrect (infeasible) solution for this particular problem.\n');
end
for err_idx = 1:err_cnt
fprintf(' %s\n', errs{err_idx});
end
if bpmpd_bug_fatal
error('%s\n%s', ...
'To avoid this bug in BPMPD_MEX you will need', ...
'to use a different QP solver for this problem.');
end
eflag = -99;
output.message = [output.message '\nBPMPD bug: returned infeasible solution'];
end
end