qp.m

在matpower的基础上改动的、比较好的一个潮流计算程序！

function [xout, lambdaout, howout] = qp(H,f,A,b,VLB,VUB,x0,N,verbosein)
% Functionally equivalent to Optimization Toolbox QP.M, but calls
% BPMPD instead; accepts both full and sparse H and A.
% Type "help qp" before installing for help on the syntax.

%   MATPOWER Version 1.2
%   by Carlos Murillo-Sanchez, PSERC Cornell    12/10/97
%   Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/ for more info.

n = length(f);
m = length(b);

if nargin < 9
	verbosein = 0;
	if nargin < 8
		N = 0;
		if nargin < 7
			x0 = zeros(n,1);	% Until bpmpd features warm start, this is a dummy arg
			if nargin < 6
				VUB = [];
				if nargin < 5
					VLB = [];
				end
			end
		end
	end
end

if ~isempty(H)
	if ~issparse(H)
		H = sparse(H);
	end
end

if ~issparse(A)
	A = sparse(A);
end

e = -ones(m,1);
if N>0
	e(1:N,:) = zeros(N,1);
end

if ~isempty(VLB)
	llist = 1:n;
	lval = VLB;
else
	llist = [];
	lval = [];
end

if ~isempty(VUB)
	ulist = 1:n;
	uval = VUB;
else
	ulist = [];
	uval = [];
end

if verbosein == -1
	prnlev = 0;
else
	prnlev = 1;
end

myopt = bpopt;
%myopt(14)= 1e-1;	% TPIV1  first relative pivot tolerance (desired)
%myopt(20)= 1e-10;	% TOPT1  stop if feasible and rel. dual gap less than this
%myopt(22)= 1e-6;	% TFEAS1 relative primal feasibility tolerance
myopt(23)= 1e-6;	% TFEAS2 relative dual feasibility tolerance
%myopt(29)= 1e-10;	% TRESX  acceptable primal residual
%myopt(30)= 1e-10;	% TRESY  acceptable dual residual
%myopt(38)= 0;		% SMETHOD1 prescaling method

[xout,y,s,w,howout] = bp(H, A, b, f, e, llist, lval, ...
						ulist, uval, myopt, prnlev);

ilt = find(w<=0);
igt = find(w>0);
mulow = zeros(n,1);
muupp = zeros(n,1);
muupp(ilt) = -w(ilt);
mulow(igt) = w(igt);

lambdaout = -y;
if ~isempty(VLB)
	lambdaout = [lambdaout; mulow];
end
if ~isempty(VUB)
	lambdaout = [lambdaout; muupp];
end

% zero out lambdas smaller than a certain tolerance
ii = find(abs(lambdaout)<1e-9);
lambdaout(ii) = zeros(size(ii));

% The next is necessary for proper operation of constr.m
if strcmp(howout, 'infeasible primal')
	lambdaout = zeros(size(lambdaout));
end

return;