华北电力大学电力系统暂态上机计算课程设计报告2014
课程设计报告
(2013—2014年度第二学期)
名称：电力系统暂态上机计算
院系：电气与电子工程学院
班级：
学号：
学生姓名：
指导教师：刘宝柱
设计周数：两周
成绩：
日期：2014年7月11日
巩固电力系统暂态分析理论知识，使学生掌握采用计算机对电力系统电磁暂态过程和机电暂态过程进行计算的基本方法，并进一步巩固计算机编程能力，为将来从事相关的技术工作打下必要的基础。
2．设计内容
2.1计算和形成阻抗矩阵
2.2对称和不对称短路计算
2.3静态稳定性计算
2.4暂态稳定性计算
3．设计过程和成果
采用Matlab进行编程计算，按照各计算的详细流程图完成计算。
3.1阻抗矩阵计算和短路计算
3.1.1流程图
以下见阻抗矩阵计算和对称短路计算结果页：
14-7-5上午4:37MATLAB Command Window1/2；
3.1.3对称短路计算结果
以下见阻抗矩阵计算和对称短路计算结果页：
选取节点3的三相短路故障进行计算：
14-7-1下午3:27MATLAB Command Window1/2；
3.1.4不对称短路计算结果
以下见不对称短路计算结果页：
选取节点3的a相接地短路故障进行计算：
14-7-1下午3:28MATLAB Command Window1/2；
选取节点3的a相经10欧姆电阻接地短路故障进行计算：14-7-1下午3:28MATLAB Command Window1/2；
选取节点3的bc两相相间短路故障进行计算：
14-7-1下午3:28MATLAB Command Window1/2；
选取节点3的bc两相接地短路故障进行计算：
14-7-1下午3:28MATLAB Command Window1/2；
3.2.1流程图
3.2.2静态稳定性计算结果
以下见静态稳定性计算结果页：
按照《电力系统暂态分析（第三版）》第191页题目7-2进行计算：
14-7-1下午3:30MATLAB Command Window1/2；
3.3.1流程图
3.3.2暂态稳定性计算结果
以下见暂态稳定性计算结果页：
按照《电力系统暂态分析（第三版）》第164页题目6-1进行计算：
14-7-1下午3:30MATLAB Command Window1/2；
4.1计算短路电流，书中给出的手算方法与计算机编程方法有何区别？
答：书中的手算方法是将网络逐步进行化简到短路点并计算短路电流，进而再由短路点电流返回去重新推算各支路电流；而计算机编程则是通过节点导纳矩阵来计算阻抗从而求得电流的。
4.2用计算机方法，进行电力系统潮流计算和短路电流计算，有哪些区别？
答：已知条件不同，电力系统潮流计算需要根据已知条件将节点分为PV节点和PQ节点，通过不同节点的不同条件来计算各节点的潮流数据；而短路计算则不需要区分PV和PQ节点，而是根据短路节点位置指定该节点电压为0去计算。
4.3如果交给你一个任务，请你用已有的短路电流计算软件计算某地区电网某点发生短路后的短路电流，你应该做哪些工作？（收集哪些数据，如何整理，计算结果如何分析）
答：
收集数据：网路结构图（包含节点数，支路数，支路的首末节点等）、线路，变压器，发电机的正负序阻抗等等值参数、短路点位置、短路类型等数据。
整理方法：按照程序的格式要求整理数据，由于所编写Matlab程序采用了将数据定义在源码头部的方式组织数据，对于整理后的数据还要按格式写进源码文件中（而不是在程序运行时载入）才能计算。
结果分析：运行程序，观察得到的结果数据，数据符合数理常识且在大致误差范围内则认定分析正确。
根据所得的正确结果，结合正常运行状态下的潮流分析结果，分析差异特征，得出电网短路理论分析的文本结论。
4.4静态稳定的概念，暂态稳定的概念？
答：电力系统静态稳定是指电力系统受到小扰动后，不发生自发振荡或非周期性失步，自动恢复到初始运行状态的能力。
静态稳定是指系统受到小干扰影响后，不发生自发振荡或非周期性失步，自动回到初始运行状态的能力；暂态稳定是指系统受到打干扰后，能否经过暂态过程达到新的稳定运行的状态或回到原来的运行状态的能力。
4.5提高静态稳定、暂态稳定的措施？
提高静态稳定：
（1）采用自动调节励磁装置
（2）减小原件电抗
①采用分裂导线
②提高线路额定电压等级
③采用串联补偿装置
（3）改善系统结构和采用中间补偿设备
提高暂态稳定：
（1）故障的快速切除和自动重合闸的应用
（2）提高发电机输出的电磁功率
①对发电机施行强行励磁
②电气制动
③变压器中性点经小电阻接地
（3）减少原动机输出的机械功率
6.参考文献
1.《电力系统暂态过程》，常鲜戎、赵书强编，机械工业出版社，2010年1月，第1版；
2.《电力系统暂态分析》，李光琦，中国出版社，2007年，第三版。
3.《精通Matlab6.5》北京航空航天大学出版社，张志涌
附录：程序源码
见源码页：。
新版电力系统暂态上机计算课程设计报告完整无误版
课程设计报告( 2014—2015年度第二学期)名称：电力系统暂态上机计算院系：电气及电子工程学院班级：电气1211学号：1121210205学生姓名：郝阳指导教师：陈艳波设计周数：两周成绩：日期：2015年7月4日一、课程设计的目的及要求巩固电力系统暂态分析理论知识，使学生掌握采用计算机对电力系统电磁暂态过程和机电暂态过程进行计算的基本方法，并进一步巩固计算机编程能力，为将来从事相关的技术工作打下必要的基础。
二、设计正文（详细内容见附录,用A4纸，页数不限）1．对称短路计算过程流程图和计算结果2．不对称短路计算过程流程图和计算结果3．静态稳定性计算过程流程图和计算结果4．暂态稳定性计算过程流程图和计算结果5．思考题三、课程设计总结或结论本次电力系统暂态上机主要讨论了计算各种类型短路故障下系统网络中的电压电流分布以及电力系统的静态稳定、暂态稳定问题。
通过本次课程设计，本人对电力系统故障分析有了更深刻的理解，电力系统的故障时，大部分电磁量将随时间变化，描述其特性的是微分方程，这给分析计算带来一定困难。
在分析过程中通常尽量避免对微分方程直接求解，而是采用一定的工具和假设使问题得以简化，即把“微分方程代数化，暂态分析稳态化”。
在分析不对称故障时，各相之间电磁量的耦合使问题的分析更为复杂，此时常用的分析方法是采用对称分量法将不对称问题转化为对称问题来求解。
同时我对用来分析电力系统静态稳定的试探法，用来分析电力系统暂态稳定的改进欧拉法有了一些使用心得。
这及手算系统短路电流时使用的网络化简方法大大不同。
在学习中，参照潮流程序，我加深了对节点导纳矩阵建立方法的理解及学习，巩固了不同类型短路的短路电流计算方法，和序电压、序电流，相电压、相电流基于matlab软件的计算。
利用小干扰分析法判别静态稳定可以使用劳斯判据，也可以使用特征根判别的方法。
二者都能够判断，在书中的例题中使用的是劳斯判据，但是，特征根判别放法更适合在matlab软件的环境下使用，因为该软件提供了计算矩阵特征根的函数，使用起来简便易懂。
电力系统暂态的课程设计
电力系统暂态的课程设计一、课程目标知识目标：1. 学生能够理解电力系统暂态过程的基本概念，掌握暂态过程中的关键参数及其相互影响。
2. 学生能够描述常见电力系统暂态现象，如短路、故障、暂态稳定性问题等，并分析其原因及危害。
3. 学生能够掌握电力系统暂态稳定性分析的基本方法，并运用相关公式进行简单计算。
技能目标：1. 学生能够运用所学知识，分析实际电力系统暂态问题，提出解决策略。
2. 学生能够通过查阅资料，了解电力系统暂态过程的最新研究动态和发展趋势。
3. 学生能够利用仿真软件对电力系统暂态过程进行模拟，验证所学理论和方法。
情感态度价值观目标：1. 培养学生对电力系统暂态问题的关注和兴趣，激发其探索精神。
2. 增强学生的团队合作意识，使其在分析和解决电力系统暂态问题时，能够积极与他人交流与合作。
3. 培养学生严谨、务实的科学态度，使其在分析问题时，能够充分考虑到实际工况，注重实际应用。
课程性质：本课程为电力系统专业课程，旨在帮助学生深入了解电力系统暂态过程，掌握相关分析方法和解决策略。
学生特点：学生具备一定的电力系统基础知识，具有较强的学习能力和动手能力。
教学要求：结合学生特点，注重理论与实践相结合，提高学生的实际操作能力和解决问题的能力。
在教学过程中，注重启发式教学，引导学生主动思考，培养其创新意识。
二、教学内容1. 电力系统暂态过程基本概念：包括暂态过程的定义、分类及其在电力系统中的重要性。
- 教材章节：第二章 电力系统暂态过程概述- 内容列举：暂态过程的产生原因、暂态过程的基本类型及其特点。
2. 常见电力系统暂态现象及分析：如短路、故障、暂态稳定性问题等。
- 教材章节：第三章 常见电力系统暂态现象- 内容列举：短路故障的分类及影响、暂态稳定性问题的分析方法及其防治措施。
3. 电力系统暂态稳定性分析：介绍暂态稳定性分析的基本原理、方法和应用。
- 教材章节：第四章 电力系统暂态稳定性分析- 内容列举：暂态稳定性分析的基本理论、时域仿真法、直接法等分析方法。
电力系统暂态课程设计
电力系统暂态课程设计一、课程目标知识目标：1. 掌握电力系统暂态过程的基本概念、分类及特点；2. 理解电力系统暂态稳定性分析的基本原理和方法；3. 了解电力系统暂态过程中的主要影响因素及其作用机理。
技能目标：1. 能够运用所学知识对电力系统暂态过程进行正确分析；2. 能够运用暂态稳定性分析的方法，评估电力系统的稳定性；3. 能够根据实际案例，提出改善电力系统暂态稳定性的措施。
情感态度价值观目标：1. 培养学生对电力系统暂态过程研究的兴趣，激发学生探索精神；2. 增强学生的团队协作意识，培养合作解决问题的能力；3. 提高学生关注电力系统安全、稳定的意识，树立社会责任感。
本课程针对高年级学生，结合电力系统暂态过程的学科特点，注重理论与实践相结合。
通过本课程的学习，使学生能够深入理解电力系统暂态过程的本质，掌握相关分析方法和技能，为今后从事电力系统运行、维护和管理等工作奠定基础。
同时，培养学生关注社会热点问题，提高解决实际问题的能力。
课程目标的设定旨在使学生在知识、技能和情感态度价值观方面取得全面发展。
二、教学内容1. 电力系统暂态过程基本概念：暂态过程的定义、分类及其在电力系统中的重要性；教材章节：第一章第一节。
2. 电力系统暂态稳定性分析原理：暂态稳定性分析的基本理论、方法及其应用；教材章节：第二章。
3. 电力系统暂态过程中的影响因素：分析电力系统暂态过程中的主要因素，如发电机、变压器、线路等；教材章节：第三章。
4. 暂态稳定性分析的方法与案例：介绍暂态稳定性分析的方法，结合实际案例进行分析；教材章节：第四章。
5. 改善电力系统暂态稳定性的措施：探讨提高暂态稳定性的措施，如励磁系统、PSS、直流输电等；教材章节：第五章。
6. 电力系统暂态过程实验与仿真：组织学生进行暂态过程实验，利用仿真软件进行模拟分析；教材章节：第六章。
教学内容按照课程目标进行科学、系统地组织和安排，注重理论与实践相结合。
教学进度根据教材章节进行制定，确保学生能够在掌握基本理论和方法的基础上，通过实际案例和实验，深入理解电力系统暂态过程及其稳定性问题。
电力系统仿真实验指导书(华北电力大学2014)
华电潮流上机课程设计报告+程序
院系：电气与电子工程学院班级：学号：学生姓名：指导教师：刘宝柱设计周数：成绩：日期：2012年1月5日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细内容见附录）1．手算2．计算机计算3．思考题3.1潮流计算的方法有哪些？各有何特点？答：潮流计算的方法主要有高斯-赛德尔迭代法、牛顿-拉夫逊迭代法和P-Q分解法。
它们各自的特点如下：（1）高斯-赛德尔迭代法分为以节点导纳矩阵为基础的高斯-赛德尔迭代法和以以节点阻抗矩阵为基础的高斯-赛德尔迭代法的原理比较简单，要求的数字计算机的内存量也比较小，但收敛性差，当系统规模变大时，迭代次数急剧上升，往往出现迭代不收敛的情况；而阻抗法改善了电力系统潮流计算导纳法德收敛性问题，在当时获得了广泛的应用，但是，阻抗法的主要缺点是占用计算机的内存很大，每次迭代的计算量很大。
当系统不断扩大时，这些缺点就更加突出。
（2）牛顿-拉夫逊法是数学中求解非线性方程式的典型方法，有较好的收敛性。
只要在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿潮流计算程序的计算效率，牛顿法在收敛性、内存要求、计算速度方面都超过了阻抗法，成为知道目前仍被广泛采用的方法。
（3）P-Q分解法潮流计算派生于以极坐标表示时的牛顿-拉夫逊法，它根据电力系统的特点，抓住主要矛盾，对纯数学的牛顿法进行了改造。
与牛顿法相比，P-Q分解法的修正方程的系数矩阵B’和B”分别是（n-1）和（m-1）的方阵，替代了原有的（n+m-2）阶系数矩阵J；B’、B”在迭代过程中可以保持不变且为对称的系数矩阵，提高了计算速度，降低了对存储容量的要求。
P-Q分解法在计算速度方面有显著地提高，迅速得到了推广。
3.2如果交给你一个任务，请你用已有的潮流计算软件计算北京城市电网的潮流，你应该做哪些工作？（收集哪些数据，如何整理，计算结果如何分析）答：(1)在进行北京城市电网的潮流计算之前需要了解北京城市电网中所有的节点支路的相关数据，并对节点和支路分类。
新版电力系统暂态上机计算课程设计报告完整无误版附程序
课程设计报告( 2014—2015年度第二学期)名称：电力系统暂态上机计算院系：电气与电子工程学院班级：电气1211学号：1121210205学生姓名：郝阳指导教师：陈艳波设计周数：两周成绩：日期：2015年7月4日一、课程设计的目的与要求巩固电力系统暂态分析理论知识，使学生掌握采用计算机对电力系统电磁暂态过程和机电暂态过程进行计算的基本方法，并进一步巩固计算机编程能力，为将来从事相关的技术工作打下必要的基础。
二、设计正文（详细内容见附录,用A4纸，页数不限）1．对称短路计算过程流程图和计算结果2．不对称短路计算过程流程图和计算结果3．静态稳定性计算过程流程图和计算结果4．暂态稳定性计算过程流程图和计算结果5．思考题三、课程设计总结或结论本次电力系统暂态上机主要讨论了计算各种类型短路故障下系统网络中的电压电流分布以及电力系统的静态稳定、暂态稳定问题。
通过本次课程设计，本人对电力系统故障分析有了更深刻的理解，电力系统的故障时，大部分电磁量将随时间变化，描述其特性的是微分方程，这给分析计算带来一定困难。
在分析过程中通常尽量避免对微分方程直接求解，而是采用一定的工具和假设使问题得以简化，即把“微分方程代数化，暂态分析稳态化”。
在分析不对称故障时，各相之间电磁量的耦合使问题的分析更为复杂，此时常用的分析方法是采用对称分量法将不对称问题转化为对称问题来求解。
同时我对用来分析电力系统静态稳定的试探法，用来分析电力系统暂态稳定的改进欧拉法有了一些使用心得。
这与手算系统短路电流时使用的网络化简方法大大不同。
在学习中，参照潮流程序，我加深了对节点导纳矩阵建立方法的理解与学习，巩固了不同类型短路的短路电流计算方法，和序电压、序电流，相电压、相电流基于matlab软件的计算。
利用小干扰分析法判别静态稳定可以使用劳斯判据，也可以使用特征根判别的方法。
二者都能够判断，在书中的例题中使用的是劳斯判据，但是，特征根判别放法更适合在matlab软件的环境下使用，因为该软件提供了计算矩阵特征根的函数，使用起来简便易懂。
电力系统暂态上机计算课程设计报告
-0.5489+ 0.3169i
5
0.6091
0
0
0.5977+ 0.3451i
0.0000 - 0.6901i
-0.5977+ 0.3451i
c.各支路三序电流、三相电流：
支路
I1
I2
I0
Ia
Ib
Ic
1-2
0 + 0.8451i
0
0
0 + 0.8451i
0.7319- 0.4225i
3.《电力系统故障的计算机辅助分析》 大学 米麟书等
4.《电力系统潮流计算》 大学 宋文南等
5.《电力系统故障分析》 清华大学 周荣光
6.《短路电流实用计算方法》 电力工业 交通大学等
7.《精通Matlab6.5》航空航天大学，志涌
附录（设计流程图、计算结果、思考题答案）
1、对称短路计算流程图及计算结果
-0.7319- 0.4225i
1-3
0 - 4.5070i
0
0
0 - 4.5070i
-3.9032+ 2.2535i
3.9032+ 2.2535i
2-3
0 - 5.3521i
0
0
0 - 5.3521i
-4.6351+ 2.6761i
4.6351+ 2.6761i
1-4
0 + 3.6620i
0
0
0+3.6620i
另外，通过此次课程设计，我学习并掌握了一定的Matlab知识，对简单的选择语句、循环语句、数组的使用等有了一些了解，为今后学习、工作中再次学习使用此软件打下了基础，获益匪浅。
华电潮流上机课程设计报告程序..doc
院系：电气与电子工程学院班级：学号：学生姓名：指导教师：刘宝柱设计周数：成绩：日期：2012年1月5日课程课程设计报告一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细内容见附录）1．手算2．计算机计算3．思考题3.1潮流计算的方法有哪些？各有何特点？答：潮流计算的方法主要有高斯-赛德尔迭代法、牛顿-拉夫逊迭代法和P-Q分解法。
它们各自的特点如下：（1）高斯-赛德尔迭代法分为以节点导纳矩阵为基础的高斯-赛德尔迭代法和以以节点阻抗矩阵为基础的高斯-赛德尔迭代法的原理比较简单，要求的数字计算机的内存量也比较小，但收敛性差，当系统规模变大时，迭代次数急剧上升，往往出现迭代不收敛的情况；而阻抗法改善了电力系统潮流计算导纳法德收敛性问题，在当时获得了广泛的应用，但是，阻抗法的主要缺点是占用计算机的内存很大，每次迭代的计算量很大。
当系统不断扩大时，这些缺点就更加突出。
（2）牛顿-拉夫逊法是数学中求解非线性方程式的典型方法，有较好的收敛性。
只要在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿潮流计算程序的计算效率，牛顿法在收敛性、内存要求、计算速度方面都超过了阻抗法，成为知道目前仍被广泛采用的方法。
（3）P-Q分解法潮流计算派生于以极坐标表示时的牛顿-拉夫逊法，它根据电力系统的特点，抓住主要矛盾，对纯数学的牛顿法进行了改造。
与牛顿法相比，P-Q分解法的修正方程的系数矩阵B’和B”分别是（n-1）和（m-1）的方阵，替代了原有的（n+m-2）阶系数矩阵J；B’、B”在迭代过程中可以保持不变且为对称的系数矩阵，提高了计算速度，降低了对存储容量的要求。
P-Q分解法在计算速度方面有显著地提高，迅速得到了推广。
3.2如果交给你一个任务，请你用已有的潮流计算软件计算北京城市电网的潮流，你应该做哪些工作？（收集哪些数据，如何整理，计算结果如何分析）答：(1)在进行北京城市电网的潮流计算之前需要了解北京城市电网中所有的节点支路的相关数据，并对节点和支路分类。
潮流上机课程设计报告华电
《电力系统潮流上机》课程设计报告院系班级：学号：学生姓名：指导教师：设计周数成绩：日期：年月日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细内容见附录）1.手算： 要求应用牛顿-拉夫逊法或P-Q 分解法手算求解，要求精度为0.001MW 。
节点1为平衡节点，电压︒∠=00.11U ，节点2为PQ 节点，负荷功率6.08.0~2j S +=，节点3是PV 节点，1.1,4.033==U P ，两条支路分别为04.001.013j Z +=，2.005.012j Z +=，对地支路33.030j y =。
2.计算机计算：编写潮流计算程序，要求如下：2.1据给定的潮流计算任务书整理潮流计算的基础数据：节点的分类，线路模型，等值变压器模型，电压等级的归算，标幺值的计算；2.2基础数据的计算机存储：节点数据，支路数据（包括变压器）； 2.3用牛顿-拉夫逊法计算；2.4根据所选潮流计算方法画流程图，划分出功能模块，有数据输入模块，导纳阵形成模块，解线性方程组模块，计算不平衡功率模块，形成雅可比矩阵模块，解修正方程模块，计算线路潮流，网损，PV 节点无功功率和平衡节点功率，数据输出模块； 2.5据上述模块编制程序并上机调试程序，得出潮流计算结果； 2.6源程序及其程序中的符号说明集、程序流图简单系统如下图所示，支路数据如下：41.01.012j z +=，3.013j z =，5.012.014j z +=，40.008.024j z +=01528.01,202,10j y y ==，0192.01,404,10j y y ==，01413.02,404,20j y y ==1.1=k节点数据如下：18.030.01j S --= ，13.055.02j S --= ， 5.03=S ，10.13=U ，o U 005.14∠= 32S 1S 4S 4,10y 1,40y 2,10y 1,20y 4,20y 2,40y 12z 24z 14z k)1(13-k k z kz -1131)节点导纳阵#include <stdio.h> #include <math.h> #include <fstream.h> #include "LF.h"//form node conductance matrixintMakeY( intnB, intnL, Line* sL, double** YG, double** YB ) {inti,j,l;double r,x,d1,g,b,t;for(i=0;i<nB;i++)for(j=0;j<nB;j++){ YG[i][j]=0.0;YB[i][j]=0.0;}for(i=0;i<nL;i++){r=sL[i].R;x=sL[i].X;g=r/(r*r+x*x);b=-x/(r*r+x*x);switch(sL[i].Type){case 1://Linebreak;case 2://Transformerg*=1/sL[i].K;b*=1/sL[i].K;break;}YG[sL[i].NumI][sL[i].NumI]+=g;YG[sL[i].NumJ][sL[i].NumJ]+=g;YG[sL[i].NumI][sL[i].NumJ]-=g;YG[sL[i].NumJ][sL[i].NumI]-=g;YB[sL[i].NumI][sL[i].NumI]+=b+sL[i].B;YB[sL[i].NumJ][sL[i].NumJ]+=b+sL[i].B;YB[sL[i].NumI][sL[i].NumJ]-=b;YB[sL[i].NumJ][sL[i].NumI]-=b; }printf("实部：\n");for(i=0;i<nB;i++){for(j=0;j<nB;j++)printf("%lf\t",YG[i][j]);printf("\n");}printf("虚部：\n");for(i=0;i<nB;i++){for(j=0;j<nB;j++)printf("%lf\t",YB[i][j]);printf("\n");}/* Check the Y matrix */ofstreamfout("out.txt");fout<< "--------------Y Matrix--------------------" <<endl;for(i=0;i<nB;i++){for(j=0;j<nB;j++)fout<< YG[i][j] << "+j" << YB[i][j] << "\t";fout<<endl;}fout.close();return 0;}2）计算功率不平衡量#include <stdio.h>#include <math.h>#include <fstream.h>#include "LF.h"//form delta p and delta qintCalDeltaPQ( intnpv, intnpq, Bus* bus, double** YG, double** YB, int* p_Jtobus, double* deltaf ){intk,i,j;for(k=0;k<npv+npq*2;k++){ i=p_Jtobus[k];if(k<npv){ deltaf[k]=bus[i].GenP-bus[i].LoadP;for(j=0;j<npv+npq+1;j++){deltaf[k]-=bus[i].Volt*bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j ]*sin(bus[i].Phase-bus[j].Phase));}printf("PV节点%d的有功功率是%lf\n",i,deltaf[k]);}if(k<npq+npv&&k>=npv){ deltaf[k]=bus[i].GenP-bus[i].LoadP;for(j=0;j<npv+npq+1;j++){deltaf[k]-=bus[i].Volt*bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j ]*sin(bus[i].Phase-bus[j].Phase));}printf("PQ节点%d的有功功率是%lf\n",i,deltaf[k]);}if(k<npq*2+npv&&k>=npv+npq){deltaf[k]=bus[i].GenQ-bus[i].LoadQ;for(j=0;j<npv+npq+1;j++){deltaf[k]-=bus[i].Volt*bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j ]*cos(bus[i].Phase-bus[j].Phase));}printf("PQ节点%d的无功功率是 %lf\n",i,deltaf[k]);}}return 0;}3）雅各比矩阵的计算/*Purpose: for undergraduate courseTask: Load FlowCopyright @ NCEPU, Liu Chongru*/#include <stdio.h>#include <math.h>#include <fstream.h>#include "LF.h"//form Jacobian matrixintFormJacobian( intnpv, intnpq, Bus* bus, double** YG, double** YB, int* p_Jtobus, double** Jac ){ intnp = npv+npq,j,k,i,m;//TODOdouble a[14],q[14];for(k=0;k<npv+npq*2;k++){i=p_Jtobus[k];a[i]=0;q[i]=0;if(k<np)//H N{for(j=0;j<np+1;j++)if(j!=i){a[i]+=bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j]*cos(bus[i].Phase-bu s[j].Phase));q[i]+=bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j]*sin(bus[i].Phase-bu s[j].Phase));}for(m=0;m<npv+npq*2;m++){j=p_Jtobus[m];if(j!=i){if(m<np)Jac[k][m]=bus[i].Volt*bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j]*cos (bus[i].Phase-bus[j].Phase));//Form HelseJac[k][m]=bus[i].Volt*bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j]*sin (bus[i].Phase-bus[j].Phase));//Form N}else if(j==i){if(m<np)Jac[k][m]=-bus[i].Volt*a[i];//Form HelseJac[k][m]=bus[i].Volt*q[i]+2*bus[i].Volt*bus[i].Volt*YG[i][j];//Form N}}}else{for(j=0;j<np+1;j++)if(j!=i){a[i]+=bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j]*cos(bus[i].Phase-bu s[j].Phase));q[i]+=bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j]*sin(bus[i].Phase-bu s[j].Phase));}for(m=0;m<npv+npq*2;m++){j=p_Jtobus[m];if(j!=i){if(m<np)Jac[k][m]=-bus[i].Volt*bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j]*si n(bus[i].Phase-bus[j].Phase)); //Form JelseJac[k][m]=bus[i].Volt*bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j] *cos(bus[i].Phase-bus[j].Phase)); //Form L }else if(j==i){ if(m<np)Jac[k][m]=bus[i].Volt*q[i];elseJac[k][m]=bus[i].Volt*a[i]-2*bus[i].Volt*bus[i].Volt*YB[i][j];} }}}for(i=0;i<np+npq;i++){for(int j=0;j<np+npq;j++){printf("%d %d %f ",i,j,Jac[i][j]);}printf("\n");}//Output the matrix to check the Jacobian matrixofstreamfout("out.txt",ios::app);fout<< "--------------Jacobian Matrix--------------------" <<endl;for(i=0; i<np+npq;i++ ){for(j=0; j<np+npq; j++ )fout<<Jac[i][j] << "\t";fout<<endl;}fout.close();return 0;}4）线路损耗//8.calculate the power flowdouble* p_Pij, *p_Qij, *p_Pji, *p_Qji;p_Pij = new double[nL];p_Qij = new double[nL];p_Pji = new double[nL];p_Qji = new double[nL];int x1,x2;for( i=0; i<nL; i++ ){ x1=line[i].NumI;x2=line[i].NumJ;if(line[i].Type==1){p_Pij[i]=bus[x1].Volt*bus[x1].Volt*(-YG[x1][x2])-bus[x1].Volt*bus[x2].Volt*((-YG[x1][x2 ])*cos(bus[x1].Phase-bus[x2].Phase)+(-YB[x1][x2])*sin(bus[x1].Phase-bus[x2].Phase));p_Qij[i]=-bus[x1].Volt*bus[x1].Volt*(line[i].B+(-YB[x1][x2]))-bus[x1].Volt*bus[x2]. Volt*((-YG[x1][x2])*sin(bus[x1].Phase-bus[x2].Phase)-(-YB[x1][x2])*cos(bus[x1].Phase-bu s[x2].Phase));p_Pji[i]=bus[x2].Volt*bus[x2].Volt*(-YG[x2][x1])-bus[x2].Volt*bus[x1].Volt*((-YG[x2][x1 ])*cos(bus[x2].Phase-bus[x1].Phase)+(-YB[x2][x1])*sin(bus[x2].Phase-bus[x1].Phase));p_Qji[i]=-bus[x2].Volt*bus[x2].Volt*(line[i].B+(-YB[x2][x1]))-bus[x2].Volt*bus[x1]. Volt*((-YG[x2][x1])*sin(bus[x2].Phase-bus[x1].Phase)-(-YB[x2][x1])*cos(bus[x2].Phase-bu s[x1].Phase));}else{p_Pij[i]=bus[x1].Volt*bus[x1].Volt*(-YG[x1][x2])/line[i].K-bus[x1].Volt*bus[x2].Vol t*((-YG[x1][x2])*cos(bus[x1].Phase-bus[x2].Phase)+(-YB[x1][x2])*sin(bus[x1].Phase-bus[x 2].Phase));p_Qij[i]=-bus[x1].Volt*bus[x1].Volt*((-YB[x1][x2])/line[i].K+line[i].B)-bus[x1].Volt*bu s[x2].Volt*((-YG[x1][x2])*sin(bus[x1].Phase-bus[x2].Phase)-(-YB[x1][x2])*cos(bus[x1].Ph ase-bus[x2].Phase));p_Pji[i]=bus[x2].Volt*bus[x2].Volt*(-YG[x2][x1]*line[i].K)-bus[x2].Volt*bus[x1].Volt*(( -YG[x2][x1])*cos(bus[x2].Phase-bus[x1].Phase)+(-YB[x2][x1])*sin(bus[x2].Phase-bus[x1].P hase));p_Qji[i]=-bus[x2].Volt*bus[x2].Volt*((-YB[x2][x1])*line[i].K+line[i].B)-bus[x2].Volt*bu s[x1].Volt*((-YG[x2][x1])*sin(bus[x2].Phase-bus[x1].Phase)-(-YB[x2][x1])*cos(bus[x2].Ph ase-bus[x1].Phase));}}//p and q of PH bus and PV busint s=0;double p[9],q[9],Ps[9],Qs[9],PS=0,QS=0;for( i=0; i<nB; i++ ){p[i]=0;q[i]=0;for(int j=0; j<nB; j++ ){p[i]+=(bus[j].Volt*(YG[i][j])*cos(bus[j].Phase)-bus[j].Volt*(YB[i][j])*sin(bus[j].Phas e));q[i]-=(bus[j].Volt*(YG[i][j])*sin(bus[j].Phase)+bus[j].Volt*(YB[i][j])*cos(bus[j].Phase ));}Ps[i]=bus[i].Volt*cos(bus[i].Phase)*p[i]-bus[i].Volt*sin(bus[i].Phase)*q[i]; Qs[i]=bus[i].Volt*cos(bus[i].Phase)*q[i]+bus[i].Volt*sin(bus[i].Phase)*p[i];} for(i=0;i<nB;i++){ PS+=Ps[i];QS+=Qs[i];}printf("PS=%7.7f,QS=%7.7f\n",PS,QS);}//lossdoublePloss=0, Qloss=0;for( i=0; i<nB; i++ ){Ploss+=p_Pij[i]+p_Pji[i];Qloss+=p_Qij[i]+p_Qji[i];}5）程序流图如下6）得到的数据（out.txt ）输入原始数据形成节点导纳矩阵设非平衡节点电压初值)0(i e ;)0(i f令迭代次数k=0对PQ 节点计算)()(.k i k i Q P ∆∆对PV 节点计算)(2)(,k i k iU P ∆∆令节点号i=1计算雅客比矩阵各元素)()()()()()(,,,,,k ijk ij k ij k ij k ij k ij S R L J N H 增加节点号i=i+1解修正方程，由)()(.k i k i Q P ∆∆)(2k i U ∆及雅客比矩阵用高斯法求各节点的电压增量)()(,k ik if e∆∆求出max)(max)(,k k f e ∆∆迭代是否收敛ε<∆∆max)(max)(,k k f e 雅客比矩阵是否形成，i>n?计算节点的新电压)()()1(k i k i k i e e e ∆+=+ )()()1(k i k i k i f f f ∆+=+增加迭代次数k=k+1计算平衡节点的功率及线路功率停止启动--------------Y Matrix--------------------0+j-17.3611 0+j0 0+j0 0+j17.3611 0+j0 0+j0 0+j0 0+j0 0+j00+j0 0+j-16 0+j0 0+j0 0+j0 0+j0 0+j16 0+j0 0+j00+j0 0+j0 0+j-17.0648 0+j0 0+j0 0+j0 0+j0 0+j0 0+j17.06480+j17.3611 0+j0 0+j0 3.30738+j-39.3089 -1.36519+j11.6041 -1.94219+j10.5107 0+j0 0+j0 0+j00+j0 0+j0 0+j0 -1.36519+j11.6041 2.55279+j-17.3382 0+j0 -1.1876+j5.97513 0+j0 0+j00+j0 0+j0 0+j0 -1.94219+j10.5107 0+j0 3.2242+j-15.8409 0+j0 0+j0 -1.28201+j5.588240+j0 0+j16 0+j0 0+j0 -1.1876+j5.97513 0+j0 2.80473+j-35.4456 -1.61712+j13.698 0+j00+j0 0+j0 0+j0 0+j0 0+j0 0+j0 -1.61712+j13.698 2.77221+j-23.3032 -1.15509+j9.784270+j0 0+j0 0+j17.0648 0+j0 0+j0 -1.28201+j5.58824 0+j0 -1.15509+j9.78427 2.4371+j-32.1539--------------Jacobian Matrix--------------------16.4 0 0 0 0 -16.4 0 0 0 0 0 0 0 00 17.4915 0 0 0 0 0 -17.4915 0 0 0 0 0 00 0 40.1703 -11.6041 -10.5107 0 0 0 3.30738 -1.36519 -1.94219 0 0 00 0 -11.6041 17.5792 0 -5.97513 0 0 -1.36519 2.55279 0 -1.1876 0 00 0 -10.5107 0 16.0989 0 0 -5.58824 -1.94219 0 3.2242 0 0 -1.28201-16.4 0 0 -5.97513 0 36.0731 -13.698 0 0 -1.1876 0 2.80473 -1.61712 00 0 0 0 0 -13.698 23.4822 -9.78427 0 0 0 -1.61712 2.77221 -1.155090 -17.4915 0 0 -5.58824 0 -9.78427 32.864 0 0 -1.28201 0 -1.15509 2.43710 0 -3.30738 1.36519 1.94219 0 0 0 38.4474 -11.6041 -10.5107 0 0 00 0 1.36519 -2.55279 0 1.1876 0 0 -11.6041 17.0972 0 -5.97513 0 00 0 1.94219 0 -3.2242 0 0 1.28201 -10.5107 0 15.5829 0 0 -5.588240 0 0 1.1876 0 -2.80473 1.61712 0 0 -5.97513 0 34.8181 -13.698 00 0 0 0 0 1.61712 -2.77221 1.15509 0 0 0 -13.698 23.1242 -9.784270 0 0 0 1.28201 0 1.15509 -2.4371 0 0 -5.58824 0 -9.78427 31.4437--------------Jacobian Matrix--------------------16.9269 0 0 0 0 -16.9269 0 0 0 0 0 1.68793 0 00 18.1691 0 0 0 0 0 -18.1691 0 0 0 0 0 0.8836270 0 41.9297 -12.1301 -11.1536 0 0 0 3.54272 -1.0628 -1.76646 0 0 00 0 -12.0455 18.0609 0 -6.01539 0 0 -1.78138 1.30819 0 -2.10262 0 00 0 -11.0484 0 16.8144 0 0 -5.76607 -2.33608 0 2.42598 0 0 -1.97778-16.9269 0 0 -6.36224 0 37.9476 -14.6585 0 0 -0.357534 0 3.05959 -0.930027 00 0 0 0 0 -14.4721 24.8873 -10.4152 0 0 0 -2.509 1.86088 -1.473890 -18.1691 0 0 -6.05157 0 -10.4721 34.6928 0 0 -0.733327 0 -0.9919732.66270 0 -3.52149 1.0628 1.76646 0 0 0 42.0299 -12.1301 -11.1536 0 0 00 0 1.78138 -3.884 0 2.10262 0 0 -12.0455 17.2037 0 -6.01539 0 00 0 2.33608 0 -4.31386 0 0 1.97778 -11.0484 0 16.2993 0 0 -5.766071.68793 0 0 0.357534 0 -2.97549 0.930027 0 0 -6.36224 0 38.3226 -14.6585 00 0 0 0 0 2.509 -3.98289 1.47389 0 0 0 -14.4721 24.2355 -10.41520 0.883627 0 0 0.733327 0 0.991973 -2.60893 0 0 -6.05157 0 -10.4721 34.8585 --------------Jacobian Matrix--------------------16.7457 0 0 0 0 -16.7457 0 0 0 0 0 1.63043 0 00 18.0388 0 0 0 0 0 -18.0388 0 0 0 0 0 0.8501960 0 41.3695 -11.8919 -10.9686 0 0 0 3.48069 -1.02775 -1.73712 0 0 00 0 -11.8057 17.6918 0 -5.8861 0 0 -1.7602 1.28091 0 -2.0217 0 00 0 -10.8651 0 16.5476 0 0 -5.68251 -2.29737 0 2.40655 0 0 -1.91027-16.7457 0 0 -6.21183 0 37.3041 -14.3465 0 0 -0.382862 0 2.95313 -0.937485 00 0 0 0 0 -14.1704 24.4052 -10.2348 0 0 0 -2.42909 1.86079 -1.433530 -18.0388 0 0 -5.94693 0 -10.2872 34.273 0 0 -0.757656 0 -0.9892082.598470 0 -3.48089 1.02775 1.73712 0 0 0 41.3703 -11.8919 -10.9686 0 0 00 0 1.7602 -3.78189 0 2.0217 0 0 -11.8057 16.6941 0 -5.8861 0 00 0 2.29737 0 -4.20764 0 0 1.91027 -10.8651 0 15.9488 0 0 -5.682511.63043 0 0 0.382862 0 -2.95077 0.937485 0 0 -6.21183 0 37.3083 -14.3465 00 0 0 0 0 2.42909 -3.86262 1.43353 0 0 0 -14.1704 23.7059 -10.23480 0.850196 0 0 0.757656 0 0.989208 -2.59706 0 0 -5.94693 0 -10.2872 34.2743 --------------Jacobian Matrix--------------------16.7435 0 0 0 0 -16.7435 0 0 0 0 0 1.63 0 00 18.0374 0 0 0 0 0 -18.0374 0 0 0 0 0 0.850 0 41.3625 -11.8888 -10.9664 0 0 0 3.48016 -1.02713 -1.73662 0 0 00 0 -11.8026 17.6871 0 -5.8845 0 0 -1.76008 1.28053 0 -2.02045 0 00 0 -10.8628 0 16.5444 0 0 -5.68158 -2.29703 0 2.40632 0 0 -1.90929-16.7435 0 0 -6.20987 0 37.296 -14.3426 0 0 -0.383399 0 2.95114 -0.937742 00 0 0 0 0 -14.1667 24.3994 -10.2326 0 0 0 -2.42794 1.86097 -1.433020 -18.0374 0 0 -5.94567 0 -10.285 34.2681 0 0 -0.758139 0 -0.9892032.597340 0 -3.48016 1.02713 1.73662 0 0 0 41.3625 -11.8888 -10.9664 0 0 00 0 1.76008 -3.78053 0 2.02045 0 0 -11.8026 16.6871 0 -5.8845 0 00 0 2.29703 0 -4.20632 0 0 1.90929 -10.8628 0 15.9444 0 0 -5.681581.63 0 0 0.383399 0 -2.95114 0.937742 0 0 -6.20987 0 37.296 -14.3426 00 0 0 0 0 2.42794 -3.86097 1.43302 0 0 0 -14.1667 23.6994 -10.23260 0.85 0 0 0.758139 0 0.989203 -2.59734 0 0 -5.94567 0 -10.285 34.2681--------------iteration---------------iteration = 4--------------voltage magnitude and angle--------------------1.04 0 01.025 0.161967 9.280011.025 0.0814153 4.664761.02579 -0.0386902 -2.216790.995631 -0.0696178 -3.988811.01265 -0.0643572 -3.68741.02577 0.064921 3.71971.01588 0.0126979 0.7275371.03235 0.0343257 1.96672-------------bus P and Q------------------1 0.71641 0.2704592 1.63 0.06653663 0.85 -0.1085974 0 05 -1.25 -0.56 -0.9 -0.37 0 08 -1 -0.359 0 0----------------line flow--------------------NUM------i-------j-----------begin------------end--------------1 4 1 -0.71641+j-0.239232 0.71641+j0.270462 7 2 -1.63+j0.0917816 1.63+j0.06653653 9 3 -0.85+j0.149553 0.85+j-0.1085974 7 8 0.763799+j-0.00797398 -0.759046+j-0.1070415 9 8 0.241834+j0.0311946 -0.240954+j-0.2429586 7 5 0.866201+j-0.0838079 -0.843202+j-0.1131287 9 6 0.608166+j-0.180748 -0.594627+j-0.1345668 5 4 -0.406798+j-0.386872 0.409373+j0.2289319 6 4 -0.305372+j-0.165433 0.307037+j0.0102993--------------------Ploss and Qloss-----------------Ploss = 0.0471901 Qloss = -0.9574833.思考题3.1潮流计算的方法有哪些？各有何特点？答：潮流计算分为手算和机算两大类，其中机算又有高斯-赛德尔迭代法、牛顿-拉夫逊迭代法、P-Q 分解法等算法。
华电潮流上机课程设计报告(程序最终无误版)
《电力系统潮流上机》课程设计报告院系：电气与电子工程学院班级：电气1108学号： 1111550112学生：龙日尚指导教师：宝柱设计周数：两周成绩：日期：2014年1月10日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细容见附录）1．手算2．计算机计算3．思考题三、课程设计总结或结论（详细容见附录）四、参考文献1.《电力系统计算：电子数字计算机的应用》，交通大学等合编。
北京：水利电力；2.《现代电力系统分析》，王锡凡主编，科学；3.《电力系统稳态分析》，珩，中国电力，1995年，第三版；附录（设计流程图、程序、表格、数据等）4.机算潮流程序及结果// dierti.cpp : 定义控制台应用程序的入口点。
//#include "stdafx.h"struct Line //线路结构体{int Num,NumI,NumJ; //线路号左节点名右节点名float R,X,B,K; //电阻电抗电纳变比（K等于1为普通支路，不等于1为变压器支路的变比） };struct Bus //节点结构体{int Num ;float Volt,Phase,GenP,GenQ,LoadP,LoadQ;int Type;};#include"stdio.h"#include"string.h"#include"math.h"#include"stdlib.h"#define NBUS 4#define NLINE 4/* Global variables */int nL,nB,nVA,nSH;float X[NBUS];int L;double def[2*NBUS];double mn[50];void Gauss(double a[50][50],double b[50], int n) /*定义高斯法 */{int JS[50];int i,j,k;float d,t,x[50];FILE *fp;int L=1;for(i=0;i<50;i++) JS[i]=0;for(k=0;k<n;k++){d=0.0;for(j=k;j<n;j++)if(fabs(a[k][j])>d){ /*在一行中找到一个最大值赋值d,并用JS[K]记住这个最大值所在的列号*/ d=fabs(a[k][j]);JS[k]=j;}if(fabs(d)<0.000001) /*如果d的数值太小，做为被除数将带来很大的误差 */L=0;else {if(JS[k]!=k)for(i=0;i<n;i++){t=a[i][k];a[i][k]=a[i][JS[k]]; /*进行列交换，让最大值始终在对角元上*/a[i][JS[k]]=t;}}if(L==0)break;for(j=k+1;j<n;j++)a[k][j]=a[k][j]/a[k][k]; /*对角元上的元素消为1*/b[k]=b[k]/a[k][k];for(i=k+1;i<n;i++){for(j=k+1;j<n;j++)a[i][j]=a[i][j]-a[i][k]*a[k][j]; /*使下三角阵的元素为0*/b[i]=b[i]-a[i][k]*b[k];}}if(fabs(a[n-1][n-1])>0.00001){ /*用追赶法，解方程组，求未知数x*/ x[n-1]=b[n-1];for(i=n-2;i>=0;i--){t=0.0;for(j=i+1;j<n;j++)t=t+a[i][j]*x[j];x[i]=(b[i]-t);}}if((fp=fopen("gauss.txt","w"))==NULL) /*将结果写到TXT文件中*/{printf("err");exit(0);}for(i=0;i<n;i++){fprintf(fp,"%lf",x[i]);mn[i]=x[i];fprintf(fp,"\n");}fclose(fp);if(fp!=NULL) fclose(fp);}int _tmain(int argc, _TCHAR* argv[]){FILE *fp;FILE *fpout;int i,j,k,l,h,n,v;int i1,i2,i3,kp,kq;float d1,d2,d3,d4,d5,d6,r,x,g,b,tt,LL,e,ps,qs,shsh,m;struct Line sL[NLINE];struct Bus sB[NBUS];float YG[NBUS+1][NBUS+1],YB[NBUS+1][NBUS+1];double u[50][2];i1=i2=i3=0;d1=d2=d3=d4=d5=d6=ps=qs=0.0;for(i=0;i<NBUS;i++)if((fp=fopen("in.txt","r"))==NULL){ printf("Can not open the file named 'in.txt' \n");exit(0);}fscanf(fp,"%d,%d,%d",&nB,&nL,&nSH);for(i=0;i<nB;i++){sB[i].Num=sB[i].Type=0;sB[i].Volt=1.0;sB[i].Phase=sB[i].GenP=sB[i].GenQ=sB[i].LoadP=sB[i].LoadQ=0.0;fscanf(fp,"%d,%f,%f,%f,%f,%f,%f,%d",&i1,&d1,&d2,&d3,&d4,&d5,&d6,&i2);sB[i].Num=i1;sB[i].Volt=d1;sB[i].Phase=d2;sB[i].GenP=d3;sB[i].GenQ=d4;sB[i].LoadP=d5,sB[i].LoadQ=d6;sB[i].T ype=i2;};for(i=0;i<nL;i++){sL[i].Num=sL[i].NumI=sL[i].NumJ=0;sL[i].R=sL[i].X=sL[i].B=0.0;sL[i].K=1.0;fscanf(fp,"%2d %3d %3d %f %f %f %f",&i1,&i2,&i3,&d1,&d2,&d3,&d4);sL[i].Num=i1;sL[i].NumI=i2;sL[i].NumJ=i3;sL[i].R=d1;sL[i].X=d2;sL[i].B=d3;sL[i].K=d4;}if(fp!=NULL) fclose(fp);/*Make Y Matrix*/for(i=1;i<nB+1;i++)for(j=1;j<nB+1;j++){YG[i][j]=0.0;YB[i][j]=0.0;};for(l=0; l<nL; l++){i=sL[l].NumI;j=sL[l].NumJ;r=sL[l].R;x=sL[l].X;d1=r*r+x*x;g=r/d1;b=-x/d1;m=sL[l].K;if(fabs(sL[l].K-1.0)<0.000001) //普通支路 {YG[i][i]=YG[i][i]+g;YG[j][j]=YG[j][j]+g;YB[i][i]=YB[i][i]+b+sL[l].B;YB[j][j]=YB[j][j]+b+sL[l].B;YG[i][j]=YG[i][j]-g;YG[j][i]=YG[j][i]-g;YB[i][j]=YB[i][j]-b;YB[j][i]=YB[j][i]-b;}else //变压器支路{YG[i][i]=YG[i][i]+g/m+g*(m-1)/m;YG[j][j]=YG[j][j]+g/m+g*(1-m)/m/m;YB[i][i]=YB[i][i]+b/m+b*(m-1)/m;YB[j][j]=YB[j][j]+b/m+b*(1-m)/m/m;YG[i][j]=YG[i][j]-g/m;YG[j][i]=YG[j][i]-g/m;YB[i][j]=YB[i][j]-b/m;YB[j][i]=YB[j][i]-b/m; }}/* Check the Y matrix */if((fp=fopen("GGBB.txt","w"))==NULL){printf("Can not open the file named 'GGBB.txt' \n");exit(0);}fprintf(fp,"---Y Matrix---\n");for(i=1;i<nB+1;i++)for(j=1;j<nB+1;j++)if(fabs(YB[i][j]-0.0)>0.000001) fprintf(fp,"Y(%3d,%-3d)=(%10.5f,%10.5f)\n",i,j,YG[i][j],YB[i][j]);if(fp!=NULL) fclose(fp);/* 节点电压附初值 */for(i=1;i<nB+1;i++){if(sB[i-1].Type==0){u[i][0]=0.0;u[i][1]=1.0;}else if(sB[i-1].Type==1){u[i][1]=sB[i-1].Volt;u[i][0]=0.0;}else if(sB[i-1].Type==2){u[i][1]=sB[i-1].Volt;u[i][0]= sB[i-1].Phase;}for(v=1;;v++)/* 迭代次数可以无限大 */{/* 节点电压附初值 */printf("迭代第%d次赋予的电压初值为e+jf:\n",v); for(i=1;i<nB+1;i++)printf("%lf,%lf\n",u[i][1],u[i][0]);printf("\n");printf("\n");/* 求偏移量 */double P_P[10];double P_Q[10];double P_UU[10];for(i=1;i<nB+1;i++){if(sB[i-1].Type==2){P_P[i]=0.0;P_Q[i]=0.0;P_UU[i]=1.05;}if(sB[i-1].Type==0){double tempP=0.0;double tempQ=0.0;for(j=1;j<nB+1;j++){tempP+=YG[i][j]*u[j][1]-YB[i][j]*u[j][0]; tempQ+=YG[i][j]*u[j][0]+YB[i][j]*u[j][1];P_P[i]=(sB[i-1].GenP-sB[i-1].LoadP)-tempP*u[i][1]-tempQ*u[i][0]; P_Q[i]=(sB[i-1].GenQ-sB[i-1].LoadQ)-tempP*u[i][0]+tempQ*u[i][1]; P_UU[i]=0.0;}if(sB[i-1].Type==1){double tempP=0.0;double tempQ=0.0;for(j=1;j<nB+1;j++){tempP+=YG[i][j]*u[j][1]-YB[i][j]*u[j][0];tempQ+=YG[i][j]*u[j][0]+YB[i][j]*u[j][1];P_P[i]=(sB[i-1].GenP-sB[i-1].LoadP)-tempP*u[i][1]-tempQ*u[i][0]; }P_UU[i]=sB[i-1].Volt*sB[i-1].Volt-u[i][1]*u[i][1]-u[i][0]*u[i][0]; P_Q[i]=0.0;}}/* 偏移量阵 */double P_PQ[6];int a=0;for(i=1;i<3;i++){P_PQ[a]=P_P[i];a=a+2;}a=1;for(i=1;i<3;i++){P_PQ[a]=P_Q[i];a=a+2;}P_PQ[4]=P_P[3];P_PQ[5]=P_UU[3];printf("迭代第%d次的偏移量为：\n",v);for(i=0;i<6;i++){printf("%f",P_PQ[i]);printf("\n");}printf("\n");printf("\n");/* 雅可比矩阵 */double H[6][6],N[6][6],J[6][6],L[6][6],R[6][6],S[6][6],aa[6],bb[6]; for(i=1;i<5;i++){ if(fabs(sB[i-1].Type-2.0)<0.000001)continue;else{for(j=1;j<5;j++)if(i!=j){H[i][j]=-YB[i][j]*u[i][1]+YG[i][j]*u[i][0];N[i][j]=YG[i][j]*u[i][1]+YB[i][j]*u[i][0];J[i][j]=-N[i][j];L[i][j]=H[i][j];R[i][j]=0;S[i][j]=0;}else{aa[i]=bb[i]=0.0;for(n=1;n<5;n++)aa[i]+=YG[i][n]*u[n][1]-YB[i][n]*u[n][0];bb[i]+=YG[i][n]*u[n][0]+YB[i][n]*u[n][1];}H[i][i]=-YB[i][i]*u[i][1]+YG[i][i]*u[i][0]+bb[i]; N[i][i]=YG[i][i]*u[i][1]+YB[i][i]*u[i][0]+aa[i]; J[i][i]=-YG[i][i]*u[i][1]-YB[i][i]*u[i][0]+aa[i]; L[i][i]=YG[i][i]*u[i][0]-YB[i][i]*u[i][1]-bb[i]; R[i][i]=2*u[i][0];S[i][i]=2*u[i][1];}}}double ss[50][50];for(i=0;i<6;i++)for(j=0;j<6;j++)ss[i][j]=0.0;for(i=1;i<3;i++)for(j=1;j<4;j++){ss[2*i-2][2*j-2]=H[i][j];ss[2*i-2][2*j-1]=N[i][j];ss[2*i-1][2*j-2]=J[i][j];ss[2*i-1][2*j-1]=L[i][j];}i=3;for(j=1;j<4;j++){ss[2*i-2][2*j-2]=H[i][j];ss[2*i-2][2*j-1]=N[i][j];ss[2*i-1][2*j-2]=R[i][j];ss[2*i-1][2*j-1]=S[i][j];}printf("迭代第%d次的雅可比矩阵为：\n",v);for(i=0;i<6;i++){for(j=0;j<6;j++)printf("%10f",ss[i][j]);printf("\n");}printf("\n");printf("\n");Gauss(ss,P_PQ,6);for(i=1;i<nB;i++){u[i][0]=u[i][0]+mn[2*(i-1)];u[i][1]=u[i][1]+mn[2*i-1];}double max;max=fabs(P_PQ[0]);for(i=0;i<=5;i++)if (max<fabs(P_PQ[i]))max=fabs(P_PQ[i]);if(fabs(max)<0.0001){printf("满足精度要求，迭代终止，迭代次数为%d\n",v); printf("\n");printf("\n");}/* 叠代循环的括号 */printf("最终求得的节点电压值为e+jf:\n");for(i=1;i<nB+1;i++)printf("%lf,%lf\n",u[i][1],u[i][0]);printf("\n");printf("\n");double uu[5],Phase[5];for(i=1;i<nB+1;i++){uu[i]=sqrt(u[i][1]*u[i][1]+u[i][0]*u[i][0]); Phase[i]=atan(u[i][0]/u[i][1]);}for(i=1;i<nB+1;i++)printf("%lf,%lf\n",uu[i],Phase[i]);*计算线路功率和平衡节点 PV节点功率*/double P[5],Q[5];double tempP=0.0;double tempQ=0.0;for(i=1;i<nB+1;i++){for(j=1;j<nB+1;j++){tempP+=YG[i][j]*u[j][1]-YB[i][j]*u[j][0]; tempQ+=YG[i][j]*u[j][0]+YB[i][j]*u[j][1];}P[i]=tempP*u[i][1]+tempQ*u[i][0];Q[i]=tempP*u[i][0]-tempQ*u[i][1];tempP=0.0;for(i=1;i<nB+1;i++)printf("节点%d注入功率为%lf,%lf\n",i,P[i],Q[i]);/* 支路功率 */double V[4][2];for(i=1;i<5;i++)for(j=0;j<3;j++)V[i][j]=u[i][j];double sP[5][5],sQ[5][5];double dsq,dsp,dp,sumgen;for(i=1;i<NBUS+1;i++){for(j=1;j<NBUS+1;j++){sP[i][j]=0.0;sQ[i][j]=0.0;}}for(l=0; l<nL; l++){i=sL[l].NumI;j=sL[l].NumJ;r=sL[l].R;x=sL[l].X;d1=r*r+x*x;g=r/d1;if(fabs(sL[l].K-1.0)<0.000001){/*Normal lines or transformers*/sP[i][j]=V[i][1]*V[i][1]*g-V[i][1]*V[j][1]*(g*cos(V[i][0]-V[j][0])+b*sin(V[i][0]-V[j][0]));sQ[i][j]=-(V[i][1]*V[i][1]*sL[l].B+V[i][1]*V[i][1]*b+V[i][1]*V[j][1]*(g*sin(V[i][0]-V[j][0])-b*cos(V[i ][0]-V[j][0])));sP[j][i]=V[j][1]*V[j][1]*g-V[i][1]*V[j][1]*(g*cos(V[j][0]-V[i][0])+b*sin(V[j][0]-V[i][0]));sQ[j][i]=-(V[j][1]*V[j][1]*sL[l].B+V[j][1]*V[j][1]*b+V[i][1]*V[j][1]*(g*sin(V[j][0]-V[i][0])-b*cos(V[j ][0]-V[i][0])));}else{/*abnormal transformer ratio*/sP[i][j]=V[i][1]*V[i][1]*g/sL[l].B/sL[l].B-V[i][1]*V[j][1]*(g*cos(V[i][0]-V[j][0])/sL[l].B+b*sin(V[i][ 0]-V[j][0])/sL[l].B);sQ[i][j]=-(V[i][1]*V[i][1]*b/sL[l].B/sL[l].B+V[i][1]*V[j][1]*(g*sin(V[i][0]-V[j][0])/sL[l].B-b*cos(V[i ][0]-V[j][0])/sL[l].B));sP[j][i]=V[j][1]*V[j][1]*g-V[i][1]*V[j][1]*(g*cos(V[j][0]-V[i][0])/sL[l].B+b*sin(V[j][0]-V[i][0])/sL[l ].B);sQ[j][i]=-(V[i][1]*V[i][1]*b+V[i][1]*V[j][1]*(g*sin(V[j][0]-V[i][0])/sL[l].B-b*cos(V[j][0]-V[i][0])/sL [l].B));}}/* 输电效率 */dsp=P[4];dsq=Q[4];for(i=0;i<NBUS;i++){dsp+=sB[i].GenP-sB[i].LoadP;dsq+=sB[i].GenQ-sB[i].LoadQ;sumgen+=sB[i].GenP;}dp=dsp/sumgen*100;/* 输出功率情况 */if((fp=fopen("功率情况.txt","w"))==NULL){printf("Can not open the file named '功率情况.txt' \n");exit(0);}fprintf(fp,"---功率情况---\n");fprintf(fp,"平衡节点功率S=%10.5f+ j%10.5f\n",P[4],Q[4]);for(i=1;i<NBUS+1;i++)for(j=1;j<NBUS+1;j++)if(fabs(sP[i][j]-0.0)>0.000001)fprintf(fp,"S(%3d,%-3d)=(%10.5f,j%10.5f)\n",i,j,sP[i][j],sQ[i][j]); fprintf(fp,"网损为%10.5f+j%10.3f,输电效率为%10.3f\n",dsp,dsq,100-dp);if(fp!=NULL) fclose(fp);return 0;}1．导纳阵Y( 1,1 )=( 1.01534, -8.19201) Y( 1,2 )=( -0.56148, 2.30208) Y( 1,3 )=( 0.00000, 3.66667) Y( 1,4 )=( -0.45386, 1.89107) Y( 2,1 )=( -0.56148, 2.30208) Y( 2,2 )=( 1.04225, -4.67651) Y( 2,4 )=( -0.48077, 2.40385) Y( 3,1 )=( 0.00000, 3.66667) Y( 3,3 )=( 0.00000, -3.33333) Y( 4,1 )=( -0.45386, 1.89107) Y( 4,2 )=( -0.48077, 2.40385) Y( 4,4 )=( 0.93463, -4.26159)2．设定电压初值01.1;01;01)0(3)0(3)0(2)0(2)0(1)0(1j jf e j jf e j jf e +=++=++=+ 3．计算功率和电压偏移；27731.0])()([41)0(11)0(1)0(141)0(11)0(1)0(11)0(11)0(1-=++--=-=∆∑∑==jj j j jj jj s s e B f G ff B e G e P P P P0.05097])()([41)0(11)0(1)0(141)0(11)0(1)0(11)0(11)0(1-=----=-=∆∑∑==j j j j j j jj s s e B f G ef B eG f Q Q Q Q同理可算出52596.0)0(22)0(2-=-=∆P P P s ,0196.0)0(22)0(2=-=∆Q Q Q s5.0)0(33)0(3=-=∆P P P s ,0.02)0(3232)0(3=-=∆U U U s4．根据求的第一次迭代时雅可比矩阵各元素的公式计算雅可比矩阵各个元素的具体值：⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤-----------⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡20000.200000066667.300003334.40052691.406629.130208.256148.00001821.182612.456148.030208.266667.3030208.256148.006298.803803.1066667.356148.030208.299265.032104.8 5．求高斯计算后的修正量：⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∆∆∆∆∆∆=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∆∆∆∆∆∆-0.0000000.1276520.023566-0.108546-0.006511-0.007919-2)0(3)0(3)0(2)0(2)0(1)0(11)0()0(3)0(3)0(2)0(2)0(1)0(1U P Q P Q P e f e f e f J 6．计算各节点电压的一次近似值：12765.010855.000792.010000.197643.099349.0)0(3)0(3)1(3)0(2)0(2)1(2)0(1)0(1)1(1)0(3)0(3)1(3)0(2)0(2)1(2)0(1)0(1)1(1=∆+=-=∆+=-=∆+==∆+==∆+==∆+=f f f f f f fffe e e e e e e e e返回第三步重新迭代，并校验收敛与否，令410-=ε。