前推回代配电网潮流计算附33节点算例

%author:Jerry Ng %time:2019.01.03 %place:Nanchang University %头结点编号要比尾节点编号大 %利用尾节点没有子节点的特点生成前推回代顺序表 %每将一个支路加入前推回代顺序表，就将该支路的尾节点删去

function BackForwardPFC_dfs33（branchdatavoltage）
%author:Jerry Ng
%time:2019.01.03
%place:Nanchang University
%头结点编号要比尾节点编号大
%利用尾节点没有子节点的特点生成前推回代顺序表
%每将一个支路加入前推回代顺序表，就将该支路的尾节点删去

%branchdata
%第一列 头节点i编号
%第二列 尾节点j编号
%第三列 支路i-j的线路阻抗
%第四列 节点j的负荷功率

%voltage（i） 节点i电压
%bus_power（i） 节点i复功率
%impedance（ij） 支路i-j的线路阻抗
%tailbus headbus 前推回代顺序表 
%busnum 节点数量
%current（ij） 支路i-j的支路电流 current上三角矩阵

	clc
	if nargin < 3
		branchdata = [
			1 	2 	0.0922+0.047j;	
			2 	3 	0.4930+0.2511j;	
			3 	4 	0.3660+0.1864j; 	
			4 	5 	0.3811+0.1941j;	
			5 	6 	0.8190+0.7070j; 	
			6 	7	0.1872+0.6188j; 
			7 	8 	0.7114+0.2351j; 
			8 	9 	1.0300+0.7400j; 	
			9 	10 	1.0440+0.7400j; 	
			10 	11 	0.1966+0.0650j; 	
			11 	12 	0.3744+0.1238j; 
			12 	13 	1.4680+1.1550j; 
			13 	14 	0.5416+0.7129j;
			14 	15 	0.5910+0.5260j; 
			15 	16 	0.7463+0.5450j; 
			16 	17 	1.2890+1.7210j; 
			17 	18 	0.3720+0.5740j;
			2 	19 	0.1640+0.1565j;
			19 	10 	1.5042+1.3554j;
			20 	21 	0.4095+0.4784j;
			21 	22 	0.7089+0.9373j;
			3 	23 	0.4512+0.3083j;
			23 	24 	0.8980+0.7091j;
			24 	25 	0.8960+0.7011j;
			6 	26 	0.2030+0.1034j;
			26 	27 	0.2842+0.1447j;
			27 	28 	1.0590+0.9337j;
			28 	29 	0.8042+0.7006j;
			29 	30	0.5075+0.2585j; 
			30 	31 	0.9744+0.9630j;
			31 	32 	0.3105+0.3619j;
			32 	33 	0.3410+0.5362j];%线路参数
			bus_power = [
						0100+60j90+40j120+80j60+30j60+20j200+100j200+100j...
						60+20j60+20j45+30j60+35j60+35j120+80j60+10j...
						60+20j60+20j90+40j90+40j90+40j90+40j90+40j...
						90+50j420+200j420+200j60+25j60+25j60+20j120+70j...
						200+600j150+70j210+100j60+40j];
		voltage = [12.66];
		%基准电压
	end
	bus_power = bus_power/1000;
	branchdata（:3） = 1*real（branchdata（:3））+1j*imag（branchdata（:3））;
	[voltageimpedancetailbusheadbusbusnum] = ...
		DataProcess（branchdatavoltage）;
	tailbus 
    headbus
	for n = 1:100		
        fprintf（‘第%d轮迭代结果\n‘n）;	
        fprintf（‘各支路的首段电流 A\n‘）;
		current = Backward（bus_powervoltagetailbusheadbusbusnum）;
		
		voltage_last = voltage;
		fprintf（‘各节点电压 kV\n‘）;
		voltage = Forward（currentvoltageimpedancetailbusheadbusbusnum）; 
		
		%收敛判据
		if max（abs（voltage-voltage_last）/voltage（1））<1e-9
			break;
		end	
	end
end

function [voltageimpedancetailbusheadbusbusnum] = ...
			DataProcess（branchdatavoltage）
	%数据预处理
    busnum = max（[branchdata（:1）‘branchdata（:2）‘]）;
    if length（voltage）<2
		%当U只给出了电源点的电压时，
		%根据电源点电压高于其他节点电压5%，计算其他节点的初始电压
		voltage（2:busnum） = ones（1busnum-1）*（voltage（1）/1.05）;
	end
	[branchdatavoltage] = BusRename（branchdatavoltagebusnum）;
	
	impedance = 100001*ones（busnum）;%impedance（ij）表示节点i j之间的阻抗
	for n = 1:busnum-1
		temp1 = branchdata（n1）;
		temp2 = branchdata（n2）;
		%将各支路参数保存到一个矩阵中
		impedance（temp1temp2） = branchdata（n3）;
		impedance（temp2temp1） = impedance（temp1temp2）;
	end
	
	%[tailbusheadbus] = PathDesign（impe