选用牛顿-拉佛森方法，利用matlab 软件计算基于PQ 节点情况下的潮流计算。
一．所用公式
j =n
⎧
⎪∆P i =P i -∑[e i (G ij e j -B ij f j ) +f i (G ij f j +B ij e j )]
j =1
⎪
j =n ⎪⎪
⎨∆Q i =Q i -∑[f i (G ij e j -B ij f j ) -e i (G ij f j +B ij e j )]
j =1⎪
⎪∆U 2=U 2-(e 2+f 2)
i i i
⎪i ⎪⎩
i ≠j
∂P i ⎧H =⎪ij ∂f =-B ij e i +G ij f i
i
⎪⎪∂P i N =⎪ij ∂e =G ij e i +B ij f i
i ⎪
⎪∂Q i J ==-G ij e i -B ij f i =-N ij ⎪ij
∂f i
⎪
⎨∂Q ⎪L ij =i =-B ij e i +G ij f i =H ij
∂e i ⎪⎪2
∂U i ⎪R ==0ij
⎪∂f i ⎪
∂U i 2⎪
⎪S ij =∂e =0
i ⎩
i =j
∂P i ⎧H =⎪ii ∂f =-B ij e i +G ij f i +b ii
i
⎪⎪∂P i
⎪N ii =∂e =G ij e i +B ij f i +a ii
i ⎪
⎪∂Q i
=-G ij e i -B ij f i +a ii ⎪J ii =∂f i
⎪⎨
⎪L ii =∂Q i =-B ij e i +G ij f i -b ii
∂e i ⎪⎪2
∂U i ⎪R ==2f i ii
⎪∂f i ⎪
∂U i 2⎪
⎪S ii =∂e =2e i
i ⎩
其中
j =n
⎧
⎪a ii =G ii e i -B ii f i +∑(G ij e j -B ij f j )
j =1
⎪j ≠i ⎪⎨j =n
⎪b =G f +B e +(G f +B e )
∑ii i ii i ij j ij j
⎪ii
j =1
⎪j ≠i ⎩
二、程序流程图
牛顿-拉佛森例题中的简单模型系统
一、 系统单线图
图1 简单模型系统
二、系统参数
节点矩阵
%(bus#)(volt)(ang)(p)(q)(bus type) bus=[
1 1.00 0.00 -1.60 -0.80 1; 2 1.00 0.00 -2.00 -1.00 1; 3 1.00 0.00 -3.70 -1.30 1; 4 1.05 0.00 5.00 0.00 2; 5 1.05 0.00 0.00 0.00 3];
线路矩阵
%b#1 b#2 (R)(X)(G)(B)(K) line=[
1 2 0.04 0.25 0 0.25 0; 1 3 0.10 0.35 0 0.00 0; 2 3 0.08 0.30 0 0.25 0; 5 3 0.00 0.03 0 0.00 1.05; 4 2 0.00 0.015 0 0.00 1.05] ;
三、计算结果:
牛顿－拉夫逊法潮流计算结果 节点计算结果：
n 节点 节点电压 节点相角（角度） 节点注入功率 1 0.862150 -4.778511 -1.600000 + j -0.800000 2 1.077916 17.853530 -2.000000 + j -1.000000 3 1.036411 -4.281930 -3.700000 + j -1.300000 4 1.050000 21.843319 5.000000 + j 1.813084 5 1.050000 0.000000 2.579427 + j 2.299402 n 线路计算结果：
n 节点I 节点J 线路功率S(I,J) 线路功率S(J,I) dS(I,J)
1 2 -1.466181 + j -0.409076 1.584546 + j 0.672556 0.263480
1 3 -0.133819 + j -0.390924 0.156788 + j 0.471315 0.080391
线路损耗0.118365 + j 0.022969 + j
2 3 1.415454 + j -0.244333 -1.277360 + j 0.203170 0.138093 + j -0.041163 5 3 2.579427 + j 2.299402 -2.579427 + j -1.974485 0.000000 + j 0.324917
4 2 5.000000 + j 1.813084 -5.000000 + j -1.428223 0.000000 + j 0.384861
导纳矩阵： Y =
[ 1.3787 - 6.2917i -0.6240 + 3.9002i -0.7547 + 2.6415i 0 0 -0.6240 + 3.9002i 1.4539 -66.9808i -0.8299 + 3.1120i 0 +63.4921i 0 -0.7547 + 2.6415i -0.8299 + 3.1120i 1.5846 -35.7379i 0 0 +31.7460i 0 0 +63.4921i 0 0 -66.6667i 0 0 0 0 +31.7460i 0 0 -33.3333i ]
图2. Matlab运行结果
结果分析：此程序的运行结果和试验程序给出的结果是一致的。说明程序无误，但在精确度上有微小差异，这主要是和导纳矩阵的精确度以及显示精度有关。
心得：本程序分模块进行，先是排序，再是求导纳阵，然后求雅阁比，再进行迭代运算，程序本身很简洁明了，运行的时候只需要在matlab 里输入main 就行了，然后打开BUS 和line 所在的.m 文件，结果就会自动存在result 文件中了，通过编写牛顿拉夫逊法matlab 潮流计算程序复习了潮流计算的知识，也实现了计算机算法
附录：
实验源程序： Main 函数： clear
[dfile,pathname]=uigetfile('*.m','Select Data File'); if pathname == 0
error(' you must select a valid data file') else
lfile =length(dfile); % strip off .m
eval(dfile(1:lfile-2)); end;
global bus; global line;
[nb,mb]=size(bus);
[nl,ml]=size(line); % 计算bus 和line 矩阵的行数和列数 [bus,line,nPQ,nPV,nodenum] = Num(bus,line); % 对节点重新排序的子程序 Y = y(bus,line) % 计算节点导纳矩阵的子程序 myf = fopen('Result.m','w');
fprintf(myf,'计算结果');
fclose(myf); % 在当前目录下生成“Result.m ”文件，写入节点导纳矩阵 format long
EPS = 1.0e-10; % 设定误差精度
for t = 1:100 % 开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出
[dP,dQ] = dPQ(Y,bus,nPQ,nPV); % 计算功率偏差dP 和dQ 的子程序 J = Jac(bus,Y,nPQ); % 计算雅克比矩阵的子程序 UD = zeros(nPQ,nPQ); for i = 1:nPQ
UD(i,i) = bus(i,2); % 生成电压对角矩阵 end end
dAngU = J\[dP;dQ];
dAng = dAngU(1:nb-1,1); % 计算相角修正量 dU = UD*(dAngU(nb:nb+nPQ-1,1)); % 计算电压修正量 bus(1:nPQ,2) = bus(1:nPQ,2) - dU; % 修正电压 bus(1:nb-1,3) = bus(1:nb-1,3) - dAng; % 修正相角 if (max(abs(dU))
end % 判断是否满足精度误差，如满足则跳出，否则返回继续迭代 end
bus = PQ(bus,Y,nPQ,nPV); % 计算每个节点的有功和无功注入的子程序
[bus,line] = ReNum(bus,line,nodenum); % 对节点恢复编号的子程序
YtYm = YtYm_(line); % 计算线路的等效Yt 和Ym 的子程序，以计算线路潮流
bus_res = bus_res_(bus); % 计算节点数据结果的子程序 S_res = S_res_(bus,line,YtYm); % 计算线路潮流的子程序 myf = fopen('Result.m','a');
fprintf(myf,'牛顿－拉夫逊法潮流计算结果 节点计算结果：n 节点 节点电压 节点相角（角度） 节点注入功率\n'); for i = 1:nb
fprintf(myf,'%2.0f ',bus_res(i,1)); fprintf(myf,'%10.6f ',bus_res(i,2)); fprintf(myf,'%10.6f ',bus_res(i,3));
fprintf(myf,'%10.6f + j %10.6f\n',real(bus_res(i,4)),imag(bus_res(i,4))); end
fprintf(myf,'n线路计算结果：n 节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)\n'); for i = 1:nl
fprintf(myf,'%2.0f ',S_res(i,1)); fprintf(myf,'%2.0f ',S_res(i,2));
fprintf(myf,'%10.6f + j %10.6f ',real(S_res(i,3)),imag(S_res(i,3))); fprintf(myf,'%10.6f + j %10.6f ',real(S_res(i,4)),imag(S_res(i,4))); fprintf(myf,'%10.6f + j%10.6f\n',real(S_res(i,5)),imag(S_res(i,5))); end
fclose(myf); % 迭代结束后继续在“Result.m ”写入节点计算结果和线路计算结果 程序结束
"Num.m" 作用为对节点重排序，并修改相应的线路数据
function [bus,line,nPQ,nPV,nodenum] = Num(bus,line) [nb,mb]=size(bus); [nl,ml]=size(line);
nSW = 0; % number of swing bus counter nPV = 0; % number of PV bus counter nPQ = 0; % number of PQ bus counter for i = 1:nb, % nb为总节点数 type= bus(i,6); if type == 3,
nSW = nSW + 1; % increment swing bus counter SW(nSW,:)=bus(i,:); elseif type == 2,
nPV = nPV +1; % increment PV bus counter PV(nPV,:)=bus(i,:); else
nPQ = nPQ + 1; % increment PQ bus counter PQ(nPQ,:)=bus(i,:); end
end;
bus=[PQ;PV;SW]; newbus=[1:nb]';
nodenum=[newbus bus(:,1)]; bus(:,1)=newbus; for i=1:nl for j=1:2 for k=1:nb
if line(i,j)==nodenum(k,2) line(i,j)=nodenum(k,1); break end end end end
"y.m" 作用为计算节点导纳矩阵
function Y = y(bus,line) [nb,mb]=size(bus); [nl,ml]=size(line); Y=zeros(nb,nb); for k=1:nl
I=line(k,1); %读入线路参数 J=line(k,2);
Zt=line(k,3)+j*line(k,4); Yt=1/Zt;
Ym=line(k,5)+j*line(k,6); K=line(k,7);
if (K==0)&(J~=0) % 普通线路: K=0; Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt; Y(J,I)=Y(I,J); end
if (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0; Y(I,I)=Y(I,I)+Ym; end
if K>0 % 变压器线路: Zt和Ym 为折算到i 侧的值,K 在j 侧 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt/K/K; Y(I,J)=Y(I,J)-Yt/K; Y(J,I)=Y(I,J); end
if K
Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+K*K*Yt; Y(I,J)=Y(I,J)+K*Yt; Y(J,I)=Y(I,J); end End
"dPQ.m" 作用为计算功率偏差
function [dP,dQ] =dPQ(Y,bus,nPQ,nPV) % nPQ、nPV 为相应节点个数 n = nPQ + nPV +1; % 总节点个数 dP = bus(1:n-1,4);
dQ = bus(1:nPQ,5); % 对dP 和dQ 赋初值 PV节点不需计算dQ 平衡节点不参与计算 for i = 1:n-1 for j = 1:n dP(i,1) = dP(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3))); if i
end % 利用循环计算求取dP 和dQ
"Jac.m" 作用为计算雅克比矩阵
function J = Jac(bus,Y,nPQ) [nb,mb]=size(bus); H = zeros(nb-1,nb-1); N = zeros(nb-1,nPQ); K = zeros(nPQ,nb-1);
L = zeros(nPQ,nPQ); % 将雅克比矩阵分块，即：J = [H N; K L]，并初始化 Qi = zeros(nb-1,1); Pi = zeros(nb-1,1); for i = 1:nb-1
for j = 1:nb
Pi(i,1)=Pi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
Qi(i,1)=Qi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3))); end
end % 初始化并计算Qi 和Pi
for i = 1:nb-1
for j = 1:nb-1
if i~=j
H(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
else
H(i,j)=Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);
end % 分别计算H 矩阵的对角及非对角元素
if j
if i~=j
N(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
else
N(i,j)=-Pi(i,1)-real(Y(i,j))*((bus(i,2))^2);
end
end % 分别计算N 矩阵的对角及非对角元素
if i
if i~=j
K(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
else
K(i,j)=-Pi(i,1)+real(Y(i,j))*((bus(i,2))^2);
end % 分别计算K 矩阵的对角及非对角元素 if j
if i~=j
L(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
else
L(i,j)=-Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);
end
end % 分别计算L 矩阵的对角及非对角元素
end
end
end
J = [H N; K L]; % 生成雅克比矩阵
"PQ.m" 作用为计算每个节点的功率注入
function bus = PQ(bus,Y,nPQ,nPV)
n = nPQ+nPV+1; % n为总节点数
for i = nPQ+1:n-1
for j = 1:n
bus(i,5)=bus(i,5)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
end
end % 利用公式计算PV 节点的无功注入
for j =1:n
bus(n,4)=bus(n,4)+bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+imag(Y(n,j))*sin(bus(n,3)-bus(j,3)));
bus(n,5)=bus(n,5)+bus(n,2)*bus(j,2)*(real(Y(n,j))*sin(bus(n,3)-bus(j,3))-imag(Y(n,j))*cos(bus(n,3)-bus(j,3)));
end % 计算平衡节点的无功及有功注入
"ReNum.m" 作用为对节点和线路数据恢复编号
function [bus,line] = ReNum(bus,line,nodenum)
[nb,mb]=size(bus);
[nl,ml]=size(line);
bus_temp = zeros(nb,mb); % bus_temp矩阵用于临时存放bus 矩阵的数据 k = 1;
for i = 1 :nb
for j = 1 : nb
if bus(j,1) == k
bus_temp(k,:) = bus(j,:);
k = k + 1;
end
end
end % 利用bus 矩阵的首列编号重新对bus 矩阵排序并存入bus_temp矩阵中
bus = bus_temp; % 重新赋值回bus ，完成bus 矩阵的重新编号
for i=1:nl
for j=1:2
for k=1:nb
if line(i,j)==nodenum(k,1)
line(i,j)=nodenum(k,2);
break
end
end
end
end % 恢复line 的编号
"YtYm_.m" 作用为计算线路的等效Yt 和Ym ，以计算线路潮流
function YtYm = YtYm_(line)
[nl,ml]=size(line);
YtYm = zeros(nl,5); % 对YtYm 矩阵赋初值0
YtYm(:,1:2) = line(:,1:2); % 矩阵前两列为线路两段节点编号，后三列分别为线路等效Yt ，i 侧的等效Ym ，j 侧的等效Ym
j = sqrt(-1);
for k=1:nl
I=line(k,1);
J=line(k,2);
Zt=line(k,3)+j*line(k,4);
if Zt~=0
Yt=1/Zt;
end
Ym=line(k,5)+j*line(k,6);
K=line(k,7);
if (K==0)&(J~=0) % 普通线路: K=0
YtYm(k,3) = Yt;
YtYm(k,4) = Ym;
YtYm(k,5) = Ym;
end
if (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0
YtYm(k,4) = Ym;
end
if K>0 % 变压器线路: Zt和Ym 为折算到i 侧的值,K 在j 侧 YtYm(k,3) = Yt/K;
YtYm(k,4) = Ym+Yt*(K-1)/K;
YtYm(k,5) = Yt*(1-K)/K/K;
end
if K
YtYm(k,4) = Ym+Yt*(1+K);
YtYm(k,5) = Yt*(K^2+K);
end
end
"bus_res_.m" 计算并返回节点数据结果
function bus_res = bus_res_(bus)
[nb,mb]=size(bus);
bus_res = zeros(nb,4); % bus_res矩阵储存着节点计算结果 bus_res(:,1:2) = bus(:,1:2);
bus_res(:,3) = bus(:,3) *180 / pi; % 相角采用角度制
bus_res(:,4) = bus(:,4) + (sqrt(-1))*bus(:,5); % 注入功率
"S_res_.m" 计算并返回线路潮流
function S_res = S_res_(bus,line,YtYm)
[nl,ml]=size(line);
S_res = zeros(nl,5); % S_res矩阵储存着线路潮流计算结果
S_res(:,1:2) = line(:,1:2); % 前两列为节点编号
for k=1:nl
I = S_res(k,1);
J = S_res(k,2);
if (J~=0)&(I~=0)
S_res(k,3)=bus(I,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4)))-bus(I,2)*bus(J,2)*(cos(bus(I,3))+j*sin(bus(I,3)))*(conj(cos(bus(J,3))+j*sin(bus(J,3))))*conj(YtYm(k,3));
S_res(k,4)=bus(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))-bus(I,2)*bus(J,2)*(conj(cos(bus(I,3))+j*sin(bus(I,3))))*(cos(bus(J,3))+j*sin(bus(J,3)))*conj(YtYm(k,3));
S_res(k,5)=S_res(k,3) + S_res(k,4); % 利用公式计算非接地支路的潮流 else if(J==0 )
S_res(k,3)=bus(I,2)^2*conj(YtYm(k,4));
S_res(k,5)=S_res(k,3)+S_res(k,4);
else
S_res(k,4)=bus(J,2)^2*conj(YtYm(k,5));
S_res(k,5)=S_res(k,3)+S_res(k,4); % 利用公式计算接地支路的潮流 end
end
end
end
选用牛顿-拉佛森方法，利用matlab 软件计算基于PQ 节点情况下的潮流计算。
一．所用公式
j =n
⎧
⎪∆P i =P i -∑[e i (G ij e j -B ij f j ) +f i (G ij f j +B ij e j )]
j =1
⎪
j =n ⎪⎪
⎨∆Q i =Q i -∑[f i (G ij e j -B ij f j ) -e i (G ij f j +B ij e j )]
j =1⎪
⎪∆U 2=U 2-(e 2+f 2)
i i i
⎪i ⎪⎩
i ≠j
∂P i ⎧H =⎪ij ∂f =-B ij e i +G ij f i
i
⎪⎪∂P i N =⎪ij ∂e =G ij e i +B ij f i
i ⎪
⎪∂Q i J ==-G ij e i -B ij f i =-N ij ⎪ij
∂f i
⎪
⎨∂Q ⎪L ij =i =-B ij e i +G ij f i =H ij
∂e i ⎪⎪2
∂U i ⎪R ==0ij
⎪∂f i ⎪
∂U i 2⎪
⎪S ij =∂e =0
i ⎩
i =j
∂P i ⎧H =⎪ii ∂f =-B ij e i +G ij f i +b ii
i
⎪⎪∂P i
⎪N ii =∂e =G ij e i +B ij f i +a ii
i ⎪
⎪∂Q i
=-G ij e i -B ij f i +a ii ⎪J ii =∂f i
⎪⎨
⎪L ii =∂Q i =-B ij e i +G ij f i -b ii
∂e i ⎪⎪2
∂U i ⎪R ==2f i ii
⎪∂f i ⎪
∂U i 2⎪
⎪S ii =∂e =2e i
i ⎩
其中
j =n
⎧
⎪a ii =G ii e i -B ii f i +∑(G ij e j -B ij f j )
j =1
⎪j ≠i ⎪⎨j =n
⎪b =G f +B e +(G f +B e )
∑ii i ii i ij j ij j
⎪ii
j =1
⎪j ≠i ⎩
二、程序流程图
牛顿-拉佛森例题中的简单模型系统
一、 系统单线图
图1 简单模型系统
二、系统参数
节点矩阵
%(bus#)(volt)(ang)(p)(q)(bus type) bus=[
1 1.00 0.00 -1.60 -0.80 1; 2 1.00 0.00 -2.00 -1.00 1; 3 1.00 0.00 -3.70 -1.30 1; 4 1.05 0.00 5.00 0.00 2; 5 1.05 0.00 0.00 0.00 3];
线路矩阵
%b#1 b#2 (R)(X)(G)(B)(K) line=[
1 2 0.04 0.25 0 0.25 0; 1 3 0.10 0.35 0 0.00 0; 2 3 0.08 0.30 0 0.25 0; 5 3 0.00 0.03 0 0.00 1.05; 4 2 0.00 0.015 0 0.00 1.05] ;
三、计算结果:
牛顿－拉夫逊法潮流计算结果 节点计算结果：
n 节点 节点电压 节点相角（角度） 节点注入功率 1 0.862150 -4.778511 -1.600000 + j -0.800000 2 1.077916 17.853530 -2.000000 + j -1.000000 3 1.036411 -4.281930 -3.700000 + j -1.300000 4 1.050000 21.843319 5.000000 + j 1.813084 5 1.050000 0.000000 2.579427 + j 2.299402 n 线路计算结果：
n 节点I 节点J 线路功率S(I,J) 线路功率S(J,I) dS(I,J)
1 2 -1.466181 + j -0.409076 1.584546 + j 0.672556 0.263480
1 3 -0.133819 + j -0.390924 0.156788 + j 0.471315 0.080391
线路损耗0.118365 + j 0.022969 + j
2 3 1.415454 + j -0.244333 -1.277360 + j 0.203170 0.138093 + j -0.041163 5 3 2.579427 + j 2.299402 -2.579427 + j -1.974485 0.000000 + j 0.324917
4 2 5.000000 + j 1.813084 -5.000000 + j -1.428223 0.000000 + j 0.384861
导纳矩阵： Y =
[ 1.3787 - 6.2917i -0.6240 + 3.9002i -0.7547 + 2.6415i 0 0 -0.6240 + 3.9002i 1.4539 -66.9808i -0.8299 + 3.1120i 0 +63.4921i 0 -0.7547 + 2.6415i -0.8299 + 3.1120i 1.5846 -35.7379i 0 0 +31.7460i 0 0 +63.4921i 0 0 -66.6667i 0 0 0 0 +31.7460i 0 0 -33.3333i ]
图2. Matlab运行结果
结果分析：此程序的运行结果和试验程序给出的结果是一致的。说明程序无误，但在精确度上有微小差异，这主要是和导纳矩阵的精确度以及显示精度有关。
心得：本程序分模块进行，先是排序，再是求导纳阵，然后求雅阁比，再进行迭代运算，程序本身很简洁明了，运行的时候只需要在matlab 里输入main 就行了，然后打开BUS 和line 所在的.m 文件，结果就会自动存在result 文件中了，通过编写牛顿拉夫逊法matlab 潮流计算程序复习了潮流计算的知识，也实现了计算机算法
附录：
实验源程序： Main 函数： clear
[dfile,pathname]=uigetfile('*.m','Select Data File'); if pathname == 0
error(' you must select a valid data file') else
lfile =length(dfile); % strip off .m
eval(dfile(1:lfile-2)); end;
global bus; global line;
[nb,mb]=size(bus);
[nl,ml]=size(line); % 计算bus 和line 矩阵的行数和列数 [bus,line,nPQ,nPV,nodenum] = Num(bus,line); % 对节点重新排序的子程序 Y = y(bus,line) % 计算节点导纳矩阵的子程序 myf = fopen('Result.m','w');
fprintf(myf,'计算结果');
fclose(myf); % 在当前目录下生成“Result.m ”文件，写入节点导纳矩阵 format long
EPS = 1.0e-10; % 设定误差精度
for t = 1:100 % 开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出
[dP,dQ] = dPQ(Y,bus,nPQ,nPV); % 计算功率偏差dP 和dQ 的子程序 J = Jac(bus,Y,nPQ); % 计算雅克比矩阵的子程序 UD = zeros(nPQ,nPQ); for i = 1:nPQ
UD(i,i) = bus(i,2); % 生成电压对角矩阵 end end
dAngU = J\[dP;dQ];
dAng = dAngU(1:nb-1,1); % 计算相角修正量 dU = UD*(dAngU(nb:nb+nPQ-1,1)); % 计算电压修正量 bus(1:nPQ,2) = bus(1:nPQ,2) - dU; % 修正电压 bus(1:nb-1,3) = bus(1:nb-1,3) - dAng; % 修正相角 if (max(abs(dU))
end % 判断是否满足精度误差，如满足则跳出，否则返回继续迭代 end
bus = PQ(bus,Y,nPQ,nPV); % 计算每个节点的有功和无功注入的子程序
[bus,line] = ReNum(bus,line,nodenum); % 对节点恢复编号的子程序
YtYm = YtYm_(line); % 计算线路的等效Yt 和Ym 的子程序，以计算线路潮流
bus_res = bus_res_(bus); % 计算节点数据结果的子程序 S_res = S_res_(bus,line,YtYm); % 计算线路潮流的子程序 myf = fopen('Result.m','a');
fprintf(myf,'牛顿－拉夫逊法潮流计算结果 节点计算结果：n 节点 节点电压 节点相角（角度） 节点注入功率\n'); for i = 1:nb
fprintf(myf,'%2.0f ',bus_res(i,1)); fprintf(myf,'%10.6f ',bus_res(i,2)); fprintf(myf,'%10.6f ',bus_res(i,3));
fprintf(myf,'%10.6f + j %10.6f\n',real(bus_res(i,4)),imag(bus_res(i,4))); end
fprintf(myf,'n线路计算结果：n 节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)\n'); for i = 1:nl
fprintf(myf,'%2.0f ',S_res(i,1)); fprintf(myf,'%2.0f ',S_res(i,2));
fprintf(myf,'%10.6f + j %10.6f ',real(S_res(i,3)),imag(S_res(i,3))); fprintf(myf,'%10.6f + j %10.6f ',real(S_res(i,4)),imag(S_res(i,4))); fprintf(myf,'%10.6f + j%10.6f\n',real(S_res(i,5)),imag(S_res(i,5))); end
fclose(myf); % 迭代结束后继续在“Result.m ”写入节点计算结果和线路计算结果 程序结束
"Num.m" 作用为对节点重排序，并修改相应的线路数据
function [bus,line,nPQ,nPV,nodenum] = Num(bus,line) [nb,mb]=size(bus); [nl,ml]=size(line);
nSW = 0; % number of swing bus counter nPV = 0; % number of PV bus counter nPQ = 0; % number of PQ bus counter for i = 1:nb, % nb为总节点数 type= bus(i,6); if type == 3,
nSW = nSW + 1; % increment swing bus counter SW(nSW,:)=bus(i,:); elseif type == 2,
nPV = nPV +1; % increment PV bus counter PV(nPV,:)=bus(i,:); else
nPQ = nPQ + 1; % increment PQ bus counter PQ(nPQ,:)=bus(i,:); end
end;
bus=[PQ;PV;SW]; newbus=[1:nb]';
nodenum=[newbus bus(:,1)]; bus(:,1)=newbus; for i=1:nl for j=1:2 for k=1:nb
if line(i,j)==nodenum(k,2) line(i,j)=nodenum(k,1); break end end end end
"y.m" 作用为计算节点导纳矩阵
function Y = y(bus,line) [nb,mb]=size(bus); [nl,ml]=size(line); Y=zeros(nb,nb); for k=1:nl
I=line(k,1); %读入线路参数 J=line(k,2);
Zt=line(k,3)+j*line(k,4); Yt=1/Zt;
Ym=line(k,5)+j*line(k,6); K=line(k,7);
if (K==0)&(J~=0) % 普通线路: K=0; Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt; Y(J,I)=Y(I,J); end
if (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0; Y(I,I)=Y(I,I)+Ym; end
if K>0 % 变压器线路: Zt和Ym 为折算到i 侧的值,K 在j 侧 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt/K/K; Y(I,J)=Y(I,J)-Yt/K; Y(J,I)=Y(I,J); end
if K
Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+K*K*Yt; Y(I,J)=Y(I,J)+K*Yt; Y(J,I)=Y(I,J); end End
"dPQ.m" 作用为计算功率偏差
function [dP,dQ] =dPQ(Y,bus,nPQ,nPV) % nPQ、nPV 为相应节点个数 n = nPQ + nPV +1; % 总节点个数 dP = bus(1:n-1,4);
dQ = bus(1:nPQ,5); % 对dP 和dQ 赋初值 PV节点不需计算dQ 平衡节点不参与计算 for i = 1:n-1 for j = 1:n dP(i,1) = dP(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3))); if i
end % 利用循环计算求取dP 和dQ
"Jac.m" 作用为计算雅克比矩阵
function J = Jac(bus,Y,nPQ) [nb,mb]=size(bus); H = zeros(nb-1,nb-1); N = zeros(nb-1,nPQ); K = zeros(nPQ,nb-1);
L = zeros(nPQ,nPQ); % 将雅克比矩阵分块，即：J = [H N; K L]，并初始化 Qi = zeros(nb-1,1); Pi = zeros(nb-1,1); for i = 1:nb-1
for j = 1:nb
Pi(i,1)=Pi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
Qi(i,1)=Qi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3))); end
end % 初始化并计算Qi 和Pi
for i = 1:nb-1
for j = 1:nb-1
if i~=j
H(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
else
H(i,j)=Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);
end % 分别计算H 矩阵的对角及非对角元素
if j
if i~=j
N(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
else
N(i,j)=-Pi(i,1)-real(Y(i,j))*((bus(i,2))^2);
end
end % 分别计算N 矩阵的对角及非对角元素
if i
if i~=j
K(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
else
K(i,j)=-Pi(i,1)+real(Y(i,j))*((bus(i,2))^2);
end % 分别计算K 矩阵的对角及非对角元素 if j
if i~=j
L(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
else
L(i,j)=-Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);
end
end % 分别计算L 矩阵的对角及非对角元素
end
end
end
J = [H N; K L]; % 生成雅克比矩阵
"PQ.m" 作用为计算每个节点的功率注入
function bus = PQ(bus,Y,nPQ,nPV)
n = nPQ+nPV+1; % n为总节点数
for i = nPQ+1:n-1
for j = 1:n
bus(i,5)=bus(i,5)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
end
end % 利用公式计算PV 节点的无功注入
for j =1:n
bus(n,4)=bus(n,4)+bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+imag(Y(n,j))*sin(bus(n,3)-bus(j,3)));
bus(n,5)=bus(n,5)+bus(n,2)*bus(j,2)*(real(Y(n,j))*sin(bus(n,3)-bus(j,3))-imag(Y(n,j))*cos(bus(n,3)-bus(j,3)));
end % 计算平衡节点的无功及有功注入
"ReNum.m" 作用为对节点和线路数据恢复编号
function [bus,line] = ReNum(bus,line,nodenum)
[nb,mb]=size(bus);
[nl,ml]=size(line);
bus_temp = zeros(nb,mb); % bus_temp矩阵用于临时存放bus 矩阵的数据 k = 1;
for i = 1 :nb
for j = 1 : nb
if bus(j,1) == k
bus_temp(k,:) = bus(j,:);
k = k + 1;
end
end
end % 利用bus 矩阵的首列编号重新对bus 矩阵排序并存入bus_temp矩阵中
bus = bus_temp; % 重新赋值回bus ，完成bus 矩阵的重新编号
for i=1:nl
for j=1:2
for k=1:nb
if line(i,j)==nodenum(k,1)
line(i,j)=nodenum(k,2);
break
end
end
end
end % 恢复line 的编号
"YtYm_.m" 作用为计算线路的等效Yt 和Ym ，以计算线路潮流
function YtYm = YtYm_(line)
[nl,ml]=size(line);
YtYm = zeros(nl,5); % 对YtYm 矩阵赋初值0
YtYm(:,1:2) = line(:,1:2); % 矩阵前两列为线路两段节点编号，后三列分别为线路等效Yt ，i 侧的等效Ym ，j 侧的等效Ym
j = sqrt(-1);
for k=1:nl
I=line(k,1);
J=line(k,2);
Zt=line(k,3)+j*line(k,4);
if Zt~=0
Yt=1/Zt;
end
Ym=line(k,5)+j*line(k,6);
K=line(k,7);
if (K==0)&(J~=0) % 普通线路: K=0
YtYm(k,3) = Yt;
YtYm(k,4) = Ym;
YtYm(k,5) = Ym;
end
if (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0
YtYm(k,4) = Ym;
end
if K>0 % 变压器线路: Zt和Ym 为折算到i 侧的值,K 在j 侧 YtYm(k,3) = Yt/K;
YtYm(k,4) = Ym+Yt*(K-1)/K;
YtYm(k,5) = Yt*(1-K)/K/K;
end
if K
YtYm(k,4) = Ym+Yt*(1+K);
YtYm(k,5) = Yt*(K^2+K);
end
end
"bus_res_.m" 计算并返回节点数据结果
function bus_res = bus_res_(bus)
[nb,mb]=size(bus);
bus_res = zeros(nb,4); % bus_res矩阵储存着节点计算结果 bus_res(:,1:2) = bus(:,1:2);
bus_res(:,3) = bus(:,3) *180 / pi; % 相角采用角度制
bus_res(:,4) = bus(:,4) + (sqrt(-1))*bus(:,5); % 注入功率
"S_res_.m" 计算并返回线路潮流
function S_res = S_res_(bus,line,YtYm)
[nl,ml]=size(line);
S_res = zeros(nl,5); % S_res矩阵储存着线路潮流计算结果
S_res(:,1:2) = line(:,1:2); % 前两列为节点编号
for k=1:nl
I = S_res(k,1);
J = S_res(k,2);
if (J~=0)&(I~=0)
S_res(k,3)=bus(I,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4)))-bus(I,2)*bus(J,2)*(cos(bus(I,3))+j*sin(bus(I,3)))*(conj(cos(bus(J,3))+j*sin(bus(J,3))))*conj(YtYm(k,3));
S_res(k,4)=bus(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))-bus(I,2)*bus(J,2)*(conj(cos(bus(I,3))+j*sin(bus(I,3))))*(cos(bus(J,3))+j*sin(bus(J,3)))*conj(YtYm(k,3));
S_res(k,5)=S_res(k,3) + S_res(k,4); % 利用公式计算非接地支路的潮流 else if(J==0 )
S_res(k,3)=bus(I,2)^2*conj(YtYm(k,4));
S_res(k,5)=S_res(k,3)+S_res(k,4);
else
S_res(k,4)=bus(J,2)^2*conj(YtYm(k,5));
S_res(k,5)=S_res(k,3)+S_res(k,4); % 利用公式计算接地支路的潮流 end
end
end
end