华南理工大学广州学院电气工程学院电力系统分析课程设计报告
X(f,f)=X(f,f)-1/Trans(i).X; X(t,t)=X(t,t)-1/Trans(i).X; X(f,t)=1/Trans(i).X; X(t,f)=1/Trans(i).X; end
G=real(Y); B=imag(Y); end
%生成JP、JQ矩阵
function [JP,JQ]=FormJPQmatrix(Bus,B,B2,busnumber) JP=B2; JQ=B;
for i=1:busnumber
if Bus(i).type==1 %该循环将JP和JQ矩阵对应平衡节点的非对角元素置零,对角元素置1(相当于去掉了对应平衡节点的JP和JQ矩阵的行和列) for k=1:busnumber JQ(i,k)=0.; JQ(k,i)=0.; JP(i,k)=0.;
JP(k,i)=0.; end
JQ(i,i)=1.; JP(i,i)=1.; end
if Bus(i).type==3 %该循环将JQ矩阵对应PV节点的非对角元素置零,对角元素置1(相当于去掉了对应PV节点JQ矩阵的行和列) for k=1:busnumber JQ(i,k)=0.;
JQ(k,i)=0.; end
JQ(i,i)=1.; end end end
%解潮流函数 function
[Bus]=solve_flow_eq(Bus,B,G,busnumber,iJP,iJQ,epsilon,max_no_itearation) for i=1:busnumber
NodeV(i)=Bus(i).V; Nodea(i)=Bus(i).angle;
dQGQL(i)=Bus(i).QG-Bus(i).QL; dPGPL(i)=Bus(i).PG-Bus(i).PL;
17
华南理工大学广州学院电气工程学院电力系统分析课程设计报告
end
NodeV=NodeV'; %节点电压模值列向量初值 Nodea=Nodea'; %节点电压相角列向量初值 dQGQL=dQGQL'; %节点注入无功列向量 dPGPL=dPGPL'; %节点注入无功列向量 pause
maxdP=1.; %最大节点不平衡有功初值 maxdQ=1.; %最大节点不平衡无功初值 epsilon=0.000001; %设置潮流计算精度 noiteration=0; %迭代初始次数为零
while (maxdP>epsilon)&(maxdQ>epsilon) %条件1：(maxdP>epsilon)且条件2：(maxdP>epsilon) %都成立时继续循环
[deltaP,deltaQ,maxdP,maxdQ]=FormdPQvector(Bus,NodeV,Nodea,dQGQL,dPGPL,B,G,busnumber);
%函数FormdPQvector调用
deltaP; %节点不平衡有功向量 deltaQ; %节点不平衡有功向量 maxdP %最大节点不平衡有功 maxdQ %最大节点不平衡无功
da=iJP*deltaP; %解修正方程 dV=iJQ*deltaQ; %解修正方程 Nodea=Nodea+da; %修正节点相角 NodeV=NodeV+dV; %修正节点电压 noiteration=noiteration+1; %迭代次数累计 pause
if noiteration>max_no_itearation %if语句，当迭代次数大于最大次数终止迭代 break end end
for i=1:busnumber
Bus(i).V=NodeV(i); %计算节点电压结果存入循环体 NodeV(i)=NodeV(i)*Bus(i).Base_KV; %转换成有名值
Bus(i).angle=Nodea(i); %计算节点相角（单位弧度）结果存入循环体 Nodea(i)=Nodea(i)*180/pi; %单位弧度转换成单位度 end
%计算支路潮流程序在这里加入，利用计算出的节点电压及结构体中的信息 noiteration
Nodea=Nodea' %输出母线电压kV
NodeV=NodeV' %输出母线电压相角/度 End
%计算偏节点PQ差量
18
华南理工大学广州学院电气工程学院电力系统分析课程设计报告
function
[deltaP,deltaQ,maxdP,maxdQ]=FormdPQvector(Bus,NodeV,Nodea,dQGQL,dPGPL,B,G,busnumber)
deltaQ=dQGQL; %置初值，为无功发电减去无功负荷 deltaP=dPGPL; %置初值，为有功发电减去有功负荷 maxdP=0.; %置初值
maxdQ=0.; %置初值 for i=1:busnumber
if Bus(i).type==1 %平衡节点 deltaQ(i)=0.; deltaP(i)=0.; end
if Bus(i).type==3 deltaQ(i)=0.; %y1=0; %y2=0; y3=0;
for k=1:busnumber
if (B(i,k)~=0|G(i,k)~=0)
%y1=y1+(G(i,k)*VX(k)-B(i,k)*VY(k)); %y2=y2+(G(i,k)*VY(k)+B(i,k)*VX(k));
y3=y3+NodeV(k)*(G(i,k)*cos(Nodea(i)-Nodea(k))+B(i,k)*sin(Nodea(i)-Nodea(k)));
end end
deltaP(i)=deltaP(i)-y3*NodeV(i);
TltaP2(i)=(deltaP2(i)-(y1*VX(i)+y2*VY(i)))/Bus(i).V; end
if Bus(i).type==2 %PQ节点计算有功偏差和无功偏差 %y1=0; %y2=0; y3=0; y4=0;
for k=1:busnumber
if (B(i,k)~=0|G(i,k)~=0)
%y1=y1+(G(i,k)*VX(k)-B(i,k)*VY(k)); %y2=y2+(G(i,k)*VY(k)+B(i,k)*VX(k));
y3=y3+NodeV(k)*(G(i,k)*cos(Nodea(i)-Nodea(k))+B(i,k)*sin(Nodea(i)-Nodea(k)));
y4=y4+NodeV(k)*(G(i,k)*sin(Nodea(i)-Nodea(k))-B(i,k)*cos(Nodea(i)-Nodea(k)));
19
华南理工大学广州学院电气工程学院电力系统分析课程设计报告
end end
deltaP(i)=deltaP(i)-y3*NodeV(i);
TltaP2(i)=(deltaP2(i)-(y1*VX(i)+y2*VY(i)))/Bus(i).V; deltaQ(i)=deltaQ(i)-y4*NodeV(i);
TltaQ2(i)=(deltaQ2(i)-(y1*VY(i)-y2*VX(i)))/Bus(i).V; end
%以下程序断找出 deltaP和 deltaQ向量中的绝对值最大元素 if maxdP
if maxdQ
deltaP(i)=deltaP(i)/NodeV(i); deltaQ(i)=deltaQ(i)/NodeV(i); end end
%计算潮流分布函数 function
[Sf_T,St_T,Sf_L,St_L]=comput_flow_distribution(Bus,Acline,Trans,aclinenumber,transnumber,Sbase_MVA) if transnumber>0
for l=1:transnumber
Z=Trans(l).R+j*Trans(l).X; %变压器阻抗
f=Trans(l).fbno; % f 始节点编号
Vf=Bus(f).V*cos(Bus(f).angle)+j*Bus(f).V*sin(Bus(f).angle); %始节点电压相量
t=Trans(l).tbno; % t 终节点编号 Vt=(Bus(t).V*cos(Bus(t).angle)+j*Bus(t).V*sin(Bus(t).angle))/Trans(l).k; %终节点理想变压器内侧电压相量
I=(Vf-Vt)/Z; %变压器电流相量
Sf_T(l)=Vf*conj(I)*Sbase_MVA; %变压器始节点复功率/MVar St_T(l)=Vt*conj(I)*Sbase_MVA; %变压器终节点复功率/MVar Trans(l).fS=abs(Sf_T(l)); Trans(l).tS=abs(St_T(l)); end else
Sf_T(1)=0; St_T(1)=0; end pause
20
华南理工大学广州学院电气工程学院电力系统分析课程设计报告
for l=1:aclinenumber
Z=Acline(l).R+j*Acline(l).X %线路阻抗 f=Acline(l).fbno;
Vf=Bus(f).V*cos(Bus(f).angle)+j*Bus(f).V*sin(Bus(f).angle); %始节点电压相量
t=Acline(l).tbno;
Vt=Bus(t).V*cos(Bus(t).angle)+j*Bus(t).V*sin(Bus(t).angle); %终节点电压相量
I=(Vf-Vt)/Z; %线路电流相量
Sf_L(l)=(Vf*conj(I)-j*(Bus(f).V)^2*Acline(l).hB)*Sbase_MVA; %线路始节点复功率/MVar
St_L(l)=(Vt*conj(I)+j*(Bus(t).V)^2*Acline(l).hB)*Sbase_MVA; %线路终节点复功率/MVar
Acline(l).fS=abs(Sf_L(l)); Acline(l).tS=abs(St_L(l)); end end
六、 输入与输出结果
6.1 输入数据
节点数据（Nodedata.txt）
bus1 1 18. 0. 0. 0. 0. 0. bus2 3 18. 1.63 0. 0. 0. 0. bus3 3 18. 0.85 0. 0. 0. 0. bus4 2 230. 0. 0. 0. 0. 0. bus5 2 230. 0. 0. 1.25 0.5 0. bus6 2 230. 0. 0. 0.9 0.3 0. bus7 2 230. 0. 0. 0. 0. 0. bus8 2 230. 0. 0. 1.0 0.35 0. bus9 2 230. 0. 0. 0. 0. 0. 支路数据（Aclinedata.txt）
bus4 bus5 230. 0.01 0.085 0.088 bus4 bus6 230. 0.017 0.092 0.079 bus5 bus7 230. 0.032 0.161 0.153 bus6 bus9 230. 0.039 0.17 0.179 bus7 bus8 230. 0.0085 0.072 0.0745 bus8 bus9 230. 0.0119 0.1008 0.1045 变压器数据（Transdata.txt）
bus1 18.0 18.0 bus4 230. 230. 0.0 0.0576
21