牛顿-拉夫逊法潮流计算matlab程序.doc

牛顿-拉夫逊法潮流计算matlab程序主程序 PowerFlow_NR.mfunction bus_res,S_res = PowerFlow_NR_2 % 牛顿-拉夫逊法解潮流方程的主程序 bus,line = OpDF_; % 打开数据文件的子程序，返回bus（节点数据）和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,- by longdinhohe - /longdinhohe -nn);fclose(myf); % 在当前目录下生成“Result.m”文件，写入节点导纳矩阵format longEPS = 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 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)EPS)&(max(abs(dAng)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 K0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧 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); endend 子程序4 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 inPQ+1 dQ(i,1) = dQ(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 endend % 利用循环计算求取dP和dQ 子程序5 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); endend % 初始化并计算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 nPQ+1 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 nPQ+1 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 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 K0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧 YtYm(k,3) = -Yt*K; YtYm(k,4) = Ym+Yt*(1+K); YtYm(k,5) = Yt*(K2+K); endend 子程序9 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); % 注入功率 子程序10 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); % 利用公式计