lprelax.m

关于电力系统的计算程序!最新的改进法计算潮流程序！用matlab实现的！ 适用于电力行业的相关从业人员！

function [x2, duals_rlx, idx_workc, idx_bindc] = LPrelax(a, f, b, nequs, vlb, vub, idx_workc, mpopt);

%   MATPOWER Version 2.0
%   by Deqiang (David) Gan, PSERC Cornell    12/12/97
%   Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/ for more info.

%% options
alg		= mpopt(11);
verbose	= mpopt(31);

if opf_slvr(alg) == 1
	idx_workc = find(b < 0.001);
end


converged = 0;
while converged == 0


	atemp = a(idx_workc, :);
	btemp = b(idx_workc);



	atemp = full(atemp);

	[x2, duals] = lp(f, atemp, btemp, vlb, vub, [], nequs, -1);
	

	diffs = b - a * x2; 				% diffs should be normalized by what means? under development
	idx_bindc = find(diffs < 1.0e-8);
	


%	if verbose,       % print out something for the developer
%		fprintf('\n                                                   (%g, %g) \', ...
%									length(idx_workc), length(idx_bindc));  
%	end;


	idx_violated = find(diffs < -1.0e-8);

	if isempty(idx_violated) == 1
		converged = 1;
	else
		flag = zeros(length(b), 1);			% set up flag from scratch
		flag(idx_workc) = ones(length(idx_workc), 1);	% enforce historical working constraints


		idx_add = find(diffs < 0.001);
		flag(idx_add) = ones(length(idx_add), 1);	% enforce violating constraints



		flag(1:nequs) = ones(nequs, 1); 		% enforce original equality constraints
		idx_workc_new = find(flag);

		if length(idx_workc) == length(idx_workc_new)	% safeguard step
			if isempty(find(idx_workc - idx_workc_new)) ==1
				converged = 1;
			end
		end
		idx_workc = idx_workc_new;

	end

end

duals_rlx = zeros(length(b), 1);
duals_rlx(idx_workc) = duals(1:length(btemp));


return;