九节点系统潮流计算编程牛N-R法

合集下载

三机九节点潮流暂态MATLAB仿真

三机九节点潮流暂态MATLAB仿真

三机九节点潮流暂态MATLAB仿真院系: 自动化学院专业:电力系统及其自动化学号: 姓名: 时间: 1 研究对象1.1 三机九节点系统模型100MW35MVar7239j0.0585j0.06250.0119 + j0.10080.0085 + j0.072B/2 = j0.1045B/2 = j0.0745230/13.818/23018kV13.8kV8230kV230kVB/2 = j0.179B/2 = j0.088B/2 = j0.1530.010 + j0.0850.032 + j0.161 56125MW90MW50MVar30MVarB/2 = j0.079230kV40.017 + j0.0790.039 + j0.170j0.057616.5/23016.5kV1图1.1 WSCC-9系统模型图1.1是一个三机九节点的系统阻抗图,图中给出的阻抗参数都是以100MVA为基准的标幺值。

该图中包括三台发电机,三台双绕组变压器,九条母线(节点)和三个负荷。

本文将对该系统的动态过程进行相应的仿真分析。

1.2 系统参数1.2.1 节点参数按照节点类型,9个节点分为,给出已知参数如下表:表1.1 节点已知参数节点类型电压幅值电压角度发电机有功发电机无功负荷有功负荷无功1 Vθ 1.040 0 0.7160 0.2705 0 02 PV 1.025 1.6300 0.0665 0 03 PV 1.025 0.8500 -0.1086 0 04 PQ 0 05 PQ 1.2500 0.50006 PQ 0.9000 0.30007 PQ 0 08 PQ 1(0000 0.35009 PQ 0 0上表中发电机有功、无功出力和负荷的有功无功功率均为以100MVA为基准时的标幺值。

1.2.2 支路参数表1.2 支路参数首节点末节点电阻电抗电纳一半4 5 0.0100 0.0850 0.08804 6 0.0170 0.0920 0.07905 7 0.0320 0.1610 0.15306 9 0.0390 0.1700 0.17907 8 0.0085 0.0720 0.07458 9 0.0119 0.1008 0.10451 4 0.0000 0.0576 0.00002 7 0.0000 0.0625 0.00003 9 0.0000 0.0586 0.0000上表中所有的参数均为标幺值,对于变压器支路。

九节点电力系统分析潮流程序设计

九节点电力系统分析潮流程序设计

九节点电力系统分析潮流程序设计下载提示:该文档是本店铺精心编制而成的,希望大家下载后,能够帮助大家解决实际问题。

文档下载后可定制修改,请根据实际需要进行调整和使用,谢谢!本店铺为大家提供各种类型的实用资料,如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等,想了解不同资料格式和写法,敬请关注!Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you! In addition, this shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!1.引言随着电力系统规模的不断扩大和电力市场的逐渐开放,电力系统潮流分析成为了电力系统规划、运行和市场交易中不可或缺的一部分。

牛顿拉夫逊法计算潮流步骤

牛顿拉夫逊法计算潮流步骤

牛顿拉夫逊法计算潮流步骤牛顿拉夫逊法(Newton-Raphson method)是一种用于求解非线性方程组的迭代方法,它可以用来计算电力系统潮流的解。

潮流计算是电力系统规划和运行中的重要任务,它的目标是求解电力系统中各节点的电压幅值和相角,以及线路的功率流向等参数,用于分析电力系统的稳定性和安全性,以及进行电力系统规划和调度。

下面是使用牛顿拉夫逊法计算潮流的一般步骤:步骤1:初始化首先,需要对电力系统的各个节点(包括发电机节点和负荷节点)的电压幅值和相角进行初始化,一般可以使用其中一种估计值或者历史数据作为初始值。

步骤2:建立潮流方程根据电力系统的潮流计算模型,可以建立节点电压幅值和相角的平衡方程,一般采用节点注入功率和节点电压的关系来表示。

潮流方程一般是一个非线性方程组,包含了各个节点之间的复杂关系。

步骤3:线性化方程组将潮流方程组进行线性化处理,一般采用泰勒展开的方法,将非线性方程组变为线性方程组。

线性化的过程需要计算雅可比矩阵,即方程组中的系数矩阵。

步骤4:求解线性方程组利用线性方程组的求解方法,比如高斯消元法或LU分解法等,求解线性方程组,得到电压幅值和相角的修正量。

步骤5:更新节点电压根据线性方程组的解,更新各个节点的电压幅值和相角,得到新的节点电压。

步骤6:检查收敛性判断节点电压的修正量是否小于设定的收敛阈值,如果满足收敛条件,则停止迭代;否则,返回步骤3,循环进行线性化方程组和线性方程组的求解。

步骤7:输出结果当潮流计算收敛时,输出最终的节点电压幅值和相角,以及线路的功率流向等参数。

牛顿拉夫逊法是一种高效、快速且收敛性良好的方法,广泛应用于电力系统潮流计算。

在实际应用中,可能会遇到多次迭代或者收敛性不好的情况,此时可以采用退火技术或其他优化算法进行改进。

此外,牛顿拉夫逊法的计算也可以并行化,利用多核处理器或者分布式计算集群来加速计算过程。

总之,牛顿拉夫逊法是一种重要的潮流计算方法,通过迭代计算逼近非线性方程组的解,可以得到电力系统中各节点的电压幅值和相角,用于分析电力系统的稳定性和安全性。

牛顿拉夫逊潮流计算[整理版]

牛顿拉夫逊潮流计算[整理版]

牛顿拉夫逊潮流计算 [整理版 ]float G[N][N],B[N][N]; // struct // 阻抗参数int nl; // 左节点int nr; // 右节点牛顿拉夫逊潮流计算// 整个程序为 :// 牛拉法解潮流程序#include<stdio.h>#include<math.h>#define N 4 // 节点数#define n_PQ 2 //PQ 节点数#define n_PV 1 //PV 节点数#define n_br 5 // 串联支路数void main() void disp_matrix(float *disp_p,int disp_m,int disp_n); // 矩阵显示函 float Us[2*N]={1.0,0,1.0,0,1.05,0,1.05,0}; // 电压初值 float Ps[N]={0,-0.5,0.2}; // 有功初值float Qs[N]={0,-0.3}; // 无功初值各几点电导电纳}ydata[n_br]={ {1,2,0,0.1880,-0.6815,0.6040},{1,3,0.1302,0.2479,0.0129,0.0129},{1,4,0.1736,0.3306,0.0172,0.0172},{3,4,0.2603,0.4959,0.0259,0.0259},{2,2,0,0.05,0,0} };float Z2; //ZA2=RA2+XA2 各串 联阻抗值的平方存储电压修正值float mid1[N],mid2[N],dS[2*(N-1)]; //mid1对角线元素的中间值 ,dS 存储 PQUi 勺不平衡量float Jacob[2*(N-1)][2*(N-1)],inv_J[2*(N-1)][2*(N-1)];// 雅克比行列式float dPQU=1.0; //PQU 不平衡量最大值int kk=0; // 迭代次数int i,j,k;float t;float R; // 串联电阻值 float X; // 串联电抗值 float Bl; // 左节点并联电导 float Br; // 右节点并联电纳float e[N],f[N],dfe[2*(N-1)]; //e,f存储电压的 x 轴分量和 y 轴分量 ,dfe、mid2 存储计算雅克比行列式// 形成导纳矩阵 for(i=0;i<N;i++) for(j=0;j<N;j++)G[i][j]=0;B[i][j]=0;for(i=0;i<n_br;i++) if(ydata[i].nl!=ydata[i].nr)Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X); // 串联阻抗等效导纳值// 非对角元素G[ydata[i].nl-1][ydata[i].nr-1]=(-ydata[i].R)/Z2;B[ydata[i].nl-1][ydata[i].nr-1]=ydata[i].X/Z2;G[ydata[i].nr-1][ydata[i].nl-1]=(-ydata[i].R)/Z2;B[ydata[i].nr-1][ydata[i].nl-1]=ydata[i].X/Z2;// 对角元素 float Pij[n_br]; //存储线路 i->j 的有功 float Qij[n_br]; //存储线路 i->j 的无功 float Pji[n_br]; //存储线路 j->i 的有功 float Qji[n_br]; //存储线路 j->i 的无功 float dPij[n_br]; //存储线路 i->j 的有功损耗 float dQij[n_br]; //存储线路 i->j 的无功损耗 存储线路潮流计算时的中间float AA,BB,CC,DD; //G[ydata[i].nl-1][ydata[i].nl-1]+=ydata[i].R/Z2; G[ydata[i].nr-1][ydata[i].nr-1]+=ydata[i].R/Z2; B[ydata[i].nl-1][ydata[i].nl-1]+=(-ydata[i].X/Z2); B[ydata[i].nr-1][ydata[i].nr-1]+=(-ydata[i].X/Z2); // 并联导纳等效导纳值B[ydata[i].nl-1][ydata[i].nl-1]+=ydata[i].Bl;B[ydata[i].nr-1][ydata[i].nr-1]+=ydata[i].Br; elseG[ydata[i].nl-1][ydata[i].nr-1]+=ydata[i].R;B[ydata[i].nl-1][ydata[i].nr-1]+=ydata[i].X;printf("G=\n");disp_matrix(*G,N,N);printf("B=\n");disp_matrix(*B,N,N);// 分离e,ffor(i=0;i<N;i++) e[i]=Us[2*i];f[i]=Us[2*i+1];// 主程序while(dPQU>0.00001)// 计算功率不平衡量for(i=0;i<N-1;i++) mid1[i]=0;mid2[i]=0;for(j=0;j<N;j++) mid1[i]=mid1[i]+G[i][j]*e[j]-B[i][j]*f[j]; mid2[i]=mid2[i]+G[i][j]*f[j]+B[i][j]*e[j];dS[2*i]=Ps[i]-(e[i]*mid1[i]+f[i]*mid2[i]);if(i<n_PQ) dS[2*i+1]=Qs[i]-(f[i]*mid1[i]-e[i]*mid2[i]); else dS[2*i+1]=Us[2*i]*Us[2*i]-(e[i]*e[i]+f[i]*f[i]); dPQU=0;for(i=0;i<2*(N-1);i++) if(dS[i]<0&&dPQU<-dS[i]) dPQU=-dS[i];else if(dS[i]>0&&dPQU<dS[i]) dPQU=dS[i];if(dPQU>0.00001) kk++;// 形成雅克比行列式for(i=0;i<2*(N-1);i++) for(j=0;j<2*(N-1);j++)Jacob[i][j]=0;for(j=0;j<N-1;j++)// 求H,N for(i=0;i<N-1;i++) if(i!=j)Jacob[2*i][2*j]=B[i][j]*e[i]-G[i][j]*f[i];Jacob[2*i][2*j+1]=-G[i][j]*e[i]-B[i][j]*f[i];elseJacob[2*i][2*i]=B[i][i]*e[i]-G[i][i]*f[i]-mid2[i]; Jacob[2*i][2*i+1]=-G[i][j]*e[i]-B[i][j]*f[i]-mid1[i];// 求J,L for(i=0;i<n_PQ;i++) if(i!=j)Jacob[2*i+1][2*j]=G[i][j]*e[i]+B[i][j]*f[i];Jacob[2*i+1][2*j+1]=B[i][j]*e[i]-G[i][j]*f[i];elseJacob[2*i+1][2*i]=G[i][j]*e[i]+B[i][j]*f[i]-mid1[i]; Jacob[2*i+1][2*i+1]=B[i][j]*e[i]-G[i][j]*f[i]+mid2[i];}// 求R,S for(i=n_PQ;i<N-1;i++) if(i==j)Jacob[2*i+1][2*i]=-2*f[i];Jacob[2*i+1][2*i+1]=-2*e[i];// 雅克比行列式求逆// 初始化inv_J[N][N] for(i=0;i<2*(N-1);i++)for(j=0;j<2*(N-1);j++) if(i!=j) inv_J[i][j]=0;else inv_J[i][j]=1;// 将原矩阵化简为对角阵for(i=0;i<2*(N-1);i++)for(j=0;j<2*(N-1);j++) if(i!=j) t=Jacob[j][i]/Jacob[i][i]; for(k=0;k<2*(N-1);k++)Jacob[j][k]-=Jacob[i][k]*t;inv_J[j][k]-=inv_J[i][k]*t;// 原矩阵各对角元素化为1,画出逆矩阵for(i=0;i<2*(N-1);i++) if(Jacob[i][i]!=1) t=Jacob[i][i]; for(j=0;j<2*(N-1);j++) inv_J[i][j]=inv_J[i][j]/t;// 求电压修正值for(i=0;i<2*(N-1);i++) dfe[i]=0;for(j=0;j<2*(N-1);j++) dfe[i]-=inv_J[i][j]*dS[j];for(i=0;i<N-1;i++) e[i]+=dfe[2*i+1];f[i]+=dfe[2*i];else break;// 循环结束// 求平衡节点功率mid1[N-1]=0;mid2[N-1]=0;目〒o 〒N j++) midpN 厶Hmid2运厶+G 运厶=rf 曰+B 运厶PS2 厶 HeLN 二「mid 二N 厶+nN 二 :Qs 运厶吕 N 二「mid 二 N 二〒 〔二 fo 「(nnlpQ 天 N 丄 i ++)Qs 日吕『mid 二〒e 耳mid2wpl1nff(=kkH&2rr 』kkx prinm.pH) f o 「(i H O =A N =++)Pl1mf(=%9.4r 「ps=2 prinmvIQH) fo 「(no 天 N =++) pl1mf(=%9.4r 「QS=2 prinmvIeH) fo 「(no 天 N =++) Pl1mf(=%9.4r 「e=2 prinmvIfH) f o 「(i H O =A N =++)printf("%9.4f",f[i]);printf("\n");// 求线路上的潮流// 计算S[i][j]for(i=0;i<n_br;i++)if(ydata[i].nl!=ydata[i].nr)Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X);AA=-f[ydata[i].nl-1]*ydata[i].Bl+(e[ydata[i].nl-1]-e[ydata[i].nr- 1])*ydata[i].R/Z2+(f[ydata[i].nl-1]-f[ydata[i].nr-1])*ydata[i].X/Z2;BB=-e[ydata[i].nl-1]*ydata[i].Bl-(f[ydata[i].nl-1]-f[ydata[i].nr- 1])*ydata[i].R/Z2+(e[ydata[i].nl-1]-e[ydata[i].nr-1])*ydata[i].X/Z2;Pij[i]=e[ydata[i].nl-1]*AA-f[ydata[i].nl-1]*BB;Qij[i]=e[ydata[i].nl-1]*BB+f[ydata[i].nl-1]*AA;printf("S[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nl,ydata[i].nr,Pij[i] ,Qij [i]);dPij[i]=Pij[i]+Pji[i];dQij[i]=Qij[i]+Qji[i];printf("\n");// 计算S[j][i]for(i=0;i<n_br;i++)if(ydata[i].nl!=ydata[i].nr)Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X);CC=-f[ydata[i].nr-1]*ydata[i].Br+(e[ydata[i].nr-1]-e[ydata[i].nl- 1])*ydata[i].R/Z2+(f[ydata[i].nr-1]-f[ydata[i].nl-1])*ydata[i].X/Z2;DD=-e[ydata[i].nr-1]*ydata[i].Br-(f[ydata[i].nr-1]-f[ydata[i].nl- 1])*ydata[i].R/Z2+(e[ydata[i].nr-1]-e[ydata[i].nl-1])*ydata[i].X/Z2;Pji[i]=e[ydata[i].nr-1]*CC-f[ydata[i].nr-1]*DD;Qji[i]=e[ydata[i].nr-1]*DD+f[ydata[i].nr-1]*CC;printf("S[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nr,ydata[i].nl,Pji[i] ,Qji [i]);printf("\n");// 计算dS[i][j]for(i=0;i<n_br;i++)if(ydata[i].nl!=ydata[i].nr)dPij[i]=Pij[i]+Pji[i];dQij[i]=Qij[i]+Qji[i];printf("dS[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nl,ydata[i].nr, dPij[i],dQij[i]);printf("\n");}// 主程序结束// 矩阵显示函数void disp_matrix(float *disp_p,int disp_m,int disp_n)int i,j;for(i=0;i<disp_m;i++)for(j=0;j<disp_n;j++)printf("%9.4f",*(disp_p+i*disp_m+j));printf("\n");printf("\n");}。

三机九节点电力系统仿真matlab【精品毕设、无需降重】

三机九节点电力系统仿真matlab【精品毕设、无需降重】

电力系统仿真作业------------三机九节点电力系统暂态仿真学院:能源与动力工程学院专业:电力系统及其自动化学号:姓名:***导师:授课教师:目录一、概述 (1)二、课程主要任务 (1)1.系统数据 (1)2.潮流计算 (2)3.负荷等效和支路简化 (4)4.求解电磁功率 (5)5.求解运动方程 (5)6.程序清单 (7)(1).主程序: (7)(2).极坐标转换成直角坐标函数pol2rect(V,del) (16)(3).直角坐标转换成极坐标函数rect2pol(Z) (16)(4).求解微分方程所用的得到微分量的函数Gen_fw(t,X,Y_Gen,E,Pm0,Tj) (16)三、课程总结及心得体会 (16)四、参考文献 (17)于永生电力系统仿真作业一、概述在动态稳定分析中,系统由线性化的微分方程组和代数方程组描写,并用经典的或现代的线性系统理论来进行稳定分析,分析可以在时域或频域进行。

当用计算机和现代线性系统理论分析时,常把系统线性化的微分方程组和代数方程组消去代数变量,化为状态方程形式,并广泛采用特征分析进行稳定分析。

电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。

其暂态过渡过程不仅包括电磁方面的过渡过程,而且还有机电方面的过渡过程。

由此可见,电力系统的数学模型是一个强非线性的高维状态方程组。

在动态稳定仿真中使用简单的电力系统模型,发电机用三阶模型表示。

二、课程主要任务本次课程主要应用P. M. Anderson and A. A. Fouad编写的《Power System Control and Stability》一书中所引用的Western System Coordinated Council (WSCC)三机九节点系统模型。

1.系统数据其中,节点数据如下:%节点数据% 节点电压电压发电机发电机负荷负荷节点% 号幅值相角有功无功有功无功类型(1PQ 2PV 3平衡)N=[ 1 1.04 0 0.7164 0.2705 0 0 32 1.025 0 1.63 0.0665 0 0 23 1.025 0 0.85 -0.1086 0 0 24 1 0 0 0 0 0 15 1 0 0 0 1.25 0.5 16 1 0 0 0 0.9 0.3 17 1 0 0 0 0 0 18 1 0 0 0 1 0.35 19 1 0 0 0 0 0 1];其中,支路数据如下:% 线路数据% 首端末端电阻电抗电纳(1/2) 变压器非标准变比L=[4 5 0.01 0.085 0.088 14 6 0.017 0.092 0.079 15 7 0.032 0.161 0.153 16 9 0.039 0.17 0.179 17 8 0.0085 0.072 0.0745 11电力系统仿真作业于永生28 9 0.0119 0.1008 0.1054 11 4 0 0.0576 0 12 7 0 0.0625 0 13 9 0 0.0586 0 1]; 发电机数据如下:% 发电机母线Xd Xd' Td0' Xq Xq' Tq0’Tj XfGe=[ 1 1 0.1460 0.0608 8.96 0.0969 0.0969 0 47.28 0.05762 2 0.8958 0.1198 6.00 0.8645 0.1969 0.535 12.80 0.06253 3 1.3125 0.1813 8.59 1.2578 0.2500 0.600 6.02 0.0585];系统电路结构拓扑图如下:图1 WSCC 3机9节点系统(所有参数以100MV A为基准值的标幺值)2.潮流计算首先进行潮流计算,采用牛顿拉夫逊迭代法,电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算,其任务是在已知某些运行参数的情况下,计算出系统中全部于永生 电力系统仿真作业 3的运行参数,一般来说,各个母线所供负荷的功率是已知的,各个节点电压是未知的(平衡节点除外),可以根据网络结构形成节点导纳矩阵,然后由节点导纳矩阵和网络拓扑结构列写功率方程,由于功率方程里功率是已知的,电压的幅值和相角是未知的,这样潮流计算的问题就转化为求解非线性方程组的问题了。

潮流计算程序及计算结果

潮流计算程序及计算结果

附表1:计算机计算潮流程序:%本程序的功能是用牛顿——拉夫逊法进行潮流计算% B1矩阵:1、支路首端号;2、末端号;3、支路阻抗;4、支路对地电纳% 5、支路的变比;6、支路首端处于K侧为1,1侧为0% B2矩阵:1、该节点发电机功率;2、该节点负荷功率;3、节点电压初始值% 4、PV节点电压V的给定值;5、节点所接的无功补偿设备的容量% 6、节点分类标号clear;n=13;%input('请输入节点数:n=');nl=13;%input('请输入支路数:nl=');isb=1;%input('请输入平衡母线节点号:isb=');pr=0.00001;%input('请输入误差精度:pr=');B1=[];%input('请输入由支路参数形成的矩阵:B1=');B2=[];%input('请输入各节点参数形成的矩阵:B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl); %-------修改部分------------ym=0;SB=100;UB=220;%ym=input('您输入的参数是标么值?(若不是则输入一个不为零的数值)'); if ym~=0%SB=input('请输入功率基准值:SB=');%UB=input('请输入电压基准值:UB=');YB=SB./UB./UB;BB1=B1;BB2=B2;for i=1:nlB1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;enddisp('B1矩阵B1=');disp(B1)for i=1:nB2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;enddisp('B2矩阵B2=');disp(B2)end% % %---------------------------------------------------for i=1:nl %支路数if B1(i,6)==0 %左节点处于低压侧p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元K侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧end%求导纳矩阵disp('导纳矩阵Y=');disp(Y)%----------------------------------------------------------G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部for i=1:n %给定各节点初始电压的实部和虚部e(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4); %PV节点电压给定模值endfor i=1:n %给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(i,i)=B(i,i)+B2(i,5); %i节点无功补偿量end%=========================================================== ========P=real(S);Q=imag(S);ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;while IT2~=0IT2=0;a=a+1;for i=1:nif i~=isb %非平衡节点C(i)=0;D(i)=0;for j1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)endP1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求P',Q'V2=e(i)^2+f(i)^2; %电压模平方%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素=========if B2(i,6)~=3 %非PV节点DP=P(i)-P1; %节点有功功率差DQ=Q(i)-Q1; %节点无功功率差%=============== 以上为除平衡节点外其它节点的功率计算=================%================= 求取Jacobi矩阵===================for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df=dQ/deX3=X2; % X2=dp/df X3=dQ/deX4=-X1; % X1=dP/de X4=dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;elseif j1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△PJ(m,q)=X2;endendelse%=============== 下面是针对PV节点来求取Jacobi矩阵的元素===========DP=P(i)-P1; % PV节点有功误差DV=V(i)^2-V2; % PV节点电压误差for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/dfX5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;elseif j1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;endendendendend%========= 以上为求雅可比矩阵的各个元素===================== for k=3:N0 % N0=2*n (从第三行开始,第一、二行是平衡节点)k1=k+1;N1=N; % N=N0+1 即N=2*n+1扩展列△P、△Qfor k2=k1:N1 % 扩展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); % 非对角元规格化endJ(k,k)=1; % 对角元规格化if k~=3 % 不是第三行%======================================================== ====k4=k-1;for k3=3:k4 % 用k3行从第三行开始到当前行前的k4行消去for k2=k1:N1 % k3行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endif k==N0break;end%==========================================for k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endelsefor k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endendend%====上面是用线性变换方式将Jacobi矩阵化成单位矩阵=====for k=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N); %修改节点电压实部k1=k+1;f(L)=f(L)-J(k1,N); %修改节点电压虚部end%------修改节点电压-----------for k=3:N0DET=abs(J(k,N));if DET>=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足要求的节点数加1endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解"w=-J*V"disp('迭代次数:');disp(ICT1);disp('没有达到精度要求的个数:');disp(ICT2);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;end%=============== 计算各输出量=========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列):');disp(E);EE=E*UB;disp(EE);disp('-----------------------------------------------------');disp('各节点的电压大小V为(节点号从小到大排列):');disp(V);VV=V*UB;disp(VV);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);enddisp('各节点的功率S为(节点号从小到大排列):');disp(S);disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');SS=S*SB;disp(SS);disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致):');for i=1:nlp=B1(i,1);q=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))*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);enddisp(Si(p,q));SSi(p,q)=Si(p,q)*SB;ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);%disp(SSi(p,q));disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一致):');for i=1:nlp=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))*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);enddisp(Sj(q,p));SSj(q,p)=Sj(q,p)*SB;ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);%disp(SSj(q,p));disp('-----------------------------------------------------');enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一致):');for i=1:nlp=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i));DDS(i)=DS(i)*SB;ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);%disp(DDS(i));disp('-----------------------------------------------------');endfigure(1);subplot(2,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(2,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;subplot(2,2,3);bar(real(S));ylabel('节点注入有功');grid on;subplot(2,2,4);bar(Siz);ylabel('支路首端无功');grid on;1.冬季最大运行方式潮流计算结果:计算机运行的B1,B2阵如下:B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.88+0.545*i 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.88+0.545*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.2+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表:节点号电压值相角值支路标号首端功率末端功率支路功率损耗1 242.0000 0 1-2 148.9391+3.327637i -143-27.45476i 5.93913-24.1271i2 231.0000 -3.0341 1-11 89.1956+37.193i -88.19803-61.4687i 0.997519-24.2757i3 227.4154 -4.4383 11-3 88.19803+61.4687i -88-54.5i 0.19803+6.9687i4 226.3304 -4.9004 1-12 77.8364+36.4186i -77.2404-55.7062i 0.596039-19.2876i5 229.8318 -3.9547 12-4 77.2404+55.7062i -77-47.72i 0.24036+7.9862i6 217.7226 -8.2720 1-5 63.5438+12.2204i -61.8773-32.0321i 1.66656-19.8116i7 234.9326 -3.1047 5-6 77.2597+56.3501i -77-47.72i 0.25974+8.6301i8 226.0991 -6.2429 5-7 15.3825-24.3181i 15.5354+0.274108i 0.152961-24.044i9 231.0000 0.5851 7-8 88.20085+61.55009i -88-54.5i 0.20085+7.0501i10 231.0000 3.4950 1-7 40.806-6.05737i -40.0647-22.0555i 0.741264-28.1128i11 236.1931 -1.3348 7-13 -63.6716-39.7687i 64.0965+17.5832i 0.424934-22.1855i12 237.9346 -0.8892 13-9 -64.0965-17.5832i 64.2+21.0187i 0.10348+3.4355i13 238.1868 -2.2028 1-10 79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果图形:2.冬季最小运行方式潮流计算结果:计算机运行的B1B2矩阵如下:B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.616+0.3817*i 1 0 0 20 0.539+0.3817*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.1+0.06197*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]电压调整措施:变电所1、4变压器变比:+2.5% 水电厂变压器变比:-2.5%5 234.2241 -2.0051 (12--4) 54.0234+42.2706i -53.9-38.17i 0.12342+4.1006i6 226.2159 -4.8577 (1--5) 34.1669+4.22856 -33.6427-28.276i 0.524165-24.0474i7 235.1128 -0.4737 (5--6) 54.0179+37.3169i -53.9-33.4i 0.11789+3.9169i8 223.9941 -2.4603 (5--7) -20.3752-9.04094i 20.5371-15.4558i 0.161947-24.4968i9 231.0000 3.4429 (1--7) 11.0013+1.45186i -10.8258-31.4513i 0.175442-29.9995i10 231.0000 3.9321 (7--8) 53.9768+36.0957i -53.9-33.4i 0.076797+2.6957i11 238.4187 -0.9561 (7--13) -63.6881+10.8115i 64.0874-32.7763i 0.39932-21.9648i12 239.1742 -0.6185 (13--9) -64.0874+32.7763i 64.2-29.0386i 0.11258+3.7377i13 234.9468 0.6942 (1--10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机运行结果的图形:3.夏季最大运行方式计算机计算结果:计算机运行B1B2阵如下:B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.776+0.481*i 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.776+0.481*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表:10 231.0000 3.4950 (7,8) 77.7585+53.6634i -77.6-48.1i0.1585+5.5634i11 237.0938 -1.1858 (7,13) -113.4984+9.476406i 114.6926-28.38885i 1.19422-18.9124i12 239.1742 -0.6185 (13,9) -114.6926+28.38885 115-18.18369i0.307384+10.2052i13 233.3912 3.2303 (1,10) -79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果如图:4.夏季最小运行方式:计算机运行B1B2阵如下:B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.543+0.336*i 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.543+0.336*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]10 231.0000 3.9321 (7,8) 54.3735+36.2509i -54.3-33.67i 0.073528+2.5809i11 239.0099 -0.8518 (7,13) -113.4428+24.39707i 114.6754-43.79301i 1.23255-19.3959i12 239.8726 -0.5689 (13,9) -114.6754+43.79301i 115-33.01613i 0.324605+10.7769i13 235.9565 4.4510 (1,10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机计算结果如图:5.夏季故障运行状态:调压及无功补偿措施如下:变电所3的变压器变比为-2.5%,无功补偿容量为20Mvar。

牛顿-拉夫逊法潮流计算

牛顿-拉夫逊法潮流计算

目录摘要11.设计意义与要求2 1.1设计意义21.2设计要求32.牛顿—拉夫逊算法3 2.1牛顿算法数学原理:32.2 直角坐标系下牛顿法潮流计算的原理43 详细设计过程10 3.1节点类型103.2待求量103.3导纳矩阵103.4潮流方程113.5修正方程124.程序设计15 4.1 节点导纳矩阵的形成154.2 计算各节点不平衡量164.3 雅克比矩阵计算- 19 -4.4 LU分解法求修正方程- 22 -4.5 计算网络中功率分布- 25 -5.结果分析- 25 -6.小结- 29 -参考文献- 30 -附录:- 31 -摘要潮流计算是电力网络设计及运行中最基本的计算,对电力网络的各种设计方案及各种运行方式进行潮流计算,可以得到各种电网各节点的电压,并求得网络的潮流及网络中各元件的电力损耗,进而求得电能损耗。

在数学上是多元非线性方程组的求解问题,求解的方法有很多种。

牛顿—拉夫逊法是数学上解非线性方程式的有效方法,有较好的收敛性。

将牛顿法用于潮流计算是以导纳矩阵为基础的,由于利用了导纳矩阵的对称性、稀疏性及节点编号顺序优化等技巧,使牛顿法在收敛性、占用存、计算速度等方面都达到了一定的要求。

本文以一个具体例子分析潮流计算的具体方法,并运用牛顿—拉夫逊算法求解线性方程关键词:电力系统潮流计算牛顿—拉夫逊算法1.设计意义与要求1.1设计意义潮流计算是电力系统分析中的一种最基本的计算,他的任务是对给定运行条件确定系统运行状态,如各母线上的电压(幅值及相角)、网络中的功率分布及功率损耗等。

潮流计算的结果是电力系统稳定计算和故障分析的基础。

具体表现在以下方面:(1)在电网规划阶段,通过潮流计算,合理规划电源容量及接入点,合理规划网架,选择无功补偿方案,满足规划水平的大、小方式下潮流交换控制、调峰、调相、调压的要求。

(2)在编制年运行方式时,在预计负荷增长及新设备投运基础上,选择典型方式进行潮流计算,发现电网中薄弱环节,供调度员日常调度控制参考,并对规划、基建部门提出改进网架结构,加快基建进度的建议。

MATPOWER潮流计算使用说明

MATPOWER潮流计算使用说明

MATPOWER潮流计算使用说明一、MATPOWER安装和加载数据2.打开MATLAB软件,并在命令行上输入以下命令来加载MATPOWER软件包:```matlabaddpath(genpath('<MATPOWER安装路径>'))```注意,需要将`"<MATPOWER安装路径>"`替换为你的MATPOWER软件的安装路径。

3.加载示例数据集。

MATPOWER提供了一些示例数据集,可以直接使用这些数据集进行潮流计算。

```matlabcase9```这将加载一个名为`case9`的数据集,它包含一个9节点的电力系统。

二、设置潮流计算参数在进行潮流计算之前,需要设置一些潮流计算的参数,包括:1.潮流计算算法:MATPOWER提供了不同的潮流计算算法,如牛顿-拉夫逊法(NR)和次梯度法(SC)等。

可以使用以下命令来设置潮流计算算法:mpopt = mpoption('pf.alg', '<算法名称>')```这里`<算法名称>`可以是`'NR'`或`'SC'`。

2.潮流计算收敛条件:通过设置收敛条件,可以控制潮流计算的准确性和计算时间。

以下是一些设置收敛条件的命令:```matlabmpopt = mpoption(mpopt, 'pf.tol', <收敛容限>)```这里`<收敛容限>`是一个小的正数,用于判断潮流计算是否收敛。

默认值为1e-6```matlabmpopt = mpoption(mpopt, 'pf.nr.max_it', <最大迭代次数>)```这里`<最大迭代次数>`是一个整数,用于限制牛顿-拉夫逊法的最大迭代次数。

默认值为20。

三、执行潮流计算在设置好潮流计算参数之后,可以执行潮流计算。

牛顿拉夫逊法潮流计算

牛顿拉夫逊法潮流计算

牛顿拉夫逊法潮流计算牛顿-拉夫逊法(Newton-Raphson method)是一种用于求解非线性方程的数值方法。

它通过迭代逼近根的方式,将非线性方程转化为一系列的线性方程来求解。

在电力系统中,潮流计算用于确定电力网中节点的电压幅值和相角。

潮流计算是电力系统分析的重要基础,可以用于计算电力系统的潮流分布、功率损耗、节点电压稳定度等参数,为电力系统的规划、运行和控制提供参考依据。

牛顿-拉夫逊法是一种常用的潮流计算方法,它的基本思想是通过不断迭代来逼近电网的潮流分布,直到满足一定的收敛条件。

下面将对牛顿-拉夫逊法的具体步骤进行详细介绍。

首先,我们需要建立电力网络的节点潮流方程,即功率方程。

对于每一个节点i,其节点功率方程可以表示为:Pi - Vi * (sum(Gij * cos(θi - θj)) - sum(Bij * sin(θi -θj))) = 0Qi - Vi * (sum(Gij * sin(θi - θj)) + sum(Bij * cos(θi -θj))) = 0其中,Pi和Qi分别为节点i的有功功率和无功功率,Vi和θi分别为节点i的电压幅值和相角,Gij和Bij分别为节点i和节点j之间的导纳和电纳。

接下来,我们需要对每个节点的电压幅值和相角进行初始化。

一般情况下,可以将电压幅值设置为1,相角设置为0。

然后,我们可以开始进行迭代计算。

在每一轮迭代中,我们需要计算每个节点的雅可比矩阵和功率残差,然后更新电压幅值和相角。

雅可比矩阵可以通过对节点功率方程进行求导得到,具体如下:dPi/dVi = -sum(Vj * (Gij * sin(θi - θj) + Bij * cos(θi - θj)))dPi/dθi = sum(Vj * (Gij * Vi * cos(θi - θj) - Bij * Vi * sin(θi - θj)))dQi/dVi = sum(Vj * (Gij * cos(θi - θj) - Bij * sin(θi - θj)))dQi/dθi = sum(Vj * (Gij * Vi * sin(θi - θj) + Bij * Vi * cos(θi - θj)))功率残差可以通过将节点功率方程代入得到,如下:RPi = Pi - Vi * (sum(Gij * cos(θi - θj)) - sum(Bij *sin(θi - θj)))RQi = Qi - Vi * (sum(Gij * sin(θi - θj)) + sum(Bij *cos(θi - θj)))最后,我们可以使用牛顿-拉夫逊法的迭代公式来更新电压幅值和相角,具体如下:Vi(new) = Vi(old) + ΔViθi(new) = θi(old) + Δθi其中,ΔVi和Δθi分别为通过求解线性方程组得到的电压幅值和相角的增量。

NR法潮流计算

NR法潮流计算

初始值:
( 0) ( 0) ( 0) x1 , x2 xn
按泰勒级数展开,并忽略高次项:
(0) f1 f1 f1 (0) (0) (0) f x x xn 0 1 1 2 x1 0 x2 0 xn 0 (0) f 2 f 2 f 2 (0) (0) ( 0) f2 x1 x2 xn 0 x1 0 x2 0 xn 0 f (0) f n x ( 0 ) f n x ( 0 ) f n x ( 0 ) 0 1 2 n n x x x 1 0 2 0 n 0
n为平衡节点 。
对PQ节点,功率偏差方程:
n Pi Pis U i U j (Gij cos ij Bij sin ij ) j 1 n Qi Qis U i U j (Gij sin ij Bij cos ij ) j 1
m个方程 i 1、 2 m
节点的 U i 和 i 已知,节点的 Pi 和 Qi 待求。
2.6.4 N-R法潮流计算
一、常规的潮流算法
1)高斯-赛德尔法(GS) 内存需求量很少,但计算时间长。 2)牛顿-拉夫逊法(N-R) 具有较高的收敛可靠性和收敛速度,但是需较好的 初始值,且内存占有量大。 3)PQ快速解耦法 计算时间少,内存占用少,但是对病态潮流敏感。
n
Y1n Y2 n Y3n Ynn
自导纳
Yii yi 0 yij
ji
互导纳
1 Yij yij Y ji zij
支路两端节点号
1 2 2 3 3 3 0 0 1 0 1 2
阻抗值
导纳值
0.25j 2*0.25j

9节点电力系统潮流计算

9节点电力系统潮流计算

I力系统分析课稈设计设计题目9节点电力网络潮流计算指导教师_______________________________________ 院(系、部)电气与控制工程学院_________ 专业班级_______________________________________ 学号 __________________________________ 姓名 __________________________________ 日期 ____________________________________电气工程系课程设计标准评分模板课程设计成绩评定表不及格标准:设计内容一项否决制,即5为不及格,整个设计不及格,其他4项否决;优、良、中、及格标准:以设计内容为主体,其他项超过三分之一为评定标准,否则评定为下一等级;如优秀评定,设计内容要符合5,其余九项要有4项符合才能评定为优,否则评定为良好,以此类推。

评定教师签字:最终成绩:目录1PSAS嗽件简介 (1)1.1PSASP平台的主要功能和特点 (1)1.2PSASP勺平台组成 (2)2牛顿拉夫逊潮流计算简介 (3)2.1牛顿—拉夫逊法概要 (3)2.2直角坐标下的牛顿—拉夫逊潮流计算 (5)2.3牛顿—拉夫逊潮流计算的方法 (6)3九节点系统单线图及元件数据 (7)3.1九节点系统单线图 (7)3.2系统各项元件的数据 (8)4潮流计算的结果 (10)4.1潮流计算后的单线图 (10)4.2潮流计算结果输出表格 (10)5结论 (14)电力系统分析课程设计任务书9节点系统单线图如下:基本数据如下:表1母线数据表2交流线数据表3两绕组变压器数据两绕组变压器数据(续表)表4发电机数据负荷数据表5区域数据1 PSASP 软件简介“电力系统分析综合程序” (Power System Analysis Software Package,PSA)SP 是一套历史悠久、功能强大、使用方便的电力系统分析程序,是高度集成和开发具有我国自主知识产权的大型软件包。

牛拉法潮流计算程序(附3机9节点结果对比)

牛拉法潮流计算程序(附3机9节点结果对比)

摘要电力系统潮流计算是研究电力系统稳态运行的一种重要方法,它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态,包括各母线的电压、线路的功率分布以及功率损耗等等。

潮流计算主要用于电网规划和静态安全分析,它可为扩建电力网络,以达到规划周期内所需要的输电能力提供依据;也可以对预想事故进行模拟和分析,校核预想事故下的电力系统安全性。

本文简单介绍了牛顿-拉夫逊潮流计算的原理、模型与算法,然后用具体的实例,利用MATLAB对牛顿-拉夫逊法的算法进行了验证。

关键词:电力系统潮流计算牛顿-拉夫逊法 MATLAB一、牛拉法的数学模型对一个N 节点的电力网路,列写节点电压方程,即I =Y V(1.1)式中,I 为节点注入电流列相量,Y 为节点导纳矩阵,V 为节点电压列相量。

由于异地测量的两个电流缺少时间同步信息,以注入功率替换注入电流作为已知量。

即***1+niij j ij j i i i Y V V I V Q P ∙∙===∑(1.2)其中,Y ij =G ij +j B ij ,带入上式,得到有功功率和无功功率方程 P i =V i V j G ij cos θij +B ij sin θij n j =1 (1.3)Q i =V i V j G ij sin θij −B ij cos θij n j =1(1.4)大部分情况下,已知PQ ,求解V θ。

考虑到电网的功率平衡,至少选择一台发电机来平衡全网有功功率,即至少有一个平衡节点,常选择调频或出线较多的发电机作为平衡节点。

具有无功补偿的母线能保持电压幅值恒定,这类节点可作为PV 节点。

潮流计算中节点分类总结如下:已知电力系统有m 个PQ 节点,r 个PV 节点和1个平衡节点,则可以提取m+r 个有功功率方程和m 个无功功率方程,从而求解出m+r 个θ和m 个V ,其余节点的有功和无功可通过式(1.3)、(1.4)求得,这样就完成了潮流计算。

牛顿拉夫逊法潮流计算

牛顿拉夫逊法潮流计算

牛顿拉夫逊法潮流计算牛顿拉夫逊法是计算电力系统中电流、电压的常用方法之一,也称为牛顿-拉夫逊-里特法或简称为NR法。

资深的电力系统工程师一定对这个方法非常熟悉,但是对于刚刚接触电力系统的人来说,可能会对此感到迷惑。

本文将为大家简单介绍牛顿拉夫逊法的基本步骤,帮助大家更好地理解和使用。

在介绍牛顿拉夫逊法之前,我们需要先了解一些电力系统的基本概念。

电力系统由许多发电厂、输电线路、变电站和用户组成,其中输电线路和变电站是将电能长距离输送和转换的设备。

电力系统中的发电机、负荷和输电线路都具有电阻和电抗,它们之间的复杂相互作用决定了电力系统中的电流和电压。

牛顿拉夫逊法用于计算电力系统节点之间的电流和电压。

节点是指电力系统中有电流和电压变化的点,例如发电机和变电站。

在计算电力系统节点的电流和电压时,我们需要使用一些基本的公式和原理,比如克希荷夫定律和欧姆定律。

下面是牛顿拉夫逊法的基本步骤:1. 确定电力系统中的节点和口纳负荷在计算电力系统的电流和电压之前,我们需要先确定电力系统中所有的节点和负载。

这通常是由电网规划人员完成的。

2. 初始化电力系统中的电流和电压在计算过程中,我们需要先给电力系统中的节点和口纳负荷赋初值。

此时,我们需要假设所有节点的电压相同,即电力系统处于平衡状态。

3. 建立节点电流和电压的方程组建立节点电流和电压的方程组并对其进行求解是计算电力系统电流和电压的关键步骤。

利用克希荷夫定律和欧姆定律,可以得到关于节点电流和电压的一系列方程,这个方程组的解即为电力系统的电流和电压。

4. 更新节点电流和电压求解得到电力系统的电流和电压之后,我们需要更新节点电流和电压的值。

更新后的节点电流和电压将作为下一次计算的初值。

5. 判断计算结果收敛在使用牛顿拉夫逊法计算电力系统电流和电压时,我们需要判断计算结果是否收敛。

如果计算结果没有收敛,即结果不稳定或不趋于一个确定的值,那么我们需要重新建立方程组并进行求解。

电力系统牛拉法潮流计算

电力系统牛拉法潮流计算

电力系统牛拉法潮流计算电力系统的潮流计算是电力系统运行中一项重要的工作,它用来确定电力系统中各节点的电压和功率的分布情况。

牛拉法(Newton-Raphson)方法是一种主要的潮流计算方法,它是基于牛顿迭代法的一种改进方法,可以用来求解非线性方程组,并被广泛应用于电力系统的潮流计算。

牛拉法潮流计算的基本原理是通过不断迭代求解节点电压和相应的功率,直到收敛为止。

具体步骤如下:1.建立潮流计算的数学模型。

电力系统的潮流计算可以被建模为一个非线性方程组,其中未知量为各节点的电压和功率,方程组的解表示系统的潮流分布情况。

2.初始化节点电压。

初始时,可以假设所有节点的电压为1,并根据负荷功率和潮流方向,计算各发电机节点的功率注入。

3.计算节点电压。

利用牛拉法迭代求解非线性方程组。

首先,根据电压相角和幅值的变化情况,更新节点电压;然后,利用更新的节点电压计算各发电机节点的功率注入,以及从节点注入到节点之间的功率传输;最后,根据功率平衡方程计算支路的功率。

4.判断迭代是否收敛。

判断迭代是否收敛的常用方法有两个:一是通过计算节点电压变化量来判断,如果变化量小于一定阈值,则认为计算收敛;二是通过计算功率平衡误差来判断,如果误差小于一定阈值,则认为计算收敛。

5.如果迭代未收敛,返回步骤3;如果迭代收敛,计算结束,得到系统的潮流分布情况。

牛拉法潮流计算的优点是能够处理复杂的非线性方程组,收敛速度快,并且适用于大规模电力系统的计算。

但是,牛拉法潮流计算也存在一些问题,比如可能出现发散情况,需要进行故障处理。

牛拉法潮流计算在电力系统调度和运行中起着重要的作用。

通过潮流计算,可以确保电力系统的稳定运行,优化电力系统的运行方式,提高系统的可靠性和经济性。

总结起来,牛拉法潮流计算是电力系统潮流计算的一种重要方法,通过迭代求解非线性方程组,可以得到电力系统各节点的电压和功率的分布情况。

它在电力系统调度和运行中具有重要的应用价值,可以帮助优化电力系统的运行方式,提高系统的稳定性和经济性。

NR法潮流计算

NR法潮流计算
若修正量 ∆x
(0)
x ( 0 ) 处展开为泰勒级数: 很小, 处展开为泰勒级数: 很小,在
(∆x (0) ) 2 f ( x (0) − ∆x (0) ) = f ( x (0) ) − f ′( x (0) )∆x (0) + f ′′( x (0) ) −L 2!
忽略高次项及更高项: 忽略高次项及更高项:
二、牛顿-拉夫逊法的基本原理 牛顿1、N-R法解单变量非线性方程 、 法解单变量非线性方程 非线性方程: f ( x ) = 0
∆x ( 0 ) 设方程初始解为 x ,初解与真解的偏差为
( 0)
则真解为:
x = x ( 0) − ∆x ( 0)
f (x
(0)
− ∆x ) = 0
(0)
f ( x (0) − ∆x (0) ) = 0
∆δ1 ∆δ ∆δ = 2 L ∆δ n−1
H 雅可比矩阵: 雅可比矩阵: J = J
N L
∆P H ∆Q = J
N ∆δ ∆U U L
∂∆Pi N ij = Uj ∂U j ∂∆Qi Lij = Uj ∂U j
n Pi = U i ∑ U j (Gij cos δ ij + Bij sin δ ij ) j∈i n Qi = U i ∑ U j (Gij sin δ ij − Bij cos δ ij ) j∈i
δ ij = δ i − δ j
二、节点分类
1、PQ节点 、 节点 & 待求。 已知, 节点注入的 Pi + jQi 已知,节点电压 U i 待求。 2、PV节点 2、PV节点 & 已知, 待求。 节点注入的有功 Pi 、 U i已知,节点 δ i 待求。 3、平衡节点 、 已知, 待求。 节点的 U i 和 δ i 已知,节点的 Pi 和 Qi 待求。

3机9节点潮流、短路仿真计算课程设计总结

3机9节点潮流、短路仿真计算课程设计总结

3机9节点潮流、短路仿真计算课程设计总结以3机9节点潮流、短路仿真计算课程设计总结为标题的文章概述:本次课程设计主要涉及到3机9节点潮流和短路仿真计算。

通过对电力系统进行潮流计算和短路仿真,可以了解系统的电压、电流等重要参数,为系统的稳定运行提供参考。

本文将对本次课程设计的过程、结果和总结进行详细介绍。

一、潮流计算潮流计算是电力系统中常用的一种计算方法,用于确定系统中各节点的电压、电流等参数。

在本次课程设计中,我们使用了3台发电机和9个节点的电力系统进行潮流计算。

1.1 数据准备在进行潮流计算之前,需要准备系统的基本数据,包括发电机的有功功率、无功功率和电压,各节点的负荷功率和电压等信息。

通过收集和整理这些数据,我们可以建立电力系统的节点和支路信息。

1.2 潮流计算方法潮流计算可以使用不同的方法,如高斯-赛德尔迭代法、牛顿-拉夫逊法等。

在本次课程设计中,我们选择了高斯-赛德尔迭代法进行潮流计算。

该方法通过迭代计算各节点的电压和电流,直到满足收敛条件为止。

1.3 结果分析经过潮流计算,我们得到了系统中各节点的电压、电流等参数。

通过分析这些结果,我们可以了解系统中的电力流动情况,判断系统是否存在潮流过载、电压偏差等问题,并采取相应的措施进行调整和优化。

二、短路仿真计算短路仿真计算是针对系统发生故障时的一种计算方法,用于确定短路电流的大小和分布情况。

在本次课程设计中,我们使用了相同的3机9节点电力系统进行短路仿真计算。

2.1 短路故障类型短路故障可以分为对称短路和非对称短路两种类型。

对称短路是指系统中的故障电流对称分布,非对称短路则是指故障电流非对称分布。

在本次课程设计中,我们分别考虑了对称短路和非对称短路的情况。

2.2 短路电流计算方法短路电流的计算可以使用不同的方法,如阻抗法、对称分量法等。

在本次课程设计中,我们选择了阻抗法进行短路电流的计算。

该方法通过计算系统中各节点的阻抗和电压,确定短路电流的大小和分布情况。

牛拉法潮流计算

牛拉法潮流计算

牛拉法潮流计算%本程序的功能是用牛拉法进行潮流计算 %原理介绍详见鞠平著《电气工程》%默认数据为鞠平著《电气工程》例8.4所示数据±是支路参数矩阵%第一列和第二列是节点编号。

节点编号由小到大编写%对于含有变压器的支路,第一列为低压侧节点编号,第二列为高压侧节点编号 %第三列为支路的串列阻抗参数,含变压器支路此值为变压器短路电抗 %第四列为支路的对地导纳参数,含变压器支路此值不代入计算 %第五烈为含变压器支路的变压器的变比,变压器非标准电压比%第六列为变压器是否是否含有变压器的参数,其中“1”为含有变压器,“0”为不含有变压器2为节点参数矩阵%第一列为节点注入发电功率参数 %第二列为节点负荷功率参数 %第三列为节点电压参数 %第四列 %第五列%第六列为节点类型参数,“1”为平衡节点,“2”为PQ节点,“3”为PV节点参数%X为节点号和对地参数矩阵 %第一列为节点编号 %第二列为节点对地参数clear; clc;num=input('是否采用默认数据?(1-默认数据;2-手动输入)'); if num==1 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 40.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); end end for i=1:n1Y(i,i)=Y(i,i)+X(i,2); %计及补偿电容电纳 enddisp('导纳矩阵为:');disp(Y); %显示导纳矩阵%初始化OrgS、DetaS OrgS=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; forj=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))); end endendfor i=1:n %对PV节点的处理,注意这时不可再将h初始化为0 ifi~=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))); end end end%创建PVU 用于存储PV节点的初始电压 PVU=zeros(n-h-1,1); t=0; for i=1:nif B2(i,6)==3 t=t+1;PVU(t,1)=B2(i,3); end end%创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量 h=0;for i=1:n %对PQ节点的处理 if i~=isb&B2(i,6)==2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1); TlPiDetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1); TlQi end end t=0;for i=1:n %对PV节点的处理,注意这时不可再将h初始化为0 ifi~=isb&B2(i,6)==3 h=h+1; t=t+1;DetaS(2*h-1,1)=real(B2(i,1))-OrgS(2*h-1,1); TlPiDetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2; TlUi end end % DetaS%创建I,用于存储节点电流参数i=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(2*h,1)*sqrt(-1))/conj(B2(i,3));%conj求共轭end end%创建Jacbi(雅可比矩阵) Jacbi=zeros(2*n-2); h=0; k=0;for i=1:n %对PQ节点的处理 if B2(i,6)==2 h=h+1; forj=1:n if j~=isb k=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; end end end end end k=0;for i=1:n %对PV节点的处理 if B2(i,6)==3 h=h+1; forj=1:n if j~=isb k=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; end end end end enddisp('初始雅可比矩阵为:'); 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)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end endfor i=1:n %对PV节点的处理 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); end end % B2%开始循环********************************************************************** whileabs(max(DetaU))>pr OrgS=zeros(2*n-2,1); h=0;感谢您的阅读,祝您生活愉快。

潮流分析中的N-R法和PQ分解法

潮流分析中的N-R法和PQ分解法

PQ分解法得名! H、L随迭代而变化,如何常数化? 10
进一步简化?
非对角元:Hij Lij UiUj(GijSinij BijCosij )
第二步简化:一般线路两端 ij较小(一般小于 10o~20o),且 Gij Bij,有: cosij 1
Gij sinij Bij cosij
Hij=UiUjBij, i、j=1,2,…,n-1,ij Lij=UiUjBij, i、j=1,2,…,n-r-1,ij
2 -0.234+j0.011
-0.046-j0.136
0.5+j0.093
迭代次数
1
-0.5-j0.029
S1
S1=-0 U=0.985/-0.05
3
S3
精度
潮流流向
P3=0
U=1.1/6.73
29
作业
3-14(只要列出潮流方程、修正方程、B’和B”矩阵, 无须迭代求解)
30
U=1.1/7.0
25
PQ分解法迭代过程(r=2)
U=1.05/0 4 S4
U=0.965/-6.42
S2=-0.001-j0.0002 S2
2
3
1
S1
S3
S1=-0.0006-j0.0002 U=0.985/-0.468
P3=-0.0007
U=1.1/6.77
26
PQ分解法迭代过程(r=3)
U=1.05/0 4 S4
平衡节点1的注入功率:
j0.1 j0.1
2
U
j1
2
U
j1
P G1 P 1 PD1 24,
B Sin(
QG1 Q1 QD1 11.73

三机九节点潮流计算

三机九节点潮流计算

clc;clear;close all;n1=9;n2=9;B1=load('b1.txt','ascii');B2=load('b2.txt','ascii');m=0;r=0;for i=1:n1if B1(i,8)==1m=m+1;endendr=n1-m-1;e=B1(1:n1,2);f=B1(1:n1,3);ps=B1(1:n1,4)-B1(1:n1,6);qs=B1(1:n1,5)-B1(1:n1,7);Y=zeros(n1);for j=1:n2p=B2(j,2);q=B2(j,3);J=zeros(2*n1);Y(p,q)=-1./((B2(j,4)+1i*B2(j,5))*B2(j,7));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B2(j,4)+1i*B2(j,5))+1i*B2(j,6);Y(q,q)=Y(q,q)+1/((B2(j,4)+1i*B2(j,5))*B2(j,7)^2)+1i*B2(j,6);end%导纳矩阵YH1=zeros(n1);N1=zeros(n1);J1=zeros(n1);L1=zeros(n1);G=real(Y);B=imag(Y);dp(1:(n1-1))=1;dq(1:m)=1;dx(1:(2*m+r))=1;while(abs(max(dp)>=1.0e-5)||abs(max(dq)>=1.0e-5)||abs(max(dx(1:(m+r)))>=1.0e-5)||abs(max(dx((m +r+1):(2*m+r)).*e(1:m))>=1.0e-5))for i=1:n1for j=1:n1si(i,j)=sin(f(i)-f(j));co(i,j)=cos(f(i)-f(j));if i==jH1(i,j)=e(i)^2*B(i,i)+qs(i);N1(i,j)=-e(i)^2*G(i,i)-ps(i);J1(i,j)=e(i)^2*G(i,i)-ps(i);L1(i,j)=e(i)^2*B(i,i)-qs(i);elseH1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));N1(i,j)=-e(i)*e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));J1(i,j)=e(i)*e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));L1(i,j)=-e(i)*e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));endendendH=H1(2:n1,2:n1);N=N1(2:n1,r+2:n1);J=J1(r+2:n1,2:n1);L=L1(r+2:n1,r+2:n1);J=[H N;J L];%雅克比矩阵Jfor i=2:n1a(i)=0;b(i)=0;for j=1:n1a(i)=a(i)+e(j)*(G(i,j)*co(i,j)+B(i,j)*si(i,j));if i>=(r+2)b(i)=b(i)+e(j)*(G(i,j)*si(i,j)-B(i,j)*co(i,j));endenddp(i-1)=ps(i)-e(i)*a(i);if i>=(r+2)dq(i-(r+1))=qs(i)-e(i)*b(i);endendF=[dp';dq'];dx=-J\F;f(2:n1)=f(2:n1)+dx(1:(m+r));e1=e((r+2):n1)+dx((m+r+1):(2*m+r)).*e((r+2):n1); e2=e(1:(r+1));e=[e2;e1];enddisp('各节点电压');disp('e=');disp(e');disp('单位:度f=');disp(f'*180/3.1415926)%求出各节点电压幅值及相角a=0;for i=1:n1a=a+conj(Y(1,i))*(e(i)*cos(f(i))-1i*e(i)*sin(f(i))); enda=a*(e(1)*cos(f(1))+1i*e(1)*sin(f(1)));disp('平衡节点P+jQ=');disp(a);disp('输电线路有功无功:');for i=1:n2p=B2(i,2);q=B2(i,3);s1(p,q)=e(p)^2*(-1i*B2(i,6))+(e(p)*cos(f(p))+1i*e(p)*sin(f(p)))*(e(p)*cos(f(p))-1i*e(p)*sin(f(p))-e( q)*cos(f(q))+1i*e(q)*sin(f(q)))*conj(-Y(q,p));s1(q,p)=e(q)^2*(-1i*B2(i,6))+(e(q)*cos(f(q))+1i*e(q)*sin(f(q)))*(e(q)*cos(f(q))-1i*e(q)*sin(f(q))-e( p)*cos(f(p))+1i*e(p)*sin(f(p)))*conj(-Y(q,p));fprintf('S%d%d=',p,q);disp(s1(p,q));fprintf('S%d%d=',q,p);disp(s1(q,p));end%初值计算for j=1:n1s(j)=0;for i=1:n1s(j)=s(j)+conj(Y(j,i))*(e(i)*cos(f(i))-1i*e(i)*sin(f(i)));ends(j)=s(j)*(e(j)*cos(f(j))+1i*e(j)*sin(f(j)));endY1=conj(-s)./(e'.^2);%负荷等值导纳disp('各节点等值导纳');disp('Y1=');disp(Y1);disp('各节点功率');disp('S=');disp(s);Ra=[0 0 0];X_d=[0.0608 0.1198 0.1813];Xq=[0.0969 0.8645 1.2578];%Xq=X_d;%不计凸极效应v=e.*cos(f)+1i*e.*sin(f);for i=1:3I(i)=conj(s(i)/(e(i)*cos(f(i))+1i*e(i)*sin(f(i))));EQ(i)=e(i)*cos(f(i))+1i*e(i)*sin(f(i))+(Ra(i)+1i*Xq(i))*I(i);EQx(i)=real(EQ(i));EQy(i)=imag(EQ(i));delta(i)=atan(EQy(i)/EQx(i));enddisp('各发电机的暂态电动势,功角和输入机械功率初值');disp('delta=');disp(delta*180/3.1415926);for i=1:3Vx(i)=e(i)*cos(f(i));Vy(i)=e(i)*sin(f(i));Vdq=[sin(delta(i)) -cos(delta(i));cos(delta(i)) sin(delta(i))]*[Vx(i);Vy(i)];Vd(i)=Vdq(1);Vq(i)=Vdq(2);Ix(i)=real(I(i));Iy(i)=imag(I(i));Idq=[sin(delta(i)) -cos(delta(i));cos(delta(i)) sin(delta(i))]*[Ix(i);Iy(i)];Id(i)=Idq(1);Iq(i)=Idq(2);E_q(i)=Vq(i)+Ra(i)*Iq(i)+X_d(i)*Id(i);enddisp('E_q=');disp(E_q);for i=1:3Pm(i)=real(s(i))+(Ix(i)^2+Iy(i)^2)*Ra(i); enddisp('Pm=');disp(Pm);。

复杂网络N-R法潮流分析与计算的设计---精品模板

复杂网络N-R法潮流分析与计算的设计---精品模板

电气工程及其自动化专业课程设计复杂网络N-R法潮流分析与计算的设计学生学号:学生姓名:班级:指导教师:起止日期:哈尔滨工程大学自动化学院课程设计报告撰写内容一、设计要求(宋体,小四号字,加黑) 用matlab 编程,N_R 法计算潮流分布 具体要求为:(1)给出程序,并给出注释(2)输出迭代次数,各节点电压,各支路电流 (3)在图中标明功率流向2345T1T2节点数据如下表所示(标幺值)支路及变压器数据 精度要求:0.0001二、设计方案(要求给出详细的设计思路及其必要的论证)(1。

)潮流计算的方法(1)高斯雅克比迭代法(2)高斯-塞得尔法(对初值要求底,迭代次数多)(3)牛顿-拉夫逊法(使用广泛)(4)PQ 快速分解法(提升运算速度)目前广泛应用的潮流计算方法都是基于节点电压法的,以节点导纳矩阵Y作为电力网络的数学模型。

节点电压Ui 和节点注入电流Ii 由节点电压方程YV=I (1)根据S=VI﹡(I﹡为I 的共轭)可得非线性的节点方程YV=I=(S/V)﹡(2)在实际的电力系统中,已知的运行条件不是节点的注入电流,而是负荷和发电机的功率,而且这些功率一般不随节点电压的变化而变化。

由于各节点注入功率与注入电流的关系为Si=Pi+jQi =ViIi﹡,因此可将式(2)改写为Ii=Si/Vi=Pi+jQi/Vi (i= 1, 2,3 ⋯ n)(3)式中,Pi 和Qi 分别为节点i 向网络注入的有功功率和无功功率,当i 为发电机节点时Pi﹥0;当i 为负荷节点时Pi﹤0;当i 为无源节点Pi=0,Qi=0;Vi 和Ii 分别为节点电压相量Vi 和节点注入电流相量Ii 的共轭.式(3)亦即潮流计算的基本方程式,它可以在直角坐标也可以在极坐标上建立2n 个实数形式功率方程式。

发电机Pi、Qi 为正,负荷Pi、Qi 为负。

展开YV=I 为Ii=ΣYijVj=YiiVi+ΣYijVi( i=1 2 3 ⋯n) (4)将式(4)代入式(3),得n 维的非线性复数的电压方程组潮流计算的基本方程为(Pi-jQi)/ Vi= YiiVi+ΣYijVi (i=1, 2,⋯ n) (5)(2。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

如图所示系统,试计算潮流分布,相关数据见《版潮流计算用户手册》P121。

#include<>#include <>float divRe(float b1,float b2,float b3,float b4){float a1r;a1r=(b1*b3+b2*b4)/(b3*b3+b4*b4);return(a1r);}float divIm(float b1,float b2,float b3,float b4){float a1i;a1i=(b2*b3-b1*b4)/(b3*b3+b4*b4);return(a1i);}float mulRe(float b1,float b2,float b3,float b4){float a2r;a2r=b1*b3-b2*b4;return(a2r);}float mulIm(float b1,float b2,float b3,float b4){float a2i;a2i=b2*b3+b1*b4;return(a2i);}float Max(float a[],int n){int i;float max;max= fabs(a[0]);for(i=1;i<n;i++){if(fabs(a[i])>max)max=fabs(a[i]);}return(max);}void main(){int i,j,k,h,km;int T=16;float eps,sumpi1,sumpi2,sumqi1,sumqi2,max,sumir,sumii,I1r,I1i,t,xx,xxx; float pi0[8],qi0[8],detpi[8],detqi[8],Iir0[8],Iii0[8],J0[16][16], detsi[16],detui[16], delta_p[9][9],delta_q[9][9],a[16][32],ni[16][16],H[8][8],N[8][8],J[8][8],L[8][8],ei1[9],fi1[9],sp[9][9],sq[9][9];static float ybr[9][9]={{,,0,0,0,,0,0,0},{,,,0,0,0,0,0,0},{0,,,,0,0,0,0,0},{0,0,,,,0,0,0,0},{0,0,0,,,,0,0,0},{,0,0,0,,,0,0,0},{0},{0},{0}};static float ybi[9][9]={{,,0,0,0,,0,0,},{,,,0,0,0,0,0,0},{0,,,,0,0,16,0,0},{0,0,,,,0,0,0,0},{0,0,0,,,,0,,0},{,0,0,0,,,0,0,0},{0,0,16,0,0,0,-16,0,0},{0,0,0,0,,0,0,,0},{,0,0,0,0,0,0,0,}};static float yd[9][9]={{0,,0,0,0,,0,0,0},{,0,,0,0,0,0,0,0},{0,,0,,0,0,0,0,0},{0,0,,0,,0,0,0,0},{0,0,0,,0,,0,0,0},{,0,0,0,,0,0,0,0},{0},{0},{0}} ;float ei0[9]={,,,,,,,,};float fi0[9]={};float pi[9]={0,,0,,0,,,,0};float qi[9]={0,,0,,0,,0,0,0};h=0;km=15;eps=;do{h+=1;printf("\nNow The %dth times\n",h);for(i=0;i<8;i++){printf("ei0[%d]=%f\t",i,ei0[i]);printf("fi0[%d]=%f\t",i,fi0[i]);}for(i=0;i<8;i++){printf("pi[%d]=%f\t",i,pi[i]);printf("qi[%d]=%f\t",i,qi[i]);}sumpi2=0;sumqi2=0;for(i=0;i<8;i++){for(j=0;j<9;j++){sumpi1=(ei0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])+fi0[i]*(ybr[i][j]*fi0[j]+ybi [i][j]*ei0[j]));sumpi2+=sumpi1;}pi0[i]=sumpi2;printf("pi0[%d]=%f\t",i,pi0[i]);sumpi2=0;}for(i=0;i<8;i++){for(j=0;j<9;j++){sumqi1=(fi0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])-ei0[i]*(ybr[i][j]*fi0[j]+ybi [i][j]*ei0[j]));sumqi2+=sumqi1;}qi0[i]=sumqi2;printf("qi0[%d]=%f\t",i,qi0[i]);sumqi2=0;}for(i=0;i<8;i++){detpi[i]=pi[i]-pi0[i];detqi[i]=qi[i]-qi0[i];if(i==6||i==7){qi0[i]=ei0[i]*ei0[i]+fi0[i]*fi0[i];detqi[i]=[i];}printf("detpi[%d]=%f\t",i,detpi[i]);printf("detqi[%d]=%f\t",i,detqi[i]);}//***************************************************************************** ********************//节点的注入电流表达式for(i=0;i<8;i++){ Iii0[i]=0;Iir0[i]=0;}for(i=0;i<8;i++){for(j=0;j<9;j++){Iir0[i]+=ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j];Iii0[i]+=ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j];}}//***************************************************************************** ****************//求解NHJL矩阵for(i=0;i<8;i++){for(j=0;j<8;j++){if(i==j){if(i==6||i==7){H[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]+Iii0[i];N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]+Iir0[i];J[i][j]=2*fi0[i];L[i][j]=2*ei0[i];}else{H[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]+Iii0[i];N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]+Iir0[i];J[i][j]=-ybr[i][j]*ei0[j]-ybi[i][i]*fi0[j]+Iir0[i];L[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]-Iii0[i];}}else{if(i==6||i==7){H[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j];N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j];J[i][j]=0;L[i][j]=0;}else{ H[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j];N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j];J[i][j]=-ybi[i][j]*fi0[j]-ybr[i][j]*ei0[j];L[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j];}}}}//形成jacobian矩阵for(i=0;i<16;i++)for(j=0;j<16;j++){if(i%2==0&&j%2==0) J0[i][j]=H[i/2][j/2];else if(i%2==0&&j%2!=0) J0[i][j]=N[i/2][(j-1)/2];else if(i%2!=0&&j%2==0) J0[i][j]=J[(i-1)/2][j/2];else J0[i][j]=L[(i-1)/2][(j-1)/2];}//for(i=0;i<16;i++)//for(j=0;j<16;j++)// printf("J0[%d][%d]=%f\t",i,j,J0[i][j]);//***************************************************************************** *********//求detui//***************************************************************************** ********for(i=0;i<16;i++){if(i%2==0) detsi[i]=detpi[i/2];else detsi[i]=detqi[(i-1)/2];}//将detp和detq用一个数组表示for(i=0;i<T;i++)for(j=0;j<(2*T);j++){ if (j<16)a[i][j]=J0[i][j];else if (j==T+i)a[i][j]=;elsea[i][j]=;}for(i=0;i<T;i++){for(k=0;k<T;k++){if(k!=i){t=a[k][i]/a[i][i];for(j=0;j<(2*T);j++){xx=a[i][j]*t;a[k][j]=a[k][j]-xx;}}}}for(i=0;i<T;i++){t=a[i][i];for(j=0;j<(2*T);j++)a[i][j]=a[i][j]/t;}for(i=0;i<T;i++)for(j=0;j<T;j++)ni[i][j]=a[i][j+T];/*printf("逆矩阵为:\n");for (i=0;i<T;i++){for (j=0;j<T;j++)printf("%10.3f",ni[i][j]);printf("\n");} */xxx=;for(i=0;i<T;i++){xxx=;for(j=0;j<T;j++){xxx=xxx+ni[i][j]*detsi[j];}detui[i]=xxx;}//检测detui满足要求与否max=Max(detui,16);printf("max=%f\n",max);//****************************************************************** //算第n次迭代后的ufor(i=0;i<T;i++){if(i%2==0){ fi1[i/2]=fi0[i/2]+detui[i];}else{ei1[(i-1)/2]=ei0[(i-1)/2]+detui[i];}}//*********************************************for(i=0;i<8;i++)//下一次迭代赋初值{ei0[i]=ei1[i];fi0[i]=fi1[i];}for(i=0;i<8;i++){printf("ei0[%d]=%f\t",i,ei0[i]);printf("fi0[%d]=%f\t",i,fi0[i]);}for(i=0;i<8;i++){pi[i]=detpi[i]+pi0[i];qi[i]=detqi[i]+qi0[i];}}while(max>eps&&h<km);printf("All do %d times\n",h);sumir=0;sumii=0;//**********************************************for(i=0;i<9;i++)//平衡节点功率计算{I1r=mulRe(ybr[8][i],-ybi[8][i],ei0[i],-fi0[i]); I1i=mulIm(ybr[8][i],-ybi[8][i],ei0[i],-fi0[i]); sumir+=I1r;sumii+=I1i;}pi[8]=mulRe(ei0[8],fi0[8],sumir,sumii);qi[8]=mulIm(ei0[8],fi0[8],sumir,sumii);printf("S9=%f+j%f\n",pi[8],qi[8]);sumpi1=0;sumpi2=0;sumqi1=0;sumqi2=0;for(i=0;i<9;i++){for(j=0;j<9;j++)if(i!=j&&ybi[i][j]!=0){sumpi1=mulRe(ei0[i],fi0[i],,yd[i][j]);sumqi1=mulIm(ei0[i],fi0[i],,yd[i][j]);sumpi2=mulRe(ei0[i]-ei0[j],fi0[i]-fi0[j],-ybr[i][j],-ybi[i][j]);sumqi2=mulIm(ei0[i]-ei0[j],fi0[i]-fi0[j],-ybr[i][j],-ybi[i][j]);sumpi1+=sumpi2;sumqi1+=sumqi2;sp[i][j]=mulRe(ei0[i],fi0[i],sumpi1,-sumqi1);sq[i][j]=mulIm(ei0[i],fi0[i],sumpi1,-sumqi1);printf("S%d=%f+j%f\n",(i+1)*10+j+1,sp[i][j],sq[i][j]);}}for(i=0;i<9;i++){for(j=0;j<9;j++)if(j !=i&&ybi[i][j]!=0){delta_p[i][j]=sp[i][j]+sp[j][i];delta_q[i][j]=sq[i][j]+sq[j][i];if( (delta_p[i][j] !=0) || (delta_q[i][j] !=0)){printf("dS%d=%f+j%f\n",(i+1)*10+j+1,delta_p[i][j],delta_q[i][j]);}}}ei1[8]=ei0[8];fi1[8]=fi0[8];for(i=0;i<9;i++){printf("u%d=%f<%f\n",i+1,sqrt(ei1[i]*ei1[i]+fi1[i]*fi1[i]),atan(fi1[i]/ei1[i])*180/;} }。

相关文档
最新文档