华电潮流上机课程设计报告+程序.
- 格式:doc
- 大小:117.50 KB
- 文档页数:20
课程设计报告( 2011—2012年度第一学期)名称:电力系统潮流上机院系:电气与电子工程学院班级:学号:学生:指导教师:设计周数:两周成绩:日期: 2011年12月19日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力,掌握计算机潮流计算的相关知识二、设计正文(详细容见附录)1.手算2.计算机计算3.思考题三、课程设计总结或结论潮流计算是研究电力系统稳态运行的一种基本计算,最初求解电力系统潮流时大多使用手算,但随着电力系统结构的日趋复杂,计算量也越来越大。
复杂电力系统潮流计算中,由于节点数量巨大,所形成的修正方程已经无法通过手算方式解决,尤其是需要迭代次数较多时,手算所需要的时间太长,计算机潮流计算无疑为解决这一问题提供了极大的便利。
计算机潮流计算可以迅速解决复杂网络的潮流计算问题,这是由于无论系统的复杂程度如何,其节点与支路的类型是固定的,所以只需要输入节点与支路的数据,就可以解决任何一个复杂网络的潮流计算问题。
即只需要一次编程,就可以基本上解决所有复杂网络的计算。
需要注意的是,在使用牛顿—拉弗逊发计算潮流时,对于初值要选择比较接近它们的精确解,否则迭代过程可能不收敛。
潮流计算C语言程序编程过程中需要注意的是,C语言无法实现复数运算,需要将得到的值的实部与虚部分开储存并计算。
这个过程复杂并且容易出错,编写程序是需要注意。
另外需要注意的一点是:C语言数组的编号是从零开始的,在程序编写过程中应注意下标的对应。
通过这一次的电力系统潮流计算编程,我不仅对C语言的编程有了更深刻的理解,也对《电力系统分析》这门课程进行了查漏补缺和巩固,对电力系统的运行也有了更加深入的了解,受益匪浅。
四、参考文献1.《电力系统计算:电子数字计算机的应用》,交通大学等合编。
:水利电力;2.《现代电力系统分析》,王锡凡主编,科学;3.《电力系统稳态分析》,珩,中国电力,2007年,第三版;附录(包括:1.手算过程及结果;2.计算机计算流程图、表格、数据;3.思考题答案)附录一、手算过程及其结果附录二、计算机计算流程图、表格、数据编写潮流计算程序简单系统如下图所示,支路数据如下: 支路14,27,39为变压器支路,参数为100.1,058.0114==K X ,050.1,063.0227==K X 100.1,059.0339==K X其余支路为线路支路,参数为075.02/,072.0019.07878=+=B j Z , 105.02/,101.0012.08989=+=B j Z 153.02/,161.0032.05757=+=B j Z 179.02/,170.0039.06969=+=B j Z 088.02/,085.0010.04545=+=B j Z 079.02/,092.0017.04646=+=B j Z节点数据如下:o U 004.11∠=•025.1,63.122==U P ,025.1,85.033==U P5.025.15j S --= ,3.09.06j S --= ,35.00.18j S --=节点数据线路数据节点导纳矩阵Y B迭代次数:Count_Num=3各节点电压:平衡节点功率:节点编号1: 0.72887-j0.10579 PV节点无功功率:节点编号2: j0.30332PV节点无功功率:节点编号3: j0.17199线路功率及损耗:线路总损耗:ΔS= 0.05887-j0.78048附录三、思考题1、潮流计算的方法有哪些?各有何特点?答:潮流计算分为简单电力网络的手算和复杂电力网络的机算两大类,其中机算又有高斯-赛德尔迭代法,牛顿-拉夫逊迭代法和P-Q分解法。
电力系统潮流上机计算实验报告11.手算过程已知:节点1:PQ 节点,节点, s(1)= s(1)= -0.5000-j0.3500 节点2:PV 节点,节点, p(2)=0.4000 v(2)=1.0500 p(2)=0.4000 v(2)=1.0500 节点3:平衡节点,:平衡节点,U(3)=1.0000U(3)=1.0000U(3)=1.0000∠∠0.0000 网络的连接图:0.0500+j0.2000 1 0.0500+j0.2000231)计算节点导纳矩阵由2000.00500.012j Z +=Þ71.418.112j y -=; 2000.00500.013j Z +=Þ71.418.113j y -=;\导纳矩阵中的各元素:42.936.271.418.171.418.1131211j j j y y Y -=-+-=+=;71.418.11212jy Y +-=-=; 71.418.11313j y Y +-=-=; =21Y71.418.11212j y Y +-=-=; 71.418.12122j y Y -==; 002323j y Y +=-=;=31Y 71.418.11313j y Y +-=-=; =32Y 002323j y Y +=-=; 71.418.13133j y Y -==;\形成导纳矩阵B Y :úúúûùêêêëé-++-+-+-+-+--=71.418.10071.418.10071.418.171.418.171.418.171.418.142.936.2j j j j j j j j j Y B 2)计算各PQ PQ、、PV 节点功率的不平衡量,及PV 节点电压的不平衡量:取:000.0000.1)0(1)0(1)0(1j jf e U +=+=000.0000.1)0(2)0(2)0(2j jf e U +=+=节点3是平衡节点,保持000.0000.1333j jf e U +=+=为定值。
课程设计报告学生姓名:学号:学院:班级:题目: 电力系统潮流计算课程设计课设题目及要求一 .题目原始资料1、系统图:两个发电厂分别通过变压器和输电线路与四个变电所相连。
2、发电厂资料: 母线1和2为发电厂高压母线,发电厂一总装机容量为( 300MW ),母线3为机压母线,机压母线上装机容量为( 100MW ),最大负荷和最小负荷分别为50MW 和20MW ;发电厂二总装机容量为( 200MW )。
3、变电所资料:(一) 变电所1、2、3、4低压母线的电压等级分别为:35KV 10KV 35KV10KV(二) 变电所的负荷分别为:60MW 40MW 40MW 50MW(三) 每个变电所的功率因数均为cos φ=0.85;变电所1变电所母线 电厂一 电厂二(四) 变电所1和变电所3分别配有两台容量为75MVA 的变压器,短路损耗414KW ,短路电压(%)=16.7;变电所2和变电所4分别配有两台容量为63MVA 的变压器,短路损耗为245KW ,短路电压(%)=10.5;4、输电线路资料:发电厂和变电所之间的输电线路的电压等级及长度标于图中,单位长度的电阻为Ω17.0,单位长度的电抗为Ω0.402,单位长度的电纳为S -610*2.78。
二、 课程设计基本内容:1. 对给定的网络查找潮流计算所需的各元件等值参数,画出等值电路图。
2. 输入各支路数据,各节点数据利用给定的程序进行在变电所在某一负荷情况下的潮流计算,并对计算结果进行分析。
3. 跟随变电所负荷按一定比例发生变化,进行潮流计算分析。
1) 4个变电所的负荷同时以2%的比例增大;2) 4个变电所的负荷同时以2%的比例下降3) 1和4号变电所的负荷同时以2%的比例下降,而2和3号变电所的负荷同时以2%的比例上升;4. 在不同的负荷情况下,分析潮流计算的结果,如果各母线电压不满足要求,进行电压的调整。
(变电所低压母线电压10KV 要求调整范围在9.5-10.5之间;电压35KV 要求调整范围在35-36之间)5. 轮流断开支路双回线中的一条,分析潮流的分布。
电力系统潮流上机课程设计教学指导书(适用于VC)一、基础知识的要求1.掌握VC环境的建立,建立一个空的或者“hello world”的DOS工作平台。
2.掌握函数的定义,头文件的使用,函数调用方法。
了解类的概念,建立类对象,学会调用类成员函数。
3.掌握数据文件的打开和关闭,学会自己根据系统数据建立数据文件,并读取数据。
学会将计算结果写到文件中。
4.对VC环境的掌握:i.打开工作平台:file->open workspace(file->recent workspace)ii.打开工作平台中的文件:选择fileview页,在source file中双击想打开的.cpp文件,或者在header file中双击想打开的.h文件。
iii.保存文件:ctrl+siv.将已有文件加入当前工作平台中:如果是cpp文件,则选中source file,单击右键,选择add file to folder,加入希望加入的cpp文件。
如果是.h文件,则选中header file,单击右键,选择add file to folder,加入希望加入的头文件。
v.编译:所有的程序写完之后需要进行编译。
编译可以选择Build按钮,也可以用快捷键F7进行编译。
如果有语法错误,则会在output的build页显示出来,逐个双击定位,处理错误,解决问题。
vi.调试:设置断点(F9),进入调试状态或者执行至下一个断点(F5,或者点击界面上的“Go”按钮),单步执行(F10),进入子函数(F11),执行到鼠标所在位置(ctrl+F10),结束调试(shift+F5)。
vii.运行:运行可以点击界面上的“叹号”按钮,也可以使用快捷键ctrl+F5,使用这个命令时,所设置的断点是不起作用的。
5.《电力系统稳态分析》教材中第四章涉及到的基本知识。
二、设计注意事项1.变量的定义要写在变量被使用之前,变量不能被重复定义。
2.VC是区分大小写的。
电力系统潮流上机设计教学指导书一、综合实验的概述1.目的:培养学生对电力销售环节的分析能力,把握电力零售的相关知识。
2.使用的工具:零售市场仿真实验软件3.差不多要求:把握电力销售环节的交易模型;销售电价的分类;电费运算;分析不同负荷特性的电力用户执行不同电价类别的差异。
4.考核方法:总成绩=出勤情形+设计报告+面试成绩5.参考教材:区域电力市场电价机制张粒子、郑华北京:中国电力出版社, 2004二、综合实验的要求1.依照各类用户负荷数据,把握用户负荷特性分析方法。
2.把握销售电价中的电价结构。
3.把握各类电价制度的电费结算模型。
4.把握各类电价制度的特点、适用的用户范畴。
5.对EXCEL的把握:打开文件、储存文件、另存文件;读取运算数据;制作分析表格;依照分析数据制作展现图。
三、综合实验的内容实验一:电力用户负荷特性分析(一)基础理论知识预备:电力用户负荷特性1)电力用户用电负荷曲线:将电力用户的有功负荷,按时刻序列绘制成的图形,称为负荷曲线。
(1) 日负荷曲线:标示出一天内每小时(每半小时、每15分钟)的负荷值,反映一天内负荷动态。
(2) 年负荷曲线:标示出一年内每月的最高负荷值,反映一年内各月负荷动态。
2)要紧的负荷特性指标及运算(1) 最高负荷:报告期(日、月、季、年)内记录的负荷中,数值最大的一个。
(2) 最低负荷:报告期(日、月、季、年)内记录的负荷中,数值最小的一个。
(3) 平均负荷:报告期内瞬时负荷的平均值,即负荷时刻数列时序平均数。
(4) 负荷率:平均负荷与最高负荷的比率。
负荷率(%)=⨯)k k W W 报告期最高负荷()报告期平均负荷(100% (5) 最小负荷率:报告期最低负荷与最高负荷的比率。
最小负荷率(%)=⨯)报告期最高负荷()报告期最低负荷(W W k k 100%(6) 峰谷差:最高负荷与最低负荷之差。
(7) 峰谷差率:峰谷差与最高负荷的比率。
峰谷差率(%)=⨯)k ()k (W W 最高负荷峰谷差100% (8) 同时率:综合最高负荷与各组成用户绝对最高负荷之和的比率,说明用电负荷同时性程度。
关于潮流计算的课程设计一、课程目标知识目标:1. 理解潮流计算的基本概念,掌握电力系统潮流的基本方程和算法;2. 学习并掌握使用相关软件工具进行潮流计算的方法;3. 了解潮流计算在电力系统运行中的应用及其重要性。
技能目标:1. 能够运用所学的潮流计算方法,解决实际的电力系统平衡问题;2. 培养学生运用计算机软件进行电力系统分析的能力;3. 通过团队协作,提高学生的问题分析、数据解读和报告撰写技巧。
情感态度价值观目标:1. 培养学生对电力系统的兴趣,激发他们探索电力科学领域的热情;2. 强化学生的安全意识,认识到电力系统运行安全的重要性;3. 增强学生的环保意识,理解潮流计算在节能减排中的作用;4. 通过团队合作,培养学生的沟通能力和集体荣誉感。
课程性质分析:本课程旨在帮助学生建立电力系统潮流计算的基本理论框架,并通过实践操作,使学生掌握相关技能,为后续专业课程打下坚实基础。
学生特点分析:考虑到学生所在年级的知识深度,课程设计将结合学生的理论基础和实际应用能力,以引导式教学为主,鼓励学生主动探索和思考。
教学要求:1. 教学内容与课本紧密相关,注重理论与实践相结合;2. 采用案例教学,提高课程的实用性和针对性;3. 教学过程中注重培养学生的创新能力和实际操作能力。
二、教学内容1. 潮流计算基本概念:介绍潮流计算的定义、作用以及在电力系统运行中的重要性;教材章节:第一章第三节。
2. 潮流计算基本方程:讲解电力系统的基本潮流方程,包括功率方程、电压方程和相角方程;教材章节:第二章第一节。
3. 潮流计算算法:介绍常用的潮流计算算法,如牛顿-拉夫逊法、快速分解法和PQ分解法;教材章节:第二章第二节。
4. 潮流计算软件应用:指导学生使用相关软件(如PSS/E、DIgSILENT PowerFactory等)进行潮流计算;教材章节:第三章。
5. 潮流计算案例分析:分析实际电力系统潮流计算案例,使学生了解潮流计算在实际工程中的应用;教材章节:第四章。
院系:电气与电子工程学院班级:学号:学生姓名:指导教师:刘宝柱设计周数:成绩:日期: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)在进行北京城市电网的潮流计算之前需要了解北京城市电网中所有的节点支路的相关数据,并对节点和支路分类。
《电力系统潮流上机》课程设计报告院系:电气与电子工程学院班级:电气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-=ε。
电力系统潮流上机计算报告
系别:
班级:
姓名:
学号:
2012年×月×日
一、程序说明
包括:程序设计思想、程序流程图、程序使用说明。
二、给定题目的手算过程(迭代两次)
包括:原题目、节点导纳矩阵、雅克比矩阵、第一次和第二次迭代
结果。
三、给定题目的程序计算结果
包括:原题目、节点导纳矩阵、雅克比矩阵、程序输入和输出文件(误差0.0001)。
四、编程特色与创新
包括:程序能够完成的基本功能;程序能够完成的高级功能(如:是否包括非标准变比变压器支路,是否采用了稀疏矩阵技术,
是否增加了人机对话界面,程序的通用性和实用性如何)。
五、总结
包括:手算结果与程序计算结果的分析比较;本次上机体会,如:独立编程体会、跟踪调试技能的掌握情况、C语言中结构体、
指针、文件输入输出的掌握情况等。
报告要求:
1.报告中除上面的第“三”项外,其他部分必须手写(最好使用黑色水笔)。
2.报告统一采用A4打印纸书写(留出页边距: 1.5~2厘米)。
不使用实验报告纸。
3.封面按上述格式书写。
4.装订统一在左侧1厘米,二个钉。
5.上述五部分内容必须齐全,各部分内容可以扩充。
6.报告书写要求字迹清楚,不得潦草,流程图中的框要用尺子画。
7.报告必须与本人提交程序吻合,否则取消成绩。
8.报告不得有雷同,否则全部取消成绩。
《电力系统潮流上机》课程设计报告院系:电气与电子工程学院班级:学号:学生姓名:指导教师:刘宝柱设计周数:成绩:日期: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)在进行北京城市电网的潮流计算之前需要了解北京城市电网中所有的节点支路的相关数据,并对节点和支路分类。
处理PQ节点时要了解节点的注入有功和无功功率;PV节点要了解节点电压大小,注入有功功率及节点所能提供的最大和最小无功功率;对于平衡节点要了解节点的电压大小、相位及节点所能提供的最大和最小有功无功功率,此外还需的节点数据有:支路的电阻电抗电纳支路变压器的变比及线路所能传输的最大容量等。
(3) 计算结果分析:应考虑PQ节点的电压是否过高或过低,分析PV节点的电压幅值是否正常及无功功率是否超出范围,分析平衡节点有功无功功率是否在节点所能提供的范围之内,分析各支路的功率,看是否超出线路传输的最大容量,分析整个系统的网损是否达到标准3.3设计中遇到的问题和解决的办法。
开始时我们进行原始数据的读入和读出,这一部分的程序比较简单,所以我们都很顺利都写出程序。
接下来就是节点导纳矩阵的形成。
编写这一部分程序时我刚开始没有考虑到变压器支路,所以当原始数据改为四节点四支路的数据时,节点导纳矩阵的结果就出错了。
经过改正之后终于得出了正确的节点导纳矩阵。
然后就是计算不平衡功率和雅可比矩阵。
计算不平衡功率时候,为了方便计算,我将各节点的不平衡功率放到一个一维数组里面。
雅可比矩阵的形成这一部分上我用了很长时间才将程序写出来。
首先是对角线上的H,N,J,L的计算公式写的很复杂,与同学讨论之后我发现可以先将每个节点的注入电流放到一个数组里面,这样后面用来求H ,N ,J,L 时就比较方便。
在这一模块的过程中,我没有注意到雅克比矩阵中PQ节点和PV节点的对应关系,而且我没有将节点号取出来,以致出现了错误的编程思路和错误的结果,后来在同组同学的帮助下我又进行了第二次的编程,思路是正确的,但是运行出来的结果却是错误的,我想应该是我哪个公式出了错,最终经过反复检查后终于得到了正确结果,在编写雅克比矩阵还需要注意H、J、N、L四个矩阵的形成以及循环的顺序,然后就是如何将雅可比矩阵的各个元素放到一个二维数组里,这一部分看起来很简单,但是却很容易出错,经过和同组同学的认真思考和检查,终于得出来了正确的雅可比矩阵。
三、课程设计总结或结论两周的潮流上机计算已经结束了,虽然时间很短,可是我却学到了很多。
潮流上机计算除了需要一定的电力系统稳态知识以外,还要有良好的C语言基础。
潮流上机的过程让我们巩固了潮流计算过程和方法,以及稳态的基础知识。
此外,从潮流数据的读入和写出,节点导纳矩阵的形成,功率不平衡量的计算,到形成雅克比矩阵和功率损耗的计算,我不断地进行编程和学习,在两周的潮流上机过程中,我完成了程序的编写,掌握了计算机潮流计算的基本方法与步骤,了解了计算机潮流计算的优点与重要性。
潮流计算最关键的是进行C语言的编程,没有良好的C语言基础很难快速地编出潮流计算程序,编写时需要注意循环和嵌套的使用,我在编写的过程中由于循环没有用好以及粗心等原因,导致多次编出来的程序出现错误,这在潮流计算中都是不应该出现的。
四、参考文献1.《电力系统稳态分析》,陈珩,中国电力出版社,2007年,第三版;2. 《C 语言程序设计》,郑玲,中国电力出版社,2009年,第一版; 附录(设计流程图、程序、表格、数据等)是 否令迭代次数k=0对PQ 节点计算)()(.k i k iQ P ∆∆对PV 节点计算)(k i P ∆令节点号i=1计算雅可比矩阵各元素H ij (k)、N ij (k)、J ij (k)、L ij (k)增加节点号i=i+1解修正方程,由)()(.k i k iQ P ∆∆及雅可比矩阵用牛顿-拉夫逊法求各节点的Δe i ,Δf i求出max (|d e i|,|d f i|)判断迭代是否收敛|dei|,|dfi|<ε雅可比矩阵是否已形成,i>n? 计算节点的新电压e i (k+1)=e i (k)+δe i (k)f i (k+1)=f i (k)+δf i (k)增加迭代次数k=k+1计算平衡节点的功率及线路功率形成节点导纳矩阵设非平衡节点电压初值e i (0) ,f i (0)输入原始数据 启动1、设计流程图停止2、程序数据输入模块struct Line{int Num,NumI,NumJ;float R,X,B,K;};struct Bus{int Num ;float Volt,Phase,GenP,GenQ,LoadP,LoadQ; int Type;};struct Shunt{int Num,NumI;float G,B;};#include"stdio.h"#include"string.h"#include"math.h"#include"stdlib.h"#define NBUS 4#define NLINE 4/* Global variables */int nL,nSH,nB,nVA;float X[NBUS];int L;void main(){FILE *fp;int i,j,k,l,h,c;int i1,i2,i3;float d1,d2,d3,d4,d5,d6,r,x,g,b,e;struct Line sL[NLINE];struct Bus sB[NBUS];struct Shunt sSH[NBUS];float YG[NBUS][NBUS],YB[NBUS][NBUS];i1=i2=i3=0;d1=d2=d3=d4=d5=d6=0.0;e=0.00001;for(i=0;i<NBUS;i++)/*Read the input data*/if((fp=fopen("data.txt","r"))==NULL){ printf("Can not open the file named 'data.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].Type=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; }for(i=0;i<nSH;i++){sSH[i].Num=sSH[i].NumI=0;sSH[i].G=sSH[i].B=0.0;fscanf(fp,"%2d %3d %f",&i1,&i2,&d1);sSH[i].Num=i1;sSH[i].NumI=i2;sSH[i].B=d1;}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;if(fabs(sL[l].K-1.0)<0.000001){/*Normal lines or transformers*/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{/*abnormal transformer ratio*/if(fabs(sL[l].B)>0.000001){YG[i][i]=YG[i][i]+g/sL[l].B/sL[l].B;YG[j][j]=YG[j][j]+g;YB[i][i]=YB[i][i]+b/sL[l].B/sL[l].B;YB[j][j]=YB[j][j]+b;YG[i][j]=YG[i][j]-g/sL[l].B;YG[j][i]=YG[j][i]-g/sL[l].B;YB[i][j]=YB[i][j]-b/sL[l].B;YB[j][i]=YB[j][i]-b/sL[l].B;}}}/* 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);/* 设定电压初值*/float V[NBUS][NBUS]={0};for(i=1;i<NBUS+1;i++)if(sB[i-1].Type==0){V[i][0]=1.0;V[i][1]=0.0;}for(i=1;i<NBUS+1;i++)if(sB[i-1].Type==1){V[i][0]=sB[i-1].Volt;V[i][1]=0.0;}for(i=1;i<NBUS+1;i++)if(sB[i-1].Type==2){V[i][0]=(sB[i-1].Volt)*cos(sB[i-1].Phase);//es V[i][1]=(sB[i-1].Volt)*sin(sB[i-1].Phase);//fs}//输出电压初值if((fp=fopen("电压初值.txt","w"))==NULL){printf("Can not open the file named '电压初值.txt' \n"); exit(0);}fprintf(fp,"---电压初值---\n");for(i=1;i<NBUS+1;i++)for(j=1;j<2;j++)fprintf(fp,"Y(%2d )=(%10.5f,%10.5f)\n",i,V[i][0],V[i][1]); if(fp!=NULL) fclose(fp);for(c=1;;c++){/* 计算偏移量*/float PQ[100][2]={0},dV[100]={0},df[100]={0};for(i=1;i<NBUS;i++)if(sB[i-1].Type==0)//计算PQ节点{for(j=1;j<NBUS+1;j++){PQ[i][0]-=V[i][0]*(YG[i][j]*V[j][0]-YB[i][j]*V[j][1])+V[i][1]*(YG[i][j]*V[j][1]+YB[i][j ]*V[j][0]);PQ[i][1]-=V[i][1]*(YG[i][j]*V[j][0]-YB[i][j]*V[j][1])-V[i][0]*(YG[i][j]*V[j][1]+YB[i][j ]*V[j][0]);}}if(sB[i-1].Type==1)//计算PV节点{for(j=1;j<NBUS;j++){PQ[i][0]-=V[i][0]*(YG[i][j]*V[j][0]-YB[i][j]*V[j][1])+V[i][1]*(YG[i][j]*V[j][1]+YB[i][j ]*V[j][0]);}dV[i]=V[i][1]*V[i][1]+V[i][0]*V[i][0];}}h=1;for(i=1;i<NBUS;i++){if(sB[i-1].Type==0){df[h]=PQ[i][0]+sB[i-1].GenP-sB[i-1].LoadP;h++;df[h]=PQ[i][1]+sB[i-1].GenQ-sB[i-1].LoadQ;h++;}if(sB[i-1].Type==1){df[h]=PQ[i][0]+sB[i-1].GenP-sB[i-1].LoadP;h++;df[h]=sB[i-1].Volt*sB[i-1].Volt-dV[i];h++;}}/*输出偏移量*/if((fp=fopen("第一次迭代偏移量.txt","w"))==NULL){printf("Can not open the file named '第一次迭代偏移量.txt' \n");exit(0);}fprintf(fp,"--- 第一次迭代偏移量---\n");for(i=1;i<h;i++)fprintf(fp,"df(%2d )=(%10.5f)\n",i,df[i]);if(fp!=NULL) fclose(fp);}/* 计算雅克比矩阵*/float ga[2*(NBUS-1)]={0},gb[2*(NBUS-1)]={0};floatH[2*(NBUS-1)][2*(NBUS-1)]={0},J[2*(NBUS-1)][2*(NBUS-1)]={0},L[2*(NBUS-1)][2*(NBUS-1)]={0},N [2*(NBUS-1)][2*(NBUS-1)]={0},JJ[2*(NBUS-1)][2*(NBUS-1)]={0},R[2*(NBUS-1)][2*(NBUS-1)]={0},S [2*(NBUS-1)][2*(NBUS-1)]={0};for(i=1;i<NBUS;i++){for(j=1;j<NBUS;j++){if(i==j){ga[i]=-(PQ[i][0]*V[i][0]+PQ[i][1]*V[i][1])/(V[i][0]*V[i][0]+V[i][1]*V[i][1]); //计算a gb[i]=-(PQ[i][0]*V[i][1]-PQ[i][1]*V[i][0])/(V[i][0]*V[i][0]+V[i][1]*V[i][1]); //计算bH[i][i]=-V[i][0]*YB[i][i]+YG[i][i]*V[i][1]+gb[i];N[i][i]=V[i][0]*YG[i][i]+YB[i][i]*V[i][1]+ga[i];J[i][i]=-V[i][0]*YG[i][i]-YB[i][i]*V[i][1]+ga[i];L[i][i]=V[i][1]*YG[i][i]-YB[i][i]*V[i][0]-gb[i];R[i][i]=2*V[i][1];S[i][i]=2*V[i][0];}else{H[i][j]=-V[i][0]*YB[i][j]+YG[i][j]*V[i][1]; //非对角线元素正确。