%主程序"PowerFlow_NR.m"
function [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,'\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,'.6f ',bus_res(i,2));
fprintf(myf,'.6f ',bus_res(i,3));
fprintf(myf,'.6f + j .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,'.6f + j .6f ',real(S_res(i,3)),imag(S_res(i,3)));
fprintf(myf,'.6f + j .6f ',real(S_res(i,4)),imag(S_res(i,4)));
fprintf(myf,'.6f + j .6f\n',real(S_res(i,5)),imag(S_res(i,5)));
end
fclose(myf); % 迭代结束后继续在“Result.m”写入节点计算结果和线路计算结果程序结束
bus = [
1    1.00 0.00 -0.30 -0.18 1;
2    1.00 0.00 -0.55 -0.1
3 1;
3    1.10 0.00 0.50 0.00 2;
4    1.0
5 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.1
2 0.50 0.0 0.0192 0;
3    1 0.00 0.3 0.0 0.0 -1.1]
----------------------------------------牛顿－拉夫逊法潮流计算结果-------------------------------------------- 节点计算结果：
节点节点电压节点相角（角度）节点注入功率
1 0.984674906330845 -0.500170385513657 -0.300000000000000 - 0.180000000000000i
2 0.964797665550885 -6.450305258622626 -0.550000000000000 - 0.130000000000000i
3    1.100000000000000    6.732349388989963 0.500000000000000 + 0.093411003244513i
4    1.050000000000000 0 0.367882692523292 + 0.264698252215732i 线路计算结果：
节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)
1    2 0.246244-j0.014651 -0.239990+j0.010627 0.006254- j0.004024
1    3 -0.500001-j0.029264 0.500000+j0.093409 -0.000001+j0.064145
1    4 -0.046244-j0.136088 0.048216+j0.10452
2 0.001972-j0.031566
2    4 -0.310010-j0.140627 0.319666+j0.160176 0.009656+j0.019549
%子程序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); % 读取文件字符串长度
eval(dfile(1:lfile-2)); % 去除后缀，打开文件！注意！新浪博客中"eval"函数会自动加上"_r"后缀，此处的函数名是"eval"而不是"eval_r"，拷贝后请去除"_r"后缀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);fprintf作用
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

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。