毕业论文-潮流计算

毕业设计•目录-PAGEii-TOC\o"1-4"\h\z目录一、 电力系统潮流计算概述 1二、 节点导纳矩阵 11. 节点导纳矩阵及其各元素的物理意义 12. 节点导纳矩阵的特点 13. 算法推导 23.1 励磁支路 23.2 线路支路 23.3 变压器支路 34. 程序实现方法和技巧 34.1 变量说明 34.2 支路参数设置技巧 45. 程序框图 56. 矩阵输出格式控制 6三、 潮流计算的原理及实现方法 61. 牛顿—拉夫逊法（直角、极坐标） 61.1 牛顿—拉夫逊法简介 61.2 算法推导 8a. 直角坐标 8b. 极坐标 111.3 程序实现方法和技巧 121.4 直角坐标和极坐标的比较 131.5 程序框图 142. 高斯—塞得尔法 152.1 高斯—塞得尔法简介 152.2 算法推导 152.3 程序框图 173. P-Q分解法 183.1 P-Q分解法简介 183.2 算法推导 183.3 程序框图 214. 节点功率、功率分布及网损计算 224.1 节点功率 224.2 支路功率 22a. 变压器支路 23b. 线路支路 24c. 励磁支路 244.3 网络损耗 24四、 主程序结构图 25五、 程序清单 26六、 程序计算实例 471. 实例内容 472. 程序初值 483. 运行结果 483.1 节点导纳矩阵 483.2 牛顿极坐标 493.3 牛顿直角坐标 493.4 高斯塞得尔 503.5 P-Q分解法 514. 运行界面 52七、 毕业设计的成果、经验和不足 561. 程序设计的优点 562. 程序中的不足 563. 毕业设计成果 574. 毕业设计体会 57参考书目 58毕业设计•电力系统潮流计算-PAGE60-电力系统潮流计算电力系统潮流计算概述潮流计算是电力系统中最基本的计算，它在给定电网结构、参数及决定电力系统运行状态的边界条件的情况下，通过计算来确定电力系统的运行状态。如各母线上的电压（幅值和相角）、网络中的功率分布及功率损耗等。潮流计算是电力系统运行、规划、以及安全性、可靠性分析和优化的基础，也是各种电磁瞬时过程和机电瞬时过程分析的基础和出发点。潮流计算一般分为在线计算和离线计算。在这里主要是进行离线潮流计算的问题。潮流计算的方法通常有阻抗法和导纳法。随着计算机的应用发展，采用稀疏矩阵技术和节点优化编号技术使牛顿－拉夫逊法得到了非常广泛的应用，成为潮流计算中广泛采用的优秀算法。导纳法也就因其稀疏性的特点逐步取代了阻抗法。70年代中期，Stott在大量实践的基础上提出了快速分解潮流计算方法（P_Q分解法），使潮流计算在内存容量和计算速度上得到了大大的提高，可用于在线计算，是目前最常用的一种方法。节点导纳矩阵节点导纳矩阵及其各元素的物理意义由欧姆定律可知：，也即，其中系数矩阵Y即称为节点导纳矩阵。对于n个节点的电力网络来说，方程组的系数矩阵：（2-SEQ（2-\*ARABIC1）显然可知其中各元素的物理意义：导纳矩阵中第i列对角元素Yii，即节点i的自导纳，在数值上等于节点i加单位电压，其它节点都接地时，节点i向电力网络注入的电流。导纳矩阵中第i列非对角元素Yij，即节点i与节点j间的互导纳，在数值上等于节点i加单位电压，其它节点都接地时，节点j向电力网络注入的电流。节点导纳矩阵的特点导纳矩阵是稀疏矩阵，矩阵的阶数等于电力网的节点数。各行非对角元素的个数等于对应节点所连的不接地支路数。导纳矩阵是一个对称的方阵，非对角元素Yij=Yji。对角线元素即自导纳等于i节点所有支路导纳值的总和，即Yii=∑yii。非对角线元素即互导纳等于i、j节点间支路导纳的负值，即Yij=-yij。算法推导电力网中常见的有三种类型的支路：励磁支路、线路支路、变压器支路。在节点导纳矩阵程序中，采用逐条支路追加法。而对于各种支路处理方法略有不同，下面就分别作以推导。励磁支路如图1所示，励磁支路的另一节点是接地点0，因此当追加励磁支路时，只需要修正节点i的自导纳。（2-SEQ（2-\*ARABIC2）线路支路如图2所示，当追加线路支路时，节点i、节点j的自导纳以及它们之间的互导纳都需要修改。（2-SEQ（2-\*ARABIC3）（2-SEQ（2-\*ARABIC4）（2-SEQ（2-\*ARABIC5）图SEQ图\*ARABIC1励磁支路图SEQ图\*ARABIC2线路支路变压器支路变压器一共有四种模型，在这里采用如图3所示的模型，而其它的三种情况均可以化作此模型（详见后一节）。对于图3的模型，有：（2-SEQ（2-\*ARABIC6）（2-SEQ（2-\*ARABIC7）（2-SEQ（2-\*ARABIC8）化简后：（2-SEQ（2-\*ARABIC9）（2-SEQ（2-\*ARABIC10）图SEQ图\*ARABIC3变压器支路程序实现方法和技巧变量说明YG、YB分别存放节点导纳矩阵的实部和虚部；m存放循环数；（i、j、r、x、k）存放支路参数i、j分别存放支路的节点号，如果是接地点则存为0；r、x分别存放支路参数的实部和虚部；k存放变比（变压器支路）或接地容纳（线路支路）；支路参数设置技巧（1）如何区别支路类型当其为励磁支路时，令接地端的节点号为0，这样判断就是励磁支路；当其为变压器支路时，令靠近变压器的一端为负值，这样判断则是变压器支路；因而判断就是线路支路。（2）如何解决变压器四种模型的统一对于图4所示的模型，令，当程序判断时，节点号互换；对于图5所示的模型，令，当程序判断时，作；对于图6所示的模型，同时做上述两步；这样一来就可以把模型统一到图3所对应的计算公式。图SEQ图\*ARABIC4图SEQ图\*ARABIC5图SEQ图\*ARABIC6（3）如何让线路支路与变压器支路用同一套计算公式把线路支路的对地容纳存入变量中，当程序处理线路支路时，首先用以下两个公式（2-11）和（2-12）把容纳加入节点的自导纳中，然后置，这样即把线路支路看作变比的变压器支路。（2-SEQ（2-\*ARABIC11）（2-SEQ（2-\*ARABIC12）程序框图矩阵输出格式控制由于节点导纳矩阵的输出中，有实部和虚部之分、又有正负号，要把字符j输出到虚部的正负号之前，而且要把整个矩阵用一屏输完，就需要进行一定的处理。对于实部照直输出即可，对于虚部则先判断其正负，输出或，然后对其取绝对值输出。如果阶数比较大，还可以用printf(“%.NUMf”,变量名)；来控制输出数据的小数点位数。其中NUM是希望输出的小数点位数。潮流计算的原理及实现方法如上所述，节点电流与电压之间的关系可以通过节点方程式来表述：（3-SEQ（3-\*ARABIC1）上式写成展开的形式：（3-SEQ（3-\*ARABIC2）为了求解潮流问题，必须利用节点功率与电流之间的关系：（3-SEQ（3-\*ARABIC3）其中分别为节点向网络注入的有功功率及无功功率，当点为负荷节点时，本身应带负号。将式（3-3）代入式（3-2），可得到：（3-SEQ（3-\*ARABIC4）式（3-4）中有n个非线性复数方程式，是潮流计算问题的基本方程式，对这个方程式的不同应用和处理，就形成了不同的潮流程序。总体上可分为：牛顿—拉夫逊法、高斯—塞得尔法、P—Q分解法。以下对各种跌代方法作一一说明。牛顿—拉夫逊法（直角、极坐标）牛顿—拉夫逊法简介牛顿—拉夫逊法是把非线性方程式的求解过程变成反复对相应的线性方程式的求解过程，从几何意义上看实质是切线法、是一种逐步线性化的方法，这是牛顿—拉夫逊法的核心所在。设有变量的非线性联立方程组：（3-SEQ（3-\*ARABIC5）其近似为，设近似解与精确解分别相差，则如下的方程式成立：（3-SEQ（3-\*ARABIC6）上式中任何一式都可按泰勒级数展开，并省略二次以上高次项，可得：（3-SEQ（3-\*ARABIC7）或简写为：（3-SEQ（3-\*ARABIC8）式中，称为函数的雅可比矩阵；将代入式（3-7），可得中，然后求解，从而经过第一次迭代后的新值。再将求得的代入，如此循环，直到。（这里为预先给定的满足一定精度要求的正整数）算法推导直角坐标基本方程式：（3-SEQ（3-\*ARABIC9）将式（3-2）代入式（3-9），并展开功率方程式得：（3-SEQ（3-\*ARABIC10）节点电压以直角坐标表示时，有、，将实部和虚部分列，可得：（3-SEQ（3-\*ARABIC11）对于PQ节点，Pis，Qis是给定的，因而有：（3-SEQ（3-\*ARABIC12）对于PV节点，Pis，Vis是给定的，因而有：（3-SEQ（3-\*ARABIC13）其中令：（3-SEQ（3-\*ARABIC14）则式（3-12）和（3-13）可简化为：（3-SEQ（3-\*ARABIC15）（3-SEQ（3-\*ARABIC16）建立修正方程式：（3-SEQ（3-\*ARABIC17）其中用式（3-15）求出，用式（3-16）求出；而雅可比矩阵各元素按下式求得：时，（3-SEQ（3-\*ARABIC18）（3-SEQ（3-\*ARABIC19）（3-SEQ（3-\*ARABIC20）（3-SEQ（3-\*ARABIC21）时，（3-SEQ（3-\*ARABIC22）（3-SEQ（3-\*ARABIC23）（3-SEQ（3-\*ARABIC24）（3-SEQ（3-\*ARABIC25）（3-SEQ（3-\*ARABIC26）（3-SEQ（3-\*ARABIC27）为了求解修正方程式，必须对雅可比矩阵作一定的处理，可以用高斯消去法、三角分解法、因子表法，而在本程序中采用直接求逆的方法。求逆后有：（3-SEQ（3-\*ARABIC28）通过式（3-28）求出，再修正各节点电压，即：（3-SEQ（3-\*ARABIC29）这样节点电压就离精确值更近了一步。再运用新的电压值进行下一次迭代，直到满足精度要求。极坐标如上所述，当节点电压以极坐标表示时，有（3-SEQ（3-\*ARABIC30）把式（3-30）代入式（3-10），再将实部、虚部分开，得功率方程式：（3-SEQ（3-\*ARABIC31）式中（3-SEQ（3-\*ARABIC32）在计算中通常将式（3-31）写成：（3-SEQ（3-\*ARABIC33）在建立修正方程式时，由于PV节点的电压幅值是给定的，不在作为变量，同时，该点不能预先给定无功功率Qis，这样，方程式中也就失去了约束作用。因此，在迭代过程中应取消与PV节点有关的无功功率方程式。设节点总数为n，PV节点数为r，则：（3-SEQ（3-\*ARABIC34）为了程序上处理方便，把修正方程式排列成下列形式：（3-SEQ（3-\*ARABIC35）式中，雅可比矩阵个元素用下式计算：时（3-SEQ（3-\*ARABIC36）（3-SEQ（3-\*ARABIC37）（3-SEQ（3-\*ARABIC38）（3-SEQ（3-\*ARABIC39）时，（3-SEQ（3-\*ARABIC40）（3-SEQ（3-\*ARABIC41）（3-SEQ（3-\*ARABIC42）（3-SEQ（3-\*ARABIC43）同上，用直接求逆法处理雅可比矩阵，解出、，再修正各节点电压：（3-SEQ（3-\*ARABIC44）如此循环直到满足精度要求为止。程序实现方法和技巧直角坐标法：式（3-15）、（3-16）及式（3-22）～（3-25）中都含有节点注入电流的分量、，因而求解、及雅可比矩阵中对角线元素、、、的主要运算集中在求、上。节点注入电流的分量、只与行导纳矩阵及相应节点的电压分量有关（见式3-14），因此即可先运算得到、，从而减小运算量、提高计算精度。另外，在直角坐标和极坐标中，当时，都有、，如此就可以只作两个子函数H（i，j）、N（i，j）来进行调用，这样也减小了程序的运算量，同时增加程序的可读性。注：为了程序计算的方便，把PV节点的电压存放于无功功率Q的数组中。直角坐标和极坐标的比较表格SEQ表格\*ARABIC1计算方法比较项目直角坐标法极坐标法节点电压节点功率方程修正量方程雅可比元素见式（3-18）～(3-27)见式（3-36）～(3-43)矩阵的维数2(n-1)×2(n-1)[2(n-1)-r]×[2(n-1)-r]矩阵的处理直接求逆直接求逆节点电压修正方程程序框图高斯—塞得尔法高斯—塞得尔法简介与牛顿法不同的是高斯—塞得尔法把迭代计算所求得的最新值立即用于计算下一个变量的新值而不是等到这一轮迭代结束之后。这就是高斯—塞得尔法思想的核心所在。通常高斯—塞得尔法迭代的次数很多，为了加速收敛，可以引入加速因子。算法推导如果取，则有功方程式：（3-SEQ（3-\*ARABIC45）把实部和虚部分开：（3-SEQ（3-\*ARABIC46）令：（3-SEQ（3-\*ARABIC47）因此式（3-46）可简化为：（3-SEQ（3-\*ARABIC48）从式（3-4）可以解出：（3-SEQ（3-\*ARABIC49）将其展开，并写成高—塞法的迭代格式：（3-SEQ（3-\*ARABIC50）最后引入加速因子：（3-SEQ（3-\*ARABIC51）这就是高斯迭代的基本公式。A．对于PQ节点，节点功率是给定的，因此可以直接应用式（3-50）和式（3-51）进行迭代计算。B．对于PV节点，电压幅值是给定的，而无功功率又需要在迭代过程中逐次地求出，因此PV节点必须作以下几项计算。修正节点电压（3-SEQ（3-\*ARABIC52）用式（3-48）计算节点无功功率无功功率越限检查检验（3-SEQ（3-\*ARABIC53）作完上述三项计算后，就可以用式（3-50）和式（3-51）计算节点电压的新值。C．平衡节点不必进行迭代另外，高斯—塞得尔法的收敛判据是：（3-SEQ（3-\*ARABIC54）程序框图P-Q分解法P-Q分解法简介P-Q分解法（改进牛顿法）是从改进和简化牛顿法潮流程序的基础上提出来的。它的基本思想是：把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，把有功功率和无功功率迭代分开来进行。即：（3-SEQ（3-\*ARABIC55）这样由于把2n阶的方程组变成了二个n阶的方程组，因而计算量和内存方面都有所改善。P-Q分解法的另一个简化是将式（3-55）中的系数矩阵化作在迭代过程中不变的对称矩阵。迭代过程中不需要反复求解系数矩阵。以上两点就是P-Q分解法与其它方法的最明显的区别。算法推导如上所述，有：（3-SEQ（3-\*ARABIC56）由于线路两端电压的相角差不太大（一般10°～20°），因此近似认为、，且近似认为。作此近似后，可以将式（3-56）系数矩阵中的元素表示式从式（3-36）、（3-39）、（3-40）、（3-43）简化为：（3-SEQ（3-\*ARABIC57）将式（3-57）写成矩阵的形式有：（3-SEQ（3-\*ARABIC58）再将式（3-58）代入式（3-56）中，并利用矩阵相乘，得：（3-SEQ（3-\*ARABIC59）（3-SEQ（3-\*ARABIC60）进一步化简后可以得到：（3-SEQ（3-\*ARABIC61）（3-SEQ（3-\*ARABIC62）式（3-61）和式（3-62）也可以简写成：（3-SEQ（3-\*ARABIC63）（3-SEQ（3-\*ARABIC64）***注意：这里的与的阶数是不同的。以上两式就是P-Q分解法的修正方程式，其中系数矩阵与只不过是系统导纳矩阵的虚部，因而是对称矩阵，而且在迭代过程中维持不变。它们与功率误差方程式（3-65）：（3-SEQ（3-\*ARABIC65）构成了P-Q分解法迭代过程中基本计算公式。至于其迭代过程，与牛顿极坐标法相类似，这里就不再熬述了。程序框图节点功率、功率分布及网损计算节点功率要计算的节点功率有PV节点的无功功率，平衡节点的有功和无功功率。电压取直角坐标时，用下式计算：（3-SEQ（3-\*ARABIC66）电压取极坐标时，则有：（3-SEQ（3-\*ARABIC67）支路功率如图7所示，支路功率基本方程式：（3-SEQ（3-\*ARABIC68）也即：（3-SEQ（3-\*ARABIC69）图SEQ图\*ARABIC7图SEQ图\*ARABIC8下面对各种支路分别讨论：变压器支路在这里仍然采用节点导纳程序中的变压器模型，如图8所示，即有：（3-SEQ（3-\*ARABIC70）又由于式（2-6），所以上式可以写成：（3-SEQ（3-\*ARABIC71）当电压取直角坐标时，并把实部、虚部分开，有：（3-SEQ（3-\*ARABIC72）（3-SEQ（3-\*ARABIC73）当电压取极坐标时，并把实部、虚部分开，有：（3-SEQ（3-\*ARABIC74）（3-SEQ（3-\*ARABIC75）而又考虑到式（3-36）和（3-37），、，则式（3-74）和式（3-75）可进一步化简为：（3-SEQ（3-\*ARABIC76）（3-SEQ（3-\*ARABIC77）如此一来，就只需多次调用H1(inti,intj)和N1(inti,intj)两个子函数，这样就使程序得到了很大的简化，也增强了程序的可读性。线路支路对于线路支路，可以把接地导纳（数值存放于变量k中）先行加入线路功率之中，即：（3-SEQ（3-\*ARABIC78）然后把置变量k＝1。这样就可以把剩下的线路支路看作是变比为1的变压器支路来处理。励磁支路励磁支路的功率用下式计算：（3-SEQ（3-\*ARABIC79）网络损耗网络损耗用下式计算：（3-SEQ（3-\*ARABIC80）主程序结构图程序清单#include"stdio.h"#include"math.h"#defineN5/*定义节点数*/#definePI3.14159staticfloatYG[N+1][N+1],YB[N+1][N+1];staticfloatJ[2*N+1][2*N+1],JPQ[2*N+1];staticfloatP[N+1],Q[N+1],DP[N+1],DQ[N+1];staticfloatV[N+1],Sita[N+1],DV[N+1],DSita[N+1];staticfloate[N+1],f[N+1],De[N+1],Df[N+1];staticfloatP_x[N+1][N+1],Q_x[N+1][N+1];/*计算P、Q的子函数*/floatP1(inti){ intj; floatP1=0,tempa=0; for(j=1;j<=N;j++) { if(i==j)continue; tempa+=(V[j]*(YG[i][j]*cos(Sita[i]-Sita[j])+YB[i][j]*sin(Sita[i]-Sita[j]))); }P1=tempa*V[i]; return(P1);}floatQ1(inti){ intj; floatQ1=0,tempb=0; for(j=1;j<=N;j++) { if(i==j)continue; tempb+=(V[j]*(YG[i][j]*sin(Sita[i]-Sita[j])-YB[i][j]*cos(Sita[i]-Sita[j]))); }Q1=V[i]*tempb; return(Q1);}/*计算雅可比矩阵元素的子函数*/floatH1(inti,intj){ floatH1; H1=-V[i]*V[j]*(YG[i][j]*sin(Sita[i]-Sita[j])-YB[i][j]*cos(Sita[i]-Sita[j])); return(H1);}floatN1(inti,intj){ floatN1; N1=-V[i]*V[j]*(YG[i][j]*cos(Sita[i]-Sita[j])+YB[i][j]*sin(Sita[i]-Sita[j])); return(N1);}floatH2(inti,intj){ floatH2=0; H2=(YB[i][j]*e[i]-YG[i][j]*f[i]); return(H2);}floatN2(inti,intj){ floatN2=0;N2=(-1)*(YG[i][j]*e[i]+YB[i][j]*f[i]);return(N2);}/*主程序*/main(){ inti,j,m,a,t=0,tt=0;intflag1,flag[N+1];intKp,Kq; floatr,x,k,b;floatp1,q1,sita1,v1;floate1,f1,DV_2;floatA[N+1],B[N+1];floatE[N+1],F[N+1];floatDP_V[N+1],DQ_V[N+1],DSita_V[N+1];floattemp1,temp2,temp3,temp;floattempa1=0.0,tempa2=0.0;floatQ_max=5.0,Q_min=-5.0;/*定义无功越限检查的上、下限*/floatB1[N+1][N+1],B2[N+1][N+1]; floatmax1=0.0,max2=0.0;floatdP,dQ;/*定义网损变量*/floatspeed;/*定义加速因子*/ FILE*fp1,*fp2;intR;/*定义PV节点数变量*/intselect;/*定义程序的选择变量*//*打开支路参数原始数据文件*/ if((fp1=fopen("lqp","r"))==NULL) { printf("Thefilecannotopen"); exit(0); }/*读取数据并计算节点导纳矩阵*/ while(!feof(fp1)){ fscanf(fp1,"%d%d%f%f%f\n",&i,&j,&r,&x,&k); if(i*j==0) { if(i==0) { a=i;i=j;j=a; } YG[i][i]+=r; YB[i][i]-=x; } else { if(i*j>0) { YB[i][i]+=k; YB[j][j]+=k; k=1; } else { if(j<0) { a=i;i=j;j=a; } if(k<0) k=-1/k; } i=fabs(i); b=r*r+x*x; YG[i][j]=YG[j][i]=(-1)*r/(k*b); YB[i][j]=YB[j][i]=x/(k*b);YG[j][j]+=r/b; YB[j][j]-=x/b; YG[i][i]+=r/(k*k*b); YB[i][i]-=x/(k*k*b);}}fclose(fp1);start:clrscr();/*计算方法选择*/printf("电力系统潮流计算\n\n\n\n\n");printf("请选择要使用的计算方法:\n\n");printf("1牛顿--拉夫逊法（极坐标）\n\n");printf("2牛顿--拉夫逊法（直角坐标）\n\n");printf("3高斯--塞得尔法\n\n");printf("4P--Q分解法\n\n\n");printf("pleaseenterthenumber:");scanf("%d",&select);switch(select){case1:printf("牛顿--拉夫逊法（极坐标）运行结果:\n");gotolabel1;case2:printf("牛顿--拉夫逊法（直角坐标）运行结果:\n");gotolabel2;case3:{printf("请输入加速因子的值(1～1.5)(bestfor1.44):");scanf("%f",&speed);printf("高斯--塞得尔法运行结果:\n");gotolabel3;}case4:printf("P--Q分解法运行结果:\n");gotolabel4;}gotoend_all;/*牛顿极坐标程序入口*/label1:/*打开节点原始数据文件并读取*/ if((fp2=fopen("lqp2","r"))==NULL) { printf("Thefilecann'tbeopened"); exit(0); }R=0; while(!feof(fp2)) { fscanf(fp2,"%d%d%f%f%f%f",&i,&flag1,&p1,&q1,&v1,&sita1);flag[i]=flag1;if(flag[i]==-1)R++; P[i]=p1; Q[i]=q1; V[i]=v1; Sita[i]=sita1; } fclose(fp2);/*进行主循环*/t=0;do{/*求功率误差DP、DQ*/for(i=1;i<=N-1;i++){ DP[i]=P[i]-V[i]*V[i]*YG[i][i]-P1(i); JPQ[2*i-1]=DP[i];}for(i=1;i<=N-R-1;i++){ DQ[i]=Q[i]+V[i]*V[i]*YB[i][i]-Q1(i); JPQ[2*i]=DQ[i];}/*计算雅可比矩阵各元素*/for(i=1;i<=N-R-1;i++)for(j=1;j<=N-1;j++){ if(i!=j) { if(j<=N-R-1) { J[2*i-1][2*j-1]=H1(i,j); J[2*i-1][2*j]=N1(i,j); J[2*i][2*j-1]=(-1)*N1(i,j); J[2*i][2*j]=H1(i,j); } else { J[2*i-1][2*j-1]=H1(i,j); J[2*i][2*j-1]=(-1)*N1(i,j); }} else { J[2*i-1][2*j-1]=Q1(i); J[2*i-1][2*j]=-2*V[i]*V[i]*YG[i][i]-P1(i); J[2*i][2*j-1]=-P1(i); J[2*i][2*j]=2*V[i]*V[i]*YB[i][i]-Q1(i); }}for(i=N-R;i<=N-1;i++)for(j=1;j<=N-1;j++){ if(i!=j) { if(j<=N-R-1) { J[2*i-1][2*j-1]=H1(i,j); J[2*i-1][2*j]=N1(i,j); } else J[2*i-1][2*j-1]=H1(i,j); } else J[2*i-1][2*j-1]=Q1(i);}/*求逆雅可比矩阵*/for(m=1;m<=2*(N-1)-R;m++){for(i=1;i<=2*(N-1)-R;i++) for(j=1;j<=2*(N-1)-R;j++) if(i!=m&&j!=m) J[i][j]=J[i][j]-J[i][m]*J[m][j]/J[m][m];for(j=1;j<=2*(N-1)-R;j++) if(j!=m) J[m][j]=(-1)*J[m][j]/J[m][m];for(i=1;i<=2*(N-1)-R;i++) if(i!=m) J[i][m]=(-1)*J[i][m]/J[m][m]; J[m][m]=(-1)/J[m][m];}for(i=1;i<=2*(N-1)-R;i++)for(j=1;j<=2*(N-1)-R;j++) J[i][j]=(-1)*J[i][j]; /*解修正方程式*/for(i=1;i<=2*(N-1)-R;i++) { inttt=(i+1)/2; floattemp=0; for(j=1;j<=2*(N-1)-R;j++) { temp+=J[i][j]*JPQ[j]; if(i%2==1) DSita[tt]=-temp; else DV[tt]=-temp*V[tt]; } } max1=0; for(i=1;i<N;i++) { if(max1<fabs(DSita[i]))max1=fabs(DSita[i]); if(max1<fabs(DV[i]))max1=fabs(DV[i]); } /*修正各节点电压*/ for(i=1;i<N;i++) { Sita[i]=Sita[i]+DSita[i]; V[i]=V[i]+DV[i]; }t++;}while(max1>1e-5);gotocalculate;/*牛顿直角坐标程序入口*/label2:if((fp2=fopen("lqp2","r"))==NULL){ printf("Cannotopen"); exit(0);}while(!feof(fp2)){ fscanf(fp2,"%d%d%f%f%f%f\n",&i,&flag1,&p1,&q1,&e1,&f1); flag[i]=flag1; P[i]=p1; Q[i]=q1; if(flag[i]==-1)/*把PV节点的电压存放到无功数组中*/ Q[i]=e1; e[i]=e1; f[i]=f1;}fclose(fp2);/*进行主循环*/t=0;do{ for(i=1;i<N;i++) { floattempa=0,tempb=0; for(j=1;j<=N;j++){ tempa+=YG[i][j]*e[j]-YB[i][j]*f[j];tempb+=YG[i][j]*f[j]+YB[i][j]*e[j];} A[i]=tempa; B[i]=tempb;/*处理PQ节点*/if(flag[i]==1){ DP[i]=P[i]-(e[i]*A[i]+f[i]*B[i]);DQ[i]=Q[i]-(f[i]*A[i]-e[i]*B[i]); for(j=1;j<N;j++){ if(j==i) { J[2*i-1][2*j-1]=H2(i,j)-B[i]; J[2*i-1][2*j]=N2(i,j)-A[i]; J[2*i][2*j-1]=(-1)*N2(i,j)-A[i]; J[2*i][2*j]=H2(i,j)+B[i];}else { J[2*i-1][2*j-1]=H2(i,j); J[2*i-1][2*j]=N2(i,j); J[2*i][2*j-1]=(-1)*N2(i,j); J[2*i][2*j]=H2(i,j);}} }/*处理PV节点*/ if(flag[i]==-1) { DP[i]=P[i]-(e[i]*A[i]+f[i]*B[i]); DQ[i]=pow(Q[i],2)-pow(e[i],2)-pow(f[i],2); for(j=1;j<N;j++){ if(j==i) { J[2*i-1][2*j-1]=H2(i,j)-B[i]; J[2*i-1][2*j]=N2(i,j)-A[i];J[2*i][2*j-1]=(-2)*f[i];J[2*i][2*j]=(-2)*e[i];}else { J[2*i-1][2*j-1]=H2(i,j); J[2*i-1][2*j]=N2(i,j);J[2*i][2*j-1]=0;J[2*i][2*j]=0;}} } }/*求逆雅可比矩阵*/ for(m=1;m<2*N-1;m++){for(i=1;i<2*N-1;i++)for(j=1;j<2*N-1;j++)if(i!=m&&j!=m)J[i][j]=J[i][j]-J[i][m]*J[m][j]/J[m][m];for(j=1;j<2*N-1;j++)if(j!=m)J[m][j]=(-1)*J[m][j]/J[m][m];for(i=1;i<2*N-1;i++)if(i!=m)J[i][m]=(-1)*J[i][m]/J[m][m];J[m][m]=(-1)/J[m][m];}for(i=1;i<2*N-1;i++) for(j=1;j<2*N-1;j++)J[i][j]=(-1)*J[i][j];/*解修正方程式*/ for(i=1;i<N;i++) { Df[i]=0; De[i]=0; } for(i=1;i<N;i++)for(j=1;j<N;j++){Df[i]+=J[2*i-1][2*j-1]*DP[j]+J[2*i-1][2*j]*DQ[j];De[i]+=J[2*i][2*j-1]*DP[j]+J[2*i][2*j]*DQ[j];}/*修正各节点电压*/ max1=0;for(i=1;i<N;i++){e[i]=e[i]-De[i];f[i]=f[i]-Df[i]; if(max1<fabs(De[i]))max1=fabs(De[i]);if(max1<fabs(Df[i]))max1=fabs(Df[i]);} t++;}while(max1>1e-5);gotocalculate2;/*高斯--塞得尔法程序入口*/label3:if((fp2=fopen("lqp2","r"))==NULL){ printf("Cannotopen"); exit(0);}while(!feof(fp2)){ fscanf(fp2,"%d%d%f%f%f%f\n",&i,&flag1,&p1,&q1,&e1,&f1); flag[i]=flag1; P[i]=p1; Q[i]=q1; e[i]=e1; f[i]=f1;V[i]=sqrt(e[i]*e[i]+f[i]*f[i]);/*存放PV节点的电压幅值（常值）*/}fclose(fp2);/*进行主循环*/t=0;do{ for(i=1;i<=N;i++) {if(flag[i]==0)continue;/*如果是平衡节点则不迭代*/E[i]=e[i];/*让E、F保存电压的前一次值*/F[i]=f[i]; tempa1=0.0,tempa2=0.0; for(j=1;j<=N;j++){if(j==i)continue; tempa1+=YG[i][j]*e[j]-YB[i][j]*f[j];tempa2+=YG[i][j]*f[j]+YB[i][j]*e[j];} A[i]=tempa1; B[i]=tempa2;/*对PV节点作特殊处理*/if(flag[i]==-1){/*修正PV节点的电压（保证PV节点的幅值不变）*/temp=0;temp=sqrt(e[i]*e[i]+f[i]*f[i]);e[i]=e[i]*V[i]/temp;f[i]=f[i]*V[i]/temp;/*计算PV节点的无功功率Q*/Q[i]=(-1)*YB[i][i]*(e[i]*e[i]+f[i]*f[i])+f[i]*A[i]-e[i]*B[i];/*无功功率越限检查*/if(Q[i]>Q_max)Q[i]=Q_max;if(Q[i]<Q_min)Q[i]=Q_min;}/*计算节点电压新值*/temp1=0;temp2=0;temp3=0;temp1=(P[i]*e[i]+Q[i]*f[i])/(e[i]*e[i]+f[i]*f[i])-A[i];temp2=(P[i]*f[i]-Q[i]*e[i])/(e[i]*e[i]+f[i]*f[i])-B[i];temp3=(YG[i][i]*YG[i][i]+YB[i][i]*YB[i][i]);e[i]=(YG[i][i]*temp1+YB[i][i]*temp2)/temp3;f[i]=(YG[i][i]*temp2-YB[i][i]*temp1)/temp3;}max1=0;for(i=1;i<=N-1;i++){if(fabs(e[i]-E[i])>max1)max1=fabs(e[i]-E[i]);if(fabs(f[i]-F[i])>max1)max1=fabs(f[i]-F[i]);}/*用加速因子speed修正节点电压*/for(i=1;i<=N-1;i++) {e[i]=E[i]+speed*(e[i]-E[i]); f[i]=F[i]+speed*(f[i]-F[i]); }t++;}while(max1>1e-5);gotocalculate2;/*P--Q分解法入口*/label4:if((fp2=fopen("lqp2","r"))==NULL){ printf("Cannotopen"); exit(0);}R=0;while(!feof(fp2)){ fscanf(fp2,"%d%d%f%f%f%f\n",&i,&flag1,&p1,&q1,&v1,&sita1); flag[i]=flag1;if(flag[i]==-1)R++; P[i]=p1; Q[i]=q1; V[i]=v1; Sita[i]=sita1;}fclose(fp2);/*求出系数矩阵B1、B2*/for(a=1,i=1;a<=N;a++,i++){for(b=1,j=1;b<=N;b++,j++){B1[i][j]=YB[a][b];if(flag[b]==0)j--;}if(flag[a]==0)i--;}for(a=1,i=1;a<=N;a++,i++){for(b=1,j=1;b<=N;b++,j++){B2[i][j]=YB[a][b];if(flag[b]==0||flag[b]==-1)j--;}if(flag[a]==0||flag[a]==-1)i--;}/*求逆系数矩阵B1、B2*/for(m=1;m<=N-1;m++){for(i=1;i<=N-1;i++)for(j=1;j<=N-1;j++)if(i!=m&&j!=m)B1[i][j]=B1[i][j]-B1[i][m]*B1[m][j]/B1[m][m];for(j=1;j<=N-1;j++)if(j!=m)B1[m][j]=(-1)*B1[m][j]/B1[m][m];for(i=1;i<=N-1;i++)if(i!=m)B1[i][m]=(-1)*B1[i][m]/B1[m][m];B1[m][m]=(-1)/B1[m][m];}for(i=1;i<=N-1;i++)for(j=1;j<=N-1;j++)B1[i][j]=(-1)*B1[i][j];for(m=1;m<=N-R-1;m++){for(i=1;i<=N-R-1;i++)for(j=1;j<=N-R-1;j++)if(i!=m&&j!=m)B2[i][j]=B2[i][j]-B2[i][m]*B2[m][j]/B2[m][m];for(j=1;j<=N-R-1;j++)if(j!=m)B2[m][j]=(-1)*B2[m][j]/B2[m][m];for(i=1;i<=N-R-1;i++)if(i!=m)B2[i][m]=(-1)*B2[i][m]/B2[m][m];B2[m][m]=(-1)/B2[m][m];}for(i=1;i<=N-R-1;i++)for(j=1;j<=N-R-1;j++)B2[i][j]=(-1)*B2[i][j];/*进行主循环*/t=0;Kp=1;Kq=1;/*置收敛标志为1*/do{part1:/*进入相角迭代部分*/max1=0;for(i=1;i<=4;i++){floattemp1=0;for(j=1;j<=5;j++)temp1+=V[j]*(YG[i][j]*cos(Sita[i]-Sita[j])+YB[i][j]*s