电力系统分析潮流计算实验报告姓名：XXXXXX 学号：XXXXXXXXXX 班级：XXXXXXXX一、实验目得掌握潮流计算计算机算法得方法，熟悉MATLAB得程序调试方法。

二、实验准备根据课程内容，熟悉MATLAB软件得使用方法，自行学习MATLAB程序得基础语法，并根据所学知识编写潮流计算牛顿拉夫逊法(或PQ分解法) 得计算程序，用相应得算例在MATLAB上进行计算、调试与验证。

三、实验要求每人一组，在实验课时内，调试与修改运行程序，用算例计算输出潮流结果。

四、程序流程五、实验程序%本程序得功能就是用牛拉法进行潮流计算%原理介绍详见鞠平著《电气工程》%默认数据为鞠平著《电气工程》例8、4所示数据%B1就是支路参数矩阵%第一列与第二列就是节点编号。

节点编号由小到大编写%对于含有变压器得支路，第一列为低压侧节点编号，第二列为高压侧节点编号%第三列为支路得串列阻抗参数，含变压器支路此值为变压器短路电抗%第四列为支路得对地导纳参数，含变压器支路此值不代入计算%第五烈为含变压器支路得变压器得变比，变压器非标准电压比%第六列为变压器就是否就是否含有变压器得参数，其中“1”为含有变压器，“0”为不含有变压器%B2为节点参数矩阵%第一列为节点注入发电功率参数%第二列为节点负荷功率参数%第三列为节点电压参数%第四列%第五列%第六列为节点类型参数，“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数%X为节点号与对地参数矩阵%第一列为节点编号%第二列为节点对地参数%默认算例% n=4;% n1=4;% isb=4;% pr=0、00001;% B1=[1 2 0、1667i 0 0、8864 1;1 3 0、1302+0、2479i 0、0258i 1 0;1 4 0、1736+0、3306i 0、0344i 1 0;3 4 0、2603+0、4959i 0、0518i 1 0];% B2=[0 0 1 0 0 2;0 -0、5-0、3i 1 0 0 2;0、2 0 1、05 0 0 3;0 -0、15-0、1i 1、05 0 0 1];% X=[1 0;2 0、05i;3 0;4 0];clear;clc;num=input('就是否采用默认数据？(1-默认数据;2-手动输入)');if num==1n=4;n1=4;isb=4;pr=0、00001;B1=[1 2 0、1667i 0 0、8864 1;1 3 0、1302+0、2479i 0、0258i 1 0;1 4 0、1736+0、3306i 0、0344i 1 0;3 4 0、2603+0、4959i 0、0518i 1 0];B2=[0 0 1 0 0 2;0 -0、5-0、3i 1 0 0 2;0、2 0 1、05 0 0 3;0 -0、15-0、1i 1、05 0 0 1];X=[1 0;2 0、05i;3 0;4 0];elsen=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号与对地参数:X=');endTimes=1; %迭代次数%创建节点导纳矩阵Y=zeros(n);for i=1:n1if B1(i,6)==0 %不含变压器得支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0、5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0、5*B1(i,4);else %含有变压器得支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-B1(i,5)/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+B1(i,5)/B1(i,3)+(1-B1(i,5))/B1(i,3);Y(q,q)=Y(q,q)+B1(i,5)/B1(i,3)+(B1(i,5)*(B1(i,5)-1))/B1(i,3);endendfor i=1:n1Y(i,i)=Y(i,i)+X(i,2); %计及补偿电容电纳enddisp('导纳矩阵为：');disp(Y); %显示导纳矩阵%初始化OrgS、DetaSOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%创建OrgS，用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点得处理if i~=isb&B2(i,6)==2 %不就是平衡点&就是PQ点h=h+1;for j=1:n%公式8-74%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j ,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:n %对PV节点得处理，注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3 %不就是平衡点&就是PV点h=h+1;for j=1:n%公式8-75-a%Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej)%Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej)OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real (Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j ))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendend%创建PVU 用于存储PV节点得初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endend%创建DetaS，用于存储有功功率、无功功率与电压幅值得不平衡量h=0;for i=1:n %对PQ节点得处理if i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); %delPiDetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1); %delQiendendt=0;for i=1:n %对PV节点得处理，注意这时不可再将h初始化为0if i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1); %delPiDetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2; %delUi endend% DetaS%创建I，用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));%conj求共轭endend%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点得处理if B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素得处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else %非对角元素得处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1) %将用于内循环得指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendfor i=1:n %对PV节点得处理if B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素得处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非对角元素得处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环得指针置于初始值，以确保雅可比矩阵换行k=0;endendendendenddisp('初始雅可比矩阵为：');disp(Jacbi);%求解修正方程，获取节点电压得不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %inv矩阵求逆% DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理if B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n %对PV节点得处理if B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end% B2%开始循环**********************************************************************while abs(max(DetaU))>prOrgS=zeros(2*n-2,1);h=0;j=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real (Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j ))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:nif i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real (Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j ))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendend% OrgS%创建DetaSh=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,6)==3h=h+1;t=t+1;% DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;endend% DetaS%创建Ii=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endend% I%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendend% JacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;% DetaU%修正节点电压j=0;for i=1:nif B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:nif B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endend% B2Times=Times+1; %迭代次数加1enddisp('迭代次数为：');disp(Times);disp('收敛时电压修正量为：：');disp(DetaU);for k=1:nE(k)=B2(k,3);e(k)=real(E(k));f(k)=imag(E(k));V(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)、/e(k))*180、/pi;end%=============== 计算各输出量===========================disp('各节点得实际电压标幺值E为(节点号从小到大排列)：');disp(E); %显示各节点得实际电压标幺值E用复数表示disp('-----------------------------------------------------')disp('各节点得电压大小V为(节点号从小到大排列)：');disp(V); %显示各节点得电压大小V得模值disp('-----------------------------------------------------');disp('各节点得电压相角sida为(节点号从小到大排列)：');disp(sida); %显示各节点得电压相for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q)); %计算各节点得注入电流得共轭值endS(p)=E(p)*C(p); %计算各节点得功率S = 电压X 注入电流得共轭值enddisp('各节点得功率S为(节点号从小到大排列)：');disp(S); %显示各节点得注入功率Sline=zeros(n1,5);disp('-----------------------------------------------------');disp('各条支路得首端功率Si为(顺序同您输入B1时一致)：');for i=1:n1p=B1(i,1);q=B1(i,2);Sline(i,1)=B1(i,1);Sline(i,2)=B1(i,2);if B1(i,6)==0Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)、/2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1、/(B1(i,3)*B1(i,5))));Siz(i)=Si(p,q);elseSi(p,q)=E(p)*(conj(E(p))*((1-B1(i,5))/B1(i,3))+(conj(E(p))-conj(E(q)))*(B1(i,5)/B1(i,3)));Siz(i)=Si(p,q);endSSi(p,q)=Si(p,q);Sline(i,3)=Siz(i);ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);enddisp('-----------------------------------------------------');disp('各条支路得末端功率Sj为(顺序同您输入B1时一致)：');for i=1:n1p=B1(i,1);q=B1(i,2);if B1(i,6)==0Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)、/2)+(conj(E(q)、/B1(i,5))-conj(E(p)))*conj(1、/(B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);elseSj(q,p)=E(q)*(conj(E(q))*((B1(i,5)*(B1(i,5)-1))/B1(i,3))+(conj(E(q))-conj(E(p)))*(B1(i,5)/B1(i,3)));Sjy(i)=Sj(q,p);endSSj(q,p)=Sj(q,p);Sline(i,4)=Sjy(i);ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);enddisp('-----------------------------------------------------');disp('各条支路得功率损耗DS为(顺序同您输入B1时一致)：');for i=1:n1p=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);DDS(i)=DS(i);Sline(i,5)=DS(i);ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);enddisp('-----------------------------------------------------');disp('各支路首端编号末端编号首端功率末端功率线路损耗');disp(Sline);六、运行结果及其分析就是否采用默认数据？(1-默认数据;2-手动输入)1导纳矩阵为：2、9056 -11、5015i 0、0000 + 5、3173i -1、6606 +3、1617i -1、2450 + 2、3710i 0、0000 + 5、3173i 0、0000 -4、6633i 0、0000 + 0、0000i 0、0000 + 0、0000i -1、6606 + 3、1617i 0、0000 + 0、0000i 2、4904 - 4、7039i -0、8298 + 1、5809i -1、2450 + 2、3710i 0、0000 + 0、0000i -0、8298 + 1、5809i 2、0749 - 3、9089i 初始雅可比矩阵为：11、1267 2、7603 -5、3173 0 -3、1617 -1、6606-3、0509 11、8762 0 -5、3173 1、6606 -3、1617-5、3173 0 5、3173 0 0 00 -5、3173 0 4、0092 0 0-3、3198 -1、7436 0 0 4、8217 2、69800 0 0 0 0 2、1000迭代次数为：4收敛时电压修正量为：：1、0e-05 *0、0349-0、2445-0、0101-0、5713-0、0931-0、0073各节点得实际电压标幺值E为(节点号从小到大排列)：0、9673 - 0、0655i 1、0252 - 0、1666i 1、0495 - 0、0337i 1、0500 + 0、0000i -----------------------------------------------------各节点得电压大小V为(节点号从小到大排列)：0、9695 1、0387 1、0500 1、0500-----------------------------------------------------各节点得电压相角sida为(节点号从小到大排列)：-3、8734 -9、2315 -1、8419 0各节点得功率S为(节点号从小到大排列)：-0、0000 + 0、0000i -0、5000 - 0、3000i 0、2000 + 0、1969i 0、3277 + 0、0443i-----------------------------------------------------各条支路得首端功率Si为(顺序同您输入B1时一致)：S(1,2)=-0、5-0、30713iS(1,3)=-0、24266-0、197iS(1,4)=-0、25734-0、11013iS(3,4)=-0、055551+0、0017528i-----------------------------------------------------各条支路得末端功率Sj为(顺序同您输入B1时一致)：S(2,1)=0、5+0、24606iS(3,1)=0、25555+0、1952iS(4,1)=0、2712+0、1014iS(4,3)=0、056496-0、057061i-----------------------------------------------------各条支路得功率损耗DS为(顺序同您输入B1时一致)：DS(1,2)=0-0、06107iDS(1,3)=0、012892-0、0018014iDS(1,4)=0、013863-0、0087295iDS(3,4)=0、00094545-0、055308i-----------------------------------------------------各支路首端编号末端编号首端功率末端功率线路损耗1、0000 + 0、0000i2、0000 + 0、0000i -0、5000 - 0、3071i 0、5000 + 0、2461i 0、0000 - 0、0611i1、0000 + 0、0000i 3、0000 + 0、0000i -0、2427 - 0、1970i 0、2556 + 0、1952i 0、0129 - 0、0018i1、0000 + 0、0000i 4、0000 + 0、0000i -0、2573 - 0、1101i 0、2712 + 0、1014i 0、0139 - 0、0087i3、0000 + 0、0000i4、0000 + 0、0000i -0、0556 + 0、0018i 0、0565 - 0、0571i 0、0009 - 0、0553i七、实验体会及感悟通过这次实验，首先让我对matlab软件有了初步得了解，对它强大得矩阵运算能力有了更深得体会，同时掌握了设置断点与断点调试得一般方法，结合课本上得程序流程图与参考资料上得例子单步跟踪调试，再一次得熟悉了牛顿拉夫逊法潮流计算得一般方法与步骤，对计算机计算潮流计算有了更进一步得认识，在学习潮流计算时，虽然依次学习了节点导纳矩阵，功率方程、雅可比矩阵，但不能将它们联系起来，更不知道其中得原委，通过程序得编写，知道了其中得联系，也知道了每个方程、矩阵在计算中得作用。