潮流分析计算

%本程序的功能是用牛顿——拉夫逊法进行潮流计算
% B1矩阵：1、支路首端号；   2、末端号；   3、支路阻抗；   4、线路对地电纳 （或变压器导纳）；
%		  5、支路的变比；   6、支路首端处于K侧为1，1侧为0；
%          7、线路/变压器标识（0/1）变压器参数当支路首端处于K侧标识为1时归算至末端侧，0归算至首端侧
% B2矩阵：1、该节点发电机功率；      2、该节点负荷功率；
%         3、PQ节点电压初始值，或平衡节点及PV节点电压的给定值
%         4、节点所接无功补偿并联电容（感）的电纳
%         5、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；3为PV节点；
clear;
isb=1;      %input（‘请输入平衡母线节点号：isb=‘）;
pr=0.00001;    %input（‘请输入误差精度：pr=‘）;
%---------------------------------------------------
n=6;%input（‘请输入节点数：n=‘）;
nl=3;%input（‘请输入支路数：nl=‘）;
B1=[1 2 01i 0 ]
B2=[0     0        1.20      0   1;
   0   2.1+1i      1.00      0   2;
   0   1.8+0.40i   1.00      0   2;
   0   1.6+0.8i    1.00      0   2;
   0   3.7+1.3i    1.00      0   2;
   0   -5+0i       1.20      0   3] 
%input（‘请输入各节点参数形成的矩阵： B2=‘）;
%X=[1 0;2 0;3 0;4 0;5 0;6 0]
%-------------------------------------------------------------
%n=4;%input（‘请输入节点数：n=‘）;nl=4;%input（‘请输入支路数：nl=‘）;
%B1=[1 2 4+16i 0 1 0 0;1 3 4+16i 0 1 0 0;2 3 2+8i 0 1 0 0;2 4 1.49+48.02i 0 11/110 0 1] %input（‘请输入由支路参数形成的矩阵： B1=‘）;
%B2=[0 0 115 0 1;0 0 110 0 2;0 20+4i 110 0 2;0 10+6i 10 0 2] %input（‘请输入各节点参数形成的矩阵： B2=‘）;
%-------------------------------------------------------------
Y=zeros（n）;e=zeros（1n）;f=zeros（1n）;V=zeros（1n）;sida=zeros（1n）;S1=zeros（nl）;
% % %-----------求导纳矩阵------------------------
%for i=1:n
   % if X（i2）~=0;
      %  p=X（i1）;
       % Y（pp）=1/X（i2）;
    %end
%end
for i=1:nl								%从1到n1（总支路数）
	if B1（i7）==1        %-----------如果是变压器支路--------
    	if B1（i6）==0						%左节点（首端）处于1侧
       		p=B1（i1）;q=B1（i2）;
    	else                                %左节点（首端）处于K侧
       		p=B1（i2）;q=B1（i1）;
    	end
    	Y（pq）=Y（pq）-1./（B1（i3）*B1（i5））;	%非对角元
    	Y（qp）=Y（pq）;                      %非对角元
   		Y（qq）=Y（qq）+1./（B1（i3）*B1（i5）^2）;	%对角元K侧
    	Y（pp）=Y（pp）+1./B1（i3）+B1（i4）;		%对角元1侧+励磁导纳 
	else        %------------否则为线路支路--------------------
		p=B1（i1）;q=B1（i2）;
		Y（pq）=Y（pq）-1./B1（i3）;	%非对角元
    	Y（qp）=Y（pq）;              %非对角元
   		Y（qq）=Y（qq）+1./B1（i3）+B1（i4）./2.0000;	%对角元j侧+线路电纳的一半
    	Y（pp）=Y（pp）+1./B1（i3）+B1（i4）./2.0000;	%对角元i侧+线路电纳的一半
	end
end
disp（‘导纳矩阵 Y=‘）;disp（Y）;
%-----------给定各节点初始电压及给定各节点注入功率--------------------------
G=real（Y）;B=imag（Y）;			%分解出导纳阵的实部和虚部  
for i=1:n						%给定各节点初始电压的实部和虚部    
    e（i）=real（B2（i3））;
    f（i）=imag（B2（i3））;
    V（i）=abs（B2（i3））;				%PV、平衡节点及PQ节点电压模值 
end
for i=1:n						%给定各节点注入功率    
    S（i）=B2（i1）-B2（i2）; 		%i节点注入功率SG-SL  
    B（ii）=B（ii）+B2（i4）;		%i节点无功补偿量（电纳值）
end
%==================用牛顿-拉夫逊法迭代求解非线性代数方程（功率方程）=======================
P=real（S）;Q=imag（S）;            %分解出各节点注入的有功和无功功率
ICT1=0;IT2=1;N0=2*n;N1=N0+1;a=0; %迭代次数ICT1、a；不满足收敛要求的节点数IT2
while IT2~=0                   % N0=2*n 雅可比矩阵的阶数；N1=N0+1扩展列
      IT2=0;a=a+1;
      JZ=[‘Jacobi矩阵第（‘num2str（a）‘）次消去运算‘];JZ1=[‘Jacobi矩阵第（‘num2str（a）‘）次回代运算‘];JZ0=[‘功率方程第（‘num2str（a）‘）次差值：‘];
    %----------------求取各个节点的功率及功率偏差及PV节点的电