Matlab 牛拉法代码.doc

牛顿-拉夫逊法潮流计算的Matlab实现 (2009-09-26 13:23:26)转载标签： 电力系统 潮流计算 牛顿-拉夫逊 牛顿拉夫逊 matlab newton-raphson写了近一个礼拜的程序，藏着可惜，拿出来给需要的人转载请注明出处 有不少朋友反映程序报错： ? Output argument bus (and maybe others) not assigned during call toE:MATbestOpDF_.mOpDF_.是因为输入数据的格式问题 输入数据中是包含bus,line的，并非只有两个矩阵，如：bus = 1 1.00 0.00 -0.30 -0.18 1; 2 1.00 0.00 -0.55 -0.13 1; 3 1.10 0.00 0.50 0.00 2; 4 1.05 0.00 0.00 0.00 3;line = 1 2 0.10 0.40 0.0 0.01528 0; 4 2 0.08 0.40 0.0 0.01413 0; 1 4 0.12 0.50 0.0 0.0192 0; 3 1 0.00 0.3 0.0 0.0 -1.1;所有蓝色文字都为数据文件中的内容 -分-割-线- 【注意】 一、此程序只适用于求解节点电压以极坐标形式表示的潮流方程，没有考虑节点优化编号 二、程序在Matlab 6.5上测试通过，应该适用于目前其所有后续版本 三、程序主要通过文件方式输入输出（同时程序也返回结果向量），对输入文件格式有严格要求，具体如下： 1.输入文件可以直接在Matlab中新建m文件编写，也内容可以以文本方式编写，但最后必须改后缀名成为“.m”文件。文件名第一个字符必须是字母，后面可以跟字母、数字和下划线的任何组合，但不能和已有文件和函数冲突，不能含中文。 2.输入文件内容格式： 内容中包括bus（节点数据）和line（线路数据）两个矩阵 其中bus的格式为： bus = 节点编号 节点电压 节点相角（弧度制） 有功注入 无功注入 节点类型; line的格式为： line = 节点i 节点j 线路电阻 电抗 电导 电纳 变压器变比（普通线路为零）; 四、程序分主程序和各子程序，需全部放入Matlab当前目录下调用，计算开始后会在当前目录下生成文件Result.m保存运算结果。 -分-割-线- 以下为程序代码（红色部分，“%”后为注释）： 主程序 PowerFlow_NR.mfunction bus_res,S_res = PowerFlow_NR_2 % 牛顿-拉夫逊法解潮流方程的主程序% by longdinhohe % /longdinhohe % 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(Yt