电力系统分析潮流计算的计算机算法.doc

潮流计算的计算机算法实验报告姓名： 学号： 班级： 一、 实验目的掌握潮流计算的计算机算法。熟悉MATLAB，并掌握MATLAB程序的基本调试方法.二、 实验准备根据课程内容，熟悉MATLAB软件的使用方法,自行学习MATLAB程序的基础语法，并根据所学知识编写潮流计算牛顿拉夫逊法（或PQ分解法） 的计算程序，用相应的算例在MATLAB上进行计算、调试和验证。三、 实验要求每人一组，在实验课时内，用MATLAB调试和修改运行程序，用算例计算输出潮流结果。四、 实验程序clear; 清空内存n=input（请输入节点数：n=)；n1=input(请输入支路数:n1=)；isb=input（请输入平衡节点号：isb=）；pr=input(请输入误差精度：pr=);B1=input(请输入支路参数:B1=)；B2=input（请输入节点参数：B2=）;X=input（节点号和对地参数：X=)；Y=zeros（n）; Times=1； 一：创建节点导纳矩阵for i=1：n1 if 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。5B1（i，4）； Y(q，q)=Y（q,q）+1/B1（i，3)+0.5B1（i,4）； else 含有变压器的支路 p=B1（i,1)； q=B1(i，2); Y(p,q）=Y（p，q）1/(B1（i，3)B1(i,5)）; Y(q，p）=Y(p,q); Y(p,p)=Y（p，p）+1/B1（i，3）； Y(q,q）=Y（q,q）+1/(B1(i，5）2B1(i,3）); endendY； %将OrgS、DetaS初始化 OrgS=zeros（2*n-2，1）; DetaS=zeros（2*n2,1); %二：创建OrgS，用于存储初始功率参数h=0;j=0;for i=1：n 对PQ节点的处理 if i=isb&B2（i，6)=2 h=h+1; for j=1:n OrgS(2h1,1)=OrgS(2*h-1，1)+real（B2（i,3)）（real(Y（i,j)*real(B2（j,3）)-imag（Y(i,j)*I mag(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)）)； end endend三：对PV节点的处理，注意这时不可再将h初始化为0for i=1：n if i=isbB2（i，6)=3 h=h+1; for j=1：n 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（2h，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)）; end endendOrgS;%四：创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1，1）;t=0;for i=1:n if B2（i，6)=3 t=t+1; PVU（t，1)=B2（i,3); endendPVU；%五:创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n 对PQ节点的处理 if i=isbB2(i,6)=2 h=h+1; DetaS（2h1，1）=real（B2（i,2)OrgS（2h-1,1）; DetaS（2h，1)=imag(B2（i,2）)-OrgS（2h,1）； endendt=0；for i=1：n %六:对PV节点的处理,注意这时不可再将h初始化为0 if i=isb&B2(i，6)=3 h=h+1； t=t+1; DetaS（2*h-1,1）=real（B2（i,2）-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； endendDetaS；七：创建I，用于存储节点电流参数i=zeros(n1,1);h=0；for i=1：n if i=isb h=h+1； I(h，1）=(OrgS(2*h1，1）-OrgS（2*h,1）*sqrt（1）/conj(B2（i,3）)； endendI;八：创建Jacbi(雅可比矩阵）Jacbi=zeros（2*n2）;h=0;k=0;for i=1：n %对PQ节点的处理 if B2(i,6）=2 h=h+1; for j=1:n if j=isb k=k+1； if i=j %对角元素的处理 Jacbi（2h1,2k-1）=imag（Y（i，j）)*real（B2(i,3）)+real（Y（i，j）)imag(B2（i，3)+imag(I(h,1)）； Jacbi（2h1,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,2k-1)=-Jacbi（2*h1,2k）+2*real(I(h，1）; Jacbi(2*h,2k）=Jacbi（2h-1，2*k-1)-2*imag（I(h，1)）; else %非对角元素的处理 Jacbi（2*h-1,2*k1）=-imag(Y(i，j）real(B2（i,3)）+real（Y（i，j)*imag（B2(i,3)）; Jacbi（2*h-1，2k）=real(Y(i，j）real(B2（i，3）)+imag（Y(i,j）*imag（B2(i，3)； Jacbi（2h，2*k1)=Jacbi(2*h-1，2k）； Jacbi（2h，2k)=Jacbi（2*h-1,2*k1); end if k=(n-1） %将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0; end end end endendk=0；for i=1:n %对PV节点的处理 if B2(i,6)=3 h=h+1; for j=1：n if j=isb k=k+1； if i=j %对角元素的处理 Jacbi(2h1，2*k-1)=-imag(Y（i，j）)real（B2(i，3)）+real(Y（i，j）)*imag（B2(i，3）)+imag（I(h，1)）； Jacbi(2h1，2*k）=real(Y(i,j)*real(B2(i,3)）+imag（Y(i，j)imag（B2（i，3）+real（I(h,1）； Jacbi（2h，2*k1)=2*imag(B2（i，3）; Jacbi(2*h,2*k）=2real(B2（i，3)）; else 非对角元素的处理 Jacbi（2h1,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，2k1）=0； Jacbi（2*h,2*k)=0； end if k=(n-1） %将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0； end end end endendJacbi;九：求解修正方程，获取节点电压的不平衡量DetaU=zeros（2n-2，1）;DetaU=inv（Jacbi)DetaS;DetaU;修正节点电压j=0；for i=1：n %对PQ节点处理 if B2（i，6）=2 j=j+1; B2(i,3）=B2（i,3)+DetaU（2*j，1）+DetaU（2j-1,1)*sqrt（1）； endendfor i=1:n %对PV节点的处理 if B2(i，6)=3 j=j+1; B2（i,3)=B2(i，3)+DetaU（2*j,1)+DetaU（2j1，1)*sqrt(1); endendB2;十：开始循环*while abs(max（DetaU）prOrgS=zeros(2*n-2,1); %!！初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0；j=0；for i=1:n if i=isb&B2（i,6)=2 h=h+1； for j=1:n OrgS（2*h1，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（2h,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)）)； end endendfor i=1：n if i=isbB2（i，6）=3 h=h+1; for j=1:n OrgS(2h1，1）=OrgS（2*h1，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(2h，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)）; end endendOrgS；创建DetaSh=0；for i=1:n if i=isb&B2（i,6）=2 h=h+1； DetaS(2h1，1）=real（B2（i,2）OrgS(2h1,1）； DetaS（2*h，1）=imag（B2(i,2）OrgS（2*h，1)； endendt=0;for i=1：n if i=isbB2(i,6)=3 h=h+1； t=t+1； DetaS(2h-1，1)=real(B2(i，2）)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; endendDetaS；%创建Ii=zeros(n-1,1）；h=0;for i=1:n if i=isb h=h+1; I（h,1）=(OrgS（2*h-1,1)-OrgS(2h，1）*sqrt(1）)/conj(B2(i，3)）； endendI；创建JacbiJacbi=zeros(2*n-2)；h=0；k=0;for i=1：n if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1； if i=j Jacbi（2*h1,2*k1)=imag(Y（i,j）)*real（B2(i，3)）+real(Y（i,j）)*imag(B2(i,3)+imag(I（h，1）; Jacbi(2*h1,2*k)=real（Y（i，j）*real（B2(i，3）+imag(Y(i,j)）imag（B2（i,3）+real（I(h,1)； Jacbi（2h,2k-1）=Jacbi（2*h1,2*k)+2*real（I（h,1)）； Jacbi（2h,2*k)=Jacbi(2*h-1,2k1)2imag(I(h,1）; else Jacbi（2h-1，2*k1）=-imag（Y(i,j）real（B2(i,3）)+real(Y（i,j)）*imag(B2（i，3)）; Jacbi（2h1，2*k）=real(Y（i,j）)*real（B2（i,3)）+imag(Y（i，j)）imag（B2(i,3)； Jacbi（2h,2k1）=-Jacbi（2h1,2k）； Jacbi（2*h，2*k)=Jacbi(2h1,2k1); end if k=（n1） k=0; end end end endendk=0;for i=1:n if B2(i,6）=3 h=h+1; for j=1：n if j=isb k=k+1; if i=j Jacbi（2h1，2*k1)=-imag(Y（i,j）)real（B2（i，3)）+real（Y（i，j）)*imag（B2（i，3）)+imag(I（h，1）)； Jacbi（2h1，2k)=real（Y(i，j)*real（B2（i,3)）+imag(Y(i，j）)imag（B2(i,3)）+real（I(h，1)）； Jacbi(2h，2k1)=2*imag（B2（i，3）; Jacbi(2*h，2*k）=2real(B2（i，3）; else Jacbi（2*h-1，2k1)=imag（Y(i,j)*real(B2（i,3）)+real（Y(i，j）)*imag(B2（i，3）)； Jacbi（2*h-1,2k）=real（Y（i，j）real（B2（i,3)+imag(Y（i，j)*imag（B2(i,3）； Jacbi（2*h，2*k1）=0； Jacbi(2h，2*k)=0; end if k=(n1） k=0； end end end endendJacbi;DetaU=zeros（2n2，1）;DetaU=inv（Jacbi)*DetaS；DetaU；修正节点电压j=0；for i=1：n if B2(i,6)=2 j=j+1； B2(i,3)=B2（i,3）+DetaU（2*j,1）+DetaU（2*j1，1）*sqrt(1); endendfor i=1：n if B2（i,6)=3 j=j+1； B2（i,3）=B2（i,3）+DetaU（2*j,1)+DetaU（2*j-1,1）sqrt（1)； endendB2;Times=Times+1； %迭代次数加1endTimes；五、 实验流程六、实验结果参数输入：运行结果：七、实验体会通过这次实验，让我第一次接触到了MATLAB,并深切体会到了它的强大之处；潮流计算的计算机算法的实现不仅巩固了我的学过的知识，还让我学到一些MATLAB的编程,虽然在实验的过程中出现了很多的错误，但在老师的细心指导下，问题都解决啦；计算机为我们省去了大量的人工计算，希望在以后的学习中能接触到更多的软件，学习到更多的知识.