潮流上机课程设计报告及9节点、短路算法visual studio源程序.doc

课程设计报告( 2011- 2012年度第1学期)名 称：电力系统潮流上机 院 系：电气与电子工程学院 班 级： 学 号： 学生姓名： 指导教师： 设计周数： 两周 成 绩： 日期： 2012年1月5日 课程设计报告 电力系统综合仿真课程设计一、 目的与要求培养学生编程实现电力系统仿真计算的能力，掌握潮流计算、短路计算、暂态稳定计算的相关知识二、 主要内容1 编写潮流计算程序，要求如下：2.1据给定的潮流计算任务书整理潮流计算的基础数据：节点的分类，线路模型，等值变压器模型，电压等级的归算，标幺值的计算；2.2基础数据的计算机存储：节点数据，支路数据（包括变压器）；2.3用牛顿-拉夫逊法计算；2.4根据所选潮流计算方法画流程图，划分出功能模块，有数据输入模块，导纳阵形成模块，解线性方程组模块，计算不平衡功率模块，形成雅可比矩阵模块，解修正方程模块，计算线路潮流，网损，PV节点无功功率和平衡节点功率，数据输出模块；2.5据上述模块编制程序并上机调试程序，得出潮流计算结果；2.6源程序及其程序中的符号说明集、程序流图2 利用潮流计算结果编写故障计算程序，要求如下：2.1 发电机参数一律为PN=200MW，cos=0.85, ,2.2 变压器一律为Yd11接线，其中三角形侧接发电机，星形侧都直接接地2.3 线路零序参数一律为2.4 负荷都采用Y/Y型变压器接入，中性点与系统相连侧部接地。负荷按照对称负荷考虑，纯阻抗形式。2.5 计算5节点发生三相对称短路、单相接地短路、两相短路和两相接地短路时的故障电流，以及各支路电流和节点电压。3思考题3.1潮流计算和短路电流计算的元件模型有哪些不同点？节点导纳矩阵有哪些不同？计算过程有哪些不同点？3.2 短路电流计算中进行了哪些假设？3.3如果交给你一个任务，请你用已有的潮流计算软件计算北京城市电网的潮流，你应该做哪些工作？（收集哪些数据，如何整理，计算结果如何分析）3.4比较过去的计算机语言上机学习和本课程上机有哪些不同？设计中遇到的问题和解决的办法。三、 进度计划序号设计内容完成时间备注1教师讲解本次设计的任务及相关知识课程设计开始前2整理潮流计算的原始数据，完成数据的读入读出模块一天3编写及调试形成导纳矩阵模块一天4计算不平衡功率，形成雅可比矩阵一天5解修正方程一天6计算线路潮流，网损，PV节点无功功率和平衡节点功率一天7读入故障计算需要的数据，计算三相短路一天8计算不对称短路一天9写设计报告一天10验收一天四、 设计成果要求1 提交任务书2 计算机计算潮流程序流程图；3 计算机计算短路电流程序流程图；4 完整的仿真程序(包括变量的定义，子程序的调用等，计算结果)，并有程序注释。5 回答思考题五、 考核方式 总成绩=考勤+设计报告+面试成绩 学生姓名： 指导教师： 2012年1月5日附录： 1、计算机计算潮流程序流程图2、计算机计算三相短路电流程序流程图2、计算机计算三相短路电流程序流程图3、计算程序见附件程序包4、课程学习总结两周的时间，却恍若满满一学期。两周中学了太多东西，更确切地说，是学懂了太多东西，不论是课程知识还是学习方法。一、做题VS学懂我才终于明白，终于体会到，会做题和学会了是完全不同的概念。做题的时候只是一个两节点线路，而且省略各种支路，没有各种变压器、发电机那样复杂，感觉会做题了，貌似会了，但面对这次的9节点、变压器变比不为1的14节点和30节点潮流计算，我凌乱了认真听老师讲解却还是有地方跟不上，一边开始接触visual studio，一边重新学习理解课本上讲的东西。从数据的读入、输出到数组、负数到解方程，从节点导纳到雅可比各个元素再到正序增广网，上午在机房上完机去回寝室接着编，接着看书，一点点明白，一点点学会，程序一点点增多。每天都在体会着程序运行成功的喜悦，经历着运行失败又怎么都找不出错误的懊恼，我亦见同学看着别人都做出来了而自己还找不到错误，一边用力拍着头一边说着脏字，“怎么还做不出来！啊啊！”垂首顿足。高中以来养成的程式化做题模式和应试观念依然在影响着自己，虽然平时也有很认真的学习思考，但依然摆脱不了做题的阴影，总得做些题才放心。题是需要做的，做题也确实能找到不足，但却也容易陷入误区，因为毕竟题目是死的，题型也只有那么些，考试之前做几套往年的试卷，基本上考试就没问题了。而现在做潮流上机，一切都得从根上想明白，变压器变究竟影响哪些量，什么需要归算，究竟用有名值还是标幺值，节点导纳矩阵究竟怎样生成，哪个H、N、J、L对应Jacobi的哪个元素，错误一个一个出现，又一个一个被改正，一点点更深地明白牛顿拉夫逊法计算潮流的过程，一点点更深的理解复杂电力系统网络的模型。真的，会做题了，不等于学名白了。二、一丝不苟VS不拘小节两个意义完全相反的词，在这次潮流计算中却被我深刻的体会着。编写程序时必须一丝不苟，变量的定义和赋初值，该有和不该有的小括号，甚至一个很小的标点符号、很小的拼写错误都能让自己垂首顿足找不到错在哪。写程序真的很能锻炼人一丝不苟、认真细致的学习、做事品质。潮流计算的过程会忽略很多对结果影响不大而计算却很麻烦的数据，这就是所谓不拘小节。什么时候忽略对地支路，什么时候忽略支路电阻，认为正序阻抗与负序相等，变压器接法对零序阻抗求解的影响及其离此支路的忽略。做电力系统计算，脑子里必须清楚什么量什么时候可以忽略，什么时候不能忽略，而这前提是对概念非常清楚。三、备份的重要性编程期间受到了没及时备份的惩罚。我同时打开了两个VS程序，一个是已经做好的，又想换一种做法，结果新做法没弄好，关掉程序，第二天打开时全变成了最后的新程序，我凌乱啊然后重新写，写好了有问题怎么都找不出来，没办法又重写才好了，而至今也没有找出问题所在。如果备份了的话不仅是编程序，写文章、用CAD画工程图等这些费时费力一点点积累的东西，一定要做好备份。不仅是做这些事，其他的事也是一样，做之前要有准备，有多个准备。破釜沉舟固然使得秦关归楚，但不给自己留后路不是一个英明的决策者和实践家的作风。总之，两周的潮流上机胜似两月，明白了很多，总结了很多，获益，匪浅。谢谢这次课设，更谢谢老师！ 实验电09班 2012年1月5日9节点程序：/ 9LF2.cpp: 主项目文件。#include stdafx.h#include math.h#include #include #include using namespace System;using namespace std;#include nequation.hstruct BusDataint ID , Type;double U , A , PL , QL , PG , QG , U0 , G , B;*Bus;struct BranchDataint BusA , BusB , Type ;double R , X , B , K ;*Branch;int main(array args)/读入节点和支路数据ReadData()int BusNum,BranchNum;ifstream f1(D:1091181321Bus.txt);f1 BusNum;Bus = new BusDataBusNum;for(int i = 0;i Busi.ID Busi.Type Busi.U Busi.A Busi.PL Busi.QL Busi.PG Busi.QG Busi.U0 Busi.G Busi.B ;ifstream f2(D:1091181321Branch.txt);f2 BranchNum;Branch = new BranchDataBranchNum;for(int i = 0;i Branchi.BusA Branchi.BusB Branchi.Type Branchi.R Branchi.X Branchi.B Branchi.K ;/形成节点导纳矩阵FormY()complex *Y;Y = new complex BusNum*BusNum;for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)Yi*BusNum+j=0;for(int k=0;kBranchNum;k+)complex ZL(Branchk.R ,Branchk.X );complex YL=1.0/ZL;int i=Branchk.BusA-1;int j=Branchk.BusB-1;if(Branchk.K!=0)Yi*BusNum+i+=YL/complex(Branchk.K,0)/complex(Branchk.K,0) + complex(0,Branchk.B/2);Yi*BusNum+j-=YL/complex(Branchk.K,0);Yj*BusNum+i-=YL/complex(Branchk.K,0);Yj*BusNum+j+=YL + complex(0,Branchk.B/2);else Yi*BusNum+i+=YL + complex(0,Branchk.B/2);Yi*BusNum+j-=YL;Yj*BusNum+i-=YL;Yj*BusNum+j+=YL + complex(0,Branchk.B/2);for(int i=0;iBusNum;i+)Yi*BusNum+i+=complex( Busi.G ,Busi.B );ofstream f3(D:9LF2daona.txt);for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)f3 Yi*BusNum+j ;if(j=(BusNum-1)f3endl;f3.close();/定义第一次迭代时的初值并赋值double *e,*f,*U,*A,*P,*Q;e=new double BusNum;f=new double BusNum;U=new double BusNum;A=new double BusNum;P=new double BusNum;Q=new double BusNum;for(int i=0;iBusNum;i+)ei=Busi.U0;fi=0;Ui=sqrt(ei*ei+fi*fi); Ai=Math:Atan(fi/ei)*180/Math:PI;Pi=Busi.PG/100-Busi.PL/100;Qi=Busi.QG/100-Busi.QL/100;/计算不平衡功率for(int step=0;step10;step+)/定义并赋初值double *H,*N,*J,*L,*R,*S;double *aii,*bii;double *deltaP,*deltaQ;double *Jacobi;H=new double BusNum*BusNum;N=new double BusNum*BusNum;J=new double BusNum*BusNum;L=new double BusNum*BusNum;R=new double BusNum*BusNum;S=new double BusNum*BusNum;aii=new double BusNum;bii=new double BusNum;deltaP=new double BusNum;deltaQ=new double BusNum;int num=2*BusNum;Jacobi=new double num*num;for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)Hi*BusNum+j=0;Ni*BusNum+j=0;Ji*BusNum+j=0;Li*BusNum+j=0;Ri*BusNum+j=0;Si*BusNum+j=0;aiii=0;biii=0;deltaPi=0;deltaQi=0;Jacobii*BusNum+j=0;/计算不平衡功率for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)aiii+=Yi*BusNum+j.real()*ej-Yi*BusNum+j.imag()*fj;biii+=Yi*BusNum+j.real()*fj+Yi*BusNum+j.imag()*ej;for(int i=0;iBusNum;i+)deltaPi=Pi-(ei*aiii+fi*biii);deltaQi=Qi-(fi*aiii-ei*biii);ofstream f4(D:9LF2deltaPQ.txt);for(int i=1;iBusNum;i+)f4 deltaPindeltaQiendl;f4.close();/形成Jacobi矩阵for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)Hi*BusNum+j=-Yi*BusNum+j.imag()*ei+Yi*BusNum+j.real()*fi;Ni*BusNum+j= Yi*BusNum+j.real()*ei+Yi*BusNum+j.imag()*fi;Ji*BusNum+j=-Ni*BusNum+j;Li*BusNum+j= Hi*BusNum+j;for(int i=0;iBusNum;i+)Hi*BusNum+i+=biii;Ni*BusNum+i+=aiii;Ji*BusNum+i+=aiii;Li*BusNum+i-=biii;Ri*BusNum+i=2*fi;Si*BusNum+i=2*ei;for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)Jacobi(2*i)*num+(2*j)=Hi*BusNum+j;Jacobi(2*i)*num+(2*j+1)=Ni*BusNum+j;if(Busi.Type=2)Jacobi(2*i+1)*num+(2*j)=Ri*BusNum+j;Jacobi(2*i+1)*num+(2*j+1)=Si*BusNum+j;else Jacobi(2*i+1)*num+(2*j)=Ji*BusNum+j;Jacobi(2*i+1)*num+(2*j+1)=Li*BusNum+j;ofstream f5(D:9LF2Jacobi.txt);for(int i=2;inum;i+)for(int j=2;jnum;j+)f5Jacobii*num+j ;if(j=(num-1)f5endl;f5.close();/解方程求解新的fi,eiNEquation ob;int num2=num-2; /虽定义的Jacobi矩阵为num阶，但输出到文本的是去掉前两行的，这样解方程时就可直接读取ob.SetSize(num2);ifstream f6(D:9LF2Jacobi.txt);for(int i=0;inum2;i+)for(int j=0;job.Data(i,j);for(int i=0;iBusNum-1;i+)ob.Value(2*i) = deltaPi+1;if(Busi+1.Type = 2)ob.Value(2*i+1) = Busi+1.U * Busi+1.U - (ei+1 * ei+1 + fi+1 * fi+1);elseob.Value(2*i+1) = deltaQi+1;ob.Run();for(int i=0; iBusNum-1; i+)ei+1 += ob.Value(i*2+1);fi+1 += ob.Value(i*2);Ui+1 = sqrt(ei+1*ei+1+fi+1*fi+1); Ai+1=Math:Atan(fi+1/ei+1)*180/Math:PI;ofstream f7(D:9LF2UiAi.txt);for(int i=0;iBusNum;i+)f7Ui+1=Ui Ai+1=Aiendl;f7.close();ofstream f9(D:9LF2fiei.txt);for(int i=0;iBusNum;i+)f9fi eiendl;f9.close(); /计算线路潮流，网损，PV节点无功功率和平衡节点功率/平衡节点功率complex S1;for(int i=0;iBusNum;i+)S1+=complex(e0,f0)*complex(Y0*BusNum+i.real(),-Y0*BusNum+i.imag()*complex(ei,-fi);/线路功率complex *Sij,*Sji;Sij=new complexBranchNum;Sji=new complexBranchNum;for(int i=0;iBranchNum;i+)Siji=0;Sjii=0;for(int i=0;iBranchNum;i+)Siji=complex(eBranchi.BusA-1,fBranchi.BusA-1)*(complex(eBranchi.BusA-1,-fBranchi.BusA-1)*complex(0,-Branchi.B/2-Busi.B/2)+complex(eBranchi.BusA-1-eBranchi.BusB-1,fBranchi.BusB-1-fBranchi.BusA-1)*(1.0/complex(Branchi.R,-Branchi.X);Sjii=complex(eBranchi.BusB-1,fBranchi.BusB-1)*(complex(eBranchi.BusB-1,-fBranchi.BusB-1)*complex(0,-Branchi.B/2-Busi.B/2)+complex(eBranchi.BusB-1-eBranchi.BusA-1,fBranchi.BusA-1-fBranchi.BusB-1)*(1.0/complex(Branchi.R,-Branchi.X);/网损complex deltaS;for(int i=0;iBranchNum;i+)deltaS += Siji+Sjii;/PV节点的无功功率：前已求得，在输出时判断为PV节点并输出即可ofstream f8(D:9LF2Sij.txt);for(int i=0;iBranchNum;i+)f8 线路Branchi.BusABranchi.BusB的功率为：Siji.real()+jSiji.imag()nn 线路Branchi.BusBBranchi.BusA的功率为：Sjii.real()+jSjii.imag()nn;f8网损为：deltaSnn;for(int i=0;iBusNum;i+)if(Busi.Type=2)f8PV节点Busi.ID的无功功率为：Qinn;f8平衡节点功率为：S1endl;f8.close(); Console:WriteLine(LHello World); return 0;短路计算的程序：/ 3duanlu1.cpp: 主项目文件。#include stdafx.h#include math.h#include #include #include using namespace System;using namespace std;#include nequation.hstruct BusDataint ID , Type;double U , A , PL , QL , PG , QG , U0 , G , B;*Bus;struct BranchDataint BusA , BusB , Type ;double R , X , B , K ;*Branch;int main(array args) int BusNum=0,BranchNum=0;ifstream f1(D:3duanlu1Bus.txt);f1BusNum;Bus = new BusDataBusNum;for(int i = 0;i Busi.ID Busi.Type Busi.U Busi.A Busi.PL Busi.QL Busi.PG Busi.QG Busi.U0 Busi.G Busi.B ;ifstream f2(D:3duanlu1Branch.txt);f2 BranchNum;Branch = new BranchDataBranchNum;for(int i = 0;i Branchi.BusA Branchi.BusB Branchi.Type Branchi.R Branchi.X Branchi.B Branchi.K ;complex *Y,*U,*I,Id,*Z,*deltaU;double *f,*e;Y=new complex BusNum*BusNum;U=new complex BusNum;deltaU=new complex BusNum;f=new doubleBusNum;e=new doubleBusNum;Z=new complex BusNum;ifstream f3(D:3duanlu1fiei.txt);for(int i=0;ifiei;Ui=complex(ei,fi);ifstream f4(D:3duanlu1daona.txt);for(int i=0;iBusNum;i+)for(int j=0;j Yi*BusNum+j;/节点导纳矩阵的修改Y0*BusNum+0+=complex(0,1.0/(-0.19*100*0.85/200);Y1*BusNum+1+=complex(0,1.0/(-0.19*100*0.85/200);Y2*BusNum+2+=complex(0,1.0/(-0.19*100*0.85/200);Y4*BusNum+4+=complex(Bus4.PL/100-Bus4.PG/100,-Bus4.QL/100+Bus4.QG/100)/(U4*complex(U4.real(),-U4.imag();Y5*BusNum+5+=complex(Bus5.PL/100-Bus5.PG/100,-Bus5.QL/100+Bus5.QG/100)/(U5*complex(U5.real(),-U5.imag();Y7*BusNum+7+=complex(Bus7.PL/100-Bus7.PG/100,-Bus7.QL/100+Bus7.QG/100)/(U7*complex(U7.real(),-U7.imag();/三相对地短路故障/求解短路点的ZNEquation ob;ob.SetSize(BusNum);for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)ob.Data(i,j)=Yi*BusNum+j;for(int i=0;iBusNum;i+)ob.Value(i) = 0;ob.Value(4)=1;ob.Run();for(int i=0;iBusNum;i+)Zi=ob.Value(i);/计算故障电流IdId=-U4/Z4;/计算deltaUfor(int i=0;iBusNum;i+)deltaUi=Zi*Id;ofstream f5(D:3duanlu1deltaUi.txt);for(int i=0;iBusNum;i+)f5deltaUiendl;for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)f5 Yi*BusNum+j ;if(j=(BusNum-1)f5endl;f5.close();for(int i=0;iBusNum;i+)Ui+=deltaUi; /计算支路电流I=new complex BranchNum;for(int i=0;iBranchNum;i+)Ii=(UBranchi.BusA-1-UBranchi.BusB-1)/complex(Branchi.R,Branchi.X);ofstream f6(D:3duanlu13xiang.txt);f6短路电流幅值sqrt(Id.real()*Id.real()+Id.imag()*Id.imag() 相角Math:Atan(Id.imag()/Id.real()*180/Math:PInendl;for(int i=0;iBranchNum;i+)f6支路Branchi.BusA Branchi.BusB 电流幅值 sqrt(Ii.real()*Ii.real()+Ii.imag()*Ii.imag() 相角 Math:Atan(Ii.imag()/Ii.real()*180/Math:PIendl;f6n;for(int i=0;iBusNum;i+)f6节点i+1电压幅值sqrt(Ui.real()*Ui.real()+Ui.imag()*Ui.imag() 相角Math:Atan(Ui.imag()/Ui.real()*180/Math:PIendl;f6.close();/单相对地短路故障/正、负序阻抗与三相时相同，下面计算零序阻抗/重新形成零序导纳阵complex *Y0;Y0 = new complex BusNum*BusNum;for(int i=0;iBusNum;i+)for(int j=0;jBusNum;j+)Y0i*BusNum+j=0;for(int k=0;kBranchNum;k+)if(Branchk.Type=2)complex ZL(Branchk.R ,Branchk.X );complex YL=1.0/ZL;int i=Branchk.BusA-1;int j=Branchk.BusB-1;Y0j*BusNum+j+=complex(0.00001,0.00001) + complex(0,Branchk.B/2);Y0i*BusNum+j-=complex(0.00001,0.00001);Y0j*BusNum+i-=complex(0.00001,0.00001);Y0i*BusNum+i+=YL+complex(0,Branchk.B/2);elsecomplex ZL(Branchk.R ,3.0*Branchk.X );complex YL=1.0/ZL;int i=Branchk.BusA-1;int j=Branchk.BusB-1;Y0i*BusNum+i+=YL+complex(0,Branchk.B/2);Y0i*BusNum+j-=YL;Y0j*BusNum+i-=YL;Y0j*BusNum+j+=YL+complex(0,Branchk.B/2);for(int i=0;iBusNum;i+)Y0i*BusNum+i+=complex( Busi.