Matlab 牛拉法代码.doc

牛顿-拉夫逊法潮流计算的Matlab实现 (2009-09-26 13:23:26)转载
标签： 电力系统 潮流计算 牛顿-拉夫逊 牛顿拉夫逊 matlab newton-raphson
写了近一个礼拜的程序，藏着可惜，拿出来给需要的人
转载请注明出处

有不少朋友反映程序报错：
??? Output argument "bus" (and maybe others) not assigned during call to
"E:\MAT\best\OpDF_.m>OpDF_".
是因为输入数据的格式问题 输入数据中是包含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.m"
function [bus_res,S_res] = PowerFlow_NR_2 % 牛顿-拉夫逊法解潮流方程的主程序
%%%%%%%%%%%%%%%%%%%%% by 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 --------------- ---------------\n\n');
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
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))<EPS)
break
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 节点计算结果：\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”写入节点计算结果和线路计算结果 程序结束

子程序1 "OpDF_.m" 作用为打开数据文件
function [bus,line] = OpDF_
[dfile,pathname]=uigetfile('*.m','Select Data File'); % 数据文件类型为m文件，窗口标题为“Select Data File”
if pathname == 0
error(' you must select a valid data file') % 如果没有选择有效文件，则出现错误提示
else
lfile =length(dfile); % 去除文件后缀
(dfile(1:lfile-2)); % 打开文件
end

子程序2 "Num.m" 作用为对节点重排序，并修改相应的线路数据
function [bus,line,nPQ,nPV,nodenum] = Num_(bus,line)
[nb,mb]=size(bus);
[nl,ml]=size(line);
nSW = 0; % nSW为平衡节点个数
nPV = 0; % nPV为PV节点个数
nPQ = 0; % nPQ为PQ节点个数
for i = 1:nb, % nb为总节点数
type= bus(i,6);
if type == 3,
nSW = nSW + 1;
SW(nSW,:)=bus(i,:);% 计算并储存平衡节点
elseif type == 2,
nPV = nPV +1;
PV(nPV,:)=bus(i,:);% 计算并储存PV节点
else
nPQ = nPQ + 1;
PQ(nPQ,:)=bus(i,:);% 计算并储存PQ节点
end
end
bus=[PQ;PV;SW]; % 对bus矩阵按PQ、PV、平衡节点的顺序重新排序
nodenum=[[1:nb]' bus(:,1)];% 生成新旧节点对照表
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 % 按排序以后的节点顺序对line矩阵重新编号

子程序3 "y_.m" 作用为计算节点导纳矩阵
function Y = y_(bus,line)
[nb,mb]=size(bus);
[nl,ml]=size(line);
Y=zeros(nb,nb); % 对导纳矩阵赋初值0
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; % 接地支路不计算Yt
end
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<0 % 变压器线路: 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);
end
end

子程序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 i<nPQ+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
end
end % 利用循环计算求取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)));
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 < 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 < nPQ+1
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]; % 生成雅克比矩阵

子程序6 "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 % 计算平衡节点的无功及有功注入

子程序7 "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的编号

子程序8 "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<0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧
YtYm(k,3) = -Yt*K;
YtYm(k,4) = Ym+Yt*(1+K);
YtYm(k,5) = Yt*(K^2+K);
end
end

子程序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); % 利用公式计算非接地支路的潮流
elseif 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