直角坐标系下牛顿法潮流计算

直角坐标系下牛顿法潮流计算

1电力系统潮流计算

潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角)、网络中的功率分布及功率损耗等。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性.可靠性和经济性。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。

2节点导纳矩阵的形成

在图1（a）的简单电力系统中，若略去变压器的励磁功率和线路电容，负荷用阻抗表示，便可以得到一个有5个节点（包括零电位点）和7条支路的等值网络，如图1（b）所示。将接于节点1和4的电势源和阻抗的串联组合变换成等值的电流源和导纳的并联组合，变得到图1（c）的等值网络，其中I1y10E1和I4y40E4分别称为节点1和4的注入电流源。

1(a)234y101y122y24y233y20(b)4y34y40Ė1y244y343y40Ė41y122y23İ1y´10y20İ4(c)图1 电力系统及其网络

以零电位点作为计算节点电压的参考点，根据基尔霍夫定律，可以写出4个节点的电流平衡方程如下：

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

y12(U2U1)y20U2y23(U2U3)y24(U2U4)0y23(U3U2)y34(U3U4)0y24(U4U2)y34(U4U3)y40U4I4 （2-1） y10U1y12(U1U2)I1上述方程组经过整理可以写成

Y11U1Y12U2I1Y21U1Y22U2Y23U3Y24U40Y32U2Y33U3Y34U40Y42U2Y43U3Y44U4I4 （2-2）

式中，

Y11y10y12；

；

Y22y20y23y24y12；

；

Y33y23y34；

Y44y40y24y34；

Y12Y21y12Y23Y32y23Y24Y42y24；

Y34Y43y34。

一般的，对于有n个节点的网络，可以列写n个节点方程

Y11U1Y12U2Y21U1Y22U2Yn1U1Yn2U2也可以用矩阵写成

Y1nUnI1Y2nUnI2YnnUnIn （2-3）

Y11Y12Y1nU1I1YYY2nU2I22122Yn1Yn2YnnUnIn （2-4）

或缩写为

YUI （2-5）

矩阵Y称为节点导纳矩阵。它的对角线元素Yii称为节点i的自导纳，其值等于接于节点i的所有支路导纳之和。非对角线元素

Yij称为节点i、j 间的互导纳，它等于直接

Y0接于节点i、j间的支路导纳的负值。若节点i、j间不存在直接支路，则有ij。由此可知节点导纳矩阵是一个稀疏的对称矩阵。

3牛顿-拉夫逊法潮流计算

3.1牛顿-拉夫逊法的基本原理

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

牛顿—拉夫逊法（Newton—Raphson法）是求解非线性方程代数方程组的有效迭代计算方法。在牛顿—拉夫逊法的每一次迭代过程中，对非线性方程通过线性化处理逐步近似。下面以单变量加以说明。

设有单变量非线性方程

f(x)0 (3-1)

求解此方程时。先给出解的近似值x将满足方程，即

(0)它与真解的误差为

x(0)，则

xxx(0)(0)f(xx)0 (3-2)

将(3-8)式左边的函数在x

(0)(0)(0)附近展成泰勒级数，于是便得

f(x(0)x)f(x)(0)(0)f(x'(0))x(0)f(x''(0))(x)2!(0)2......f(n)(x)(0)(x)n!(0)n....

（3-3）

式中

f'(x),……

(0)(0)f(n)(0)(0)(x)分别为函数f(x)在x处的一阶导数，….，n阶导数。

如果差值x3-9便简化为

很小，3-9式右端x f((0)的二次及以上阶次的各项均可略去。于是，

'x(0)x)f(x)(0)(0)f(x(0))x(0)＝0 (3-4)

这是对于变量的修正量

x(0)的现行方程式，亦称修正方程式。解此方程可得修正量

x用所求的x(0)(0)f(x)(0)f(x(0)'(0) (3-5)

)去修正近似解，变得

x(1)xx(0)x(0)f(x)(0)f(x'(0) (3-6)

)(0)由于3-10是略去高次项的简化式，因此所解出的修正量x(1)也只是近似值。修

正后的近似解x同真解仍然有误差。但是，这样的迭代计算可以反复进行下去，迭代计算的通式是

(k1)(k)xxf(x)(k)f'(x)(k) (3-7)

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

迭代过程的收敛判据为

f(x)1(k) (3-8)

或

式中

x(k)2 (3-9)

，12为预先给定的小正数。

这种解法的几何意义可以从图3－1得到说明。函数y＝f(x)为图中的曲线。f(x)

＝0的解相当于曲线与x轴的交点。如果第k次迭代中得到

x(k)，则过

(k),(k)f((k))点作一切线，此切线同x轴的交点便确定了下一个近似值(k1)。

yxxx由此可见，牛顿－拉夫逊法实质上就是切线法，是一种逐步线性化的方法。

应用牛顿法求解多变量非线性方程组3-1时，假定已给出各变量的初值

(0)，x1x(0)2….

x(0)n，令

x(0)1，

x(0)2，…..

x(0)n分别为各变量的修正量，使其满足

方程3-2即

f((0)(0),(0)(0),....,(0)(0))0x1x2x2xnxn1x1((0)(0),(0)(0),....,(0)(0))0f2x1x1x2x2xnxn

......(x(0)x(0),x(0)x(0),....,x(0)x(0))0fn1122nn

(3-10)

(0)x1将上式中的n个多元函数在初始值附近分别展成泰勒级数,并略去含有，(0)x2，……，xn二次及以上阶次的各项，便得

(0)

fff111(0)(0)(0)(0)(0)f(,,...,)...x(0)0|||xxxxx12n12nx10x20xn01fff111(0)(0)(0)(0)(0),x,...,x)xx...x(0)0f2(x|||12n12n. x10x20xn0......((0),(0),...,(0))f1(0)f1(0)...f1(0)0|0x1|0x2|0xnxfnx1x2nx1x2xn

(3-11)

方程式3-17也可以写成矩阵形式

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

 ffff1|0x1((0),x(0),...,x(0))1x12nf2(x(0),x(0),...,x(0))|0212nx1.........(x(0),x(0),...,x(0))nf12nn|0x1(0)f|xf|x12220.........0...f|xn20(0)fxnn...|xn0...(0)f|x(0)xf(0)|x x...10n1202n (3-12)

方程式3-18是对于修正量x1，x2，……，xn 的线性方程组,称为牛顿法的

(0)(0)(0)x修正方程式.利用高斯消去法或三角分解法可以解出修正量1，x2，……，xn。

(0)然后对初始近似值进行修正

x(1)x(0)x(0) (i=1,2,….,n) (3-13)

iii如此反复迭代，在进行k＋1次迭代时，从求解修正方程式

f1|kx1(x(k),x(k),...,x(k))112nf2(k)(k)(k)(,,...,x)|k2x1x2nx1.........(x(k),x(k),...,x(k))nf12nn|kx1 ffff|xf|x1222k.........k...f|xn2k(k)fxnn...|xnk...f|x(k)xf(k)|x x...1kn12k2n (3-14)

得到修正量x(k)，x(k)，x(k)，并对各变量进行修正

12n

x(k1)x(k)x(k)iii (i=1,2,…,n) (3-15)

式3-20和3-21也可以缩写为

Fx(k)J(k)x(k) (3-16)

和

x(k1)x(k)x(k) (3-17)

式中的X和X分别是由n个变量和修正量组成的n维列向量；F(X)是由n个多元函数组成的n维列项量；J是n阶方阵，称为雅可比矩阵，它的第i、j个元素

Jijfxii是第n个函数

f(x,x,...,x,)对第j个变量xi12nj的偏导数；上角标(k)表示J阵的每一

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

个元素都在点

maxf(x(k)1x(k)2...,x(k)n)处取值。迭代过程一直到满足收敛判据

i,,,f(x(k),x(k),...,x(k))i12n1 (3-18)或

maxx(k)2 (3-19)

i为止。

1和

2为预先给定的小正数。

将牛顿－拉夫逊法用于潮流计算，要求将潮流方程写成形如方程式3-1的形式。

由于节点电压可以采用不同的坐标系表示，牛顿－拉夫逊法潮流计算也将相应的采用不同的计算公式。

Qifi(GijejBijfj)ei(GijfjBijej)j1j1nn3.2节点电压用直角坐标表示是的牛顿-拉夫逊法潮流计算

采用直角坐标时，节点电压可表示为

ejf Viii导纳矩阵元素则表示为

YijGijjBij

将上述表示式代入SiPijQiUiIiUi得

Piei(GijejBijfj)fi(GijfjBijej)j1j1nnYUijjinj的右端，展开并分出实部和虚部，便

--完整版学习资料分享----

（11-45）

-----WORD格式--可编辑--专业资料-----

假定系统中的第1，2，3···，m号节点为PQ节点，第i个节点的给定功率设为Pis和Qis，对对该节点可列写方程

PiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1nnPiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1nn

（i=1,2，···，m） （11-46） 假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一个节点可以列写方程

PiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1

Vi2Vis2Vi2Vis2(ei2fi2)0nn （i=m+1，m+2，···，n-1） （11-47） 第n号节点为平衡点，其电压Vnenjfn是给定的，故不参加迭代。

式（11-46）和式（11-47）总共包含了2（n-1）个方程，待求的变量有e1,f1,...,en1,fn1也是2（n-1）个。我们还可看到，方程（11-46）和式（11-47）已经具备了方程组（11-34）

的形式。因此，不难写出如下的修正方程式

WJV （11-48） 式中

22WP1Q1...PmQmPm1Vm1...Pn1Vn1

TTVe1f1...emfmem1fm1...en1fn1

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

P1P1P1P1P1P1P1P1ef1emfmem1fm1en1fn11Q1Q1Q1Q1Q1Q1Q1Q1efefefef1mmm1m1n1n11PPPPPPPPmmmmmmmme1f1emfmem1fm1en1fn1QmQmQmQmQmQmQmQmef1emfmem1fm1en1fn11 JPm1Pm1Pm1Pm1Pm1Pm1Pm1Pm1ef1emfmem1fm1en1fn1122222222Vm1Vm1Vm1Vm1Vm1Vm1Vm1Vm1efefefef11mmm1m1n1n1PPPPPPPPn1n1n1n1n1n1n1n1e1f1emfmem1fm1en1fn122222222Vn1Vn1Vn1Vn1Vn1Vn1Vn1Vn1ef1emfmem1fm1en1fn11上述方程中雅克比矩阵的各元素，可以对式（11-46）和式（11-47）求偏导数获得。当

ij时

PiQi(GijeiBijfi)ejfjPiQiBijeiGijfi)  （11-49） fjejV2iV2i0ejfj当ji时

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

nPi(GikekBikfk)GiieiBiifieik1nPi(GikfkBikek)BiieiGiififik1nQi(GikfkBikek)BiieiGiifieik1  （11-50） nQifi(GikekBikfk)GiieiBiifik1V2ie2eiiV2if2fii

修正方程式（11-48）还可以写成分块矩阵的形式

W1J11J12J1.n1 W2JJV122J2.n121V2 Wn1Jn1.1JJn1.2n1.n1Vn1式中，Wi和Vi都是二维列向量；Jij是22介方阵。

Veiif

i对于PQ节点

WPiiQ iPiPi JejfjijQ Qiiejfj对于PV节点

WPiiV2

iPiPi JejfjijV2 V2iiejfj从表达式（1-7）~（1-11）可以看到，雅克比矩阵有以下特点：

--完整版学习资料分享----

11-51） 11-52） 11-53） （（（-----WORD格式--可编辑--专业资料-----

（1）雅克比矩阵各元素都是节点电压的函数，它们的数值将在迭代过程中不断的改变。

（2）雅克比矩阵的子块Jij中的元素的表达式只用到导纳矩阵中的对应元素Yij。若

Yij0，则必有Jij0。因此，式（1-9）式中分块形式的雅克比矩阵同节点导纳矩阵一样稀疏，修正方程的求解同样可以用稀疏矩阵的求解技巧。

（3）无论在式（1-6）或式（1-9）中雅克比矩阵的元素或子块都不具有对称性。 用牛顿-拉夫逊法计算潮流的流程框图示于图11-27。首先要输入网络的原始数据以及各节点的给定值并形成节点导纳矩阵。输入节点电压初值ei(0)和fi(0)，置迭代计数k=0。然后开始进入牛顿法的迭代过程。在进行第k+1次迭代时，其计算步骤如下：

（1）按上一次迭代计算出的节点电压值e(k)和f(k)，利用式（11-46）和式（11-47）计算各类节点的不平衡量Pi(k)、Qi(k)和Vi2(k)。

（2）按条件（11-43）校验收敛，即

maxPi(k)、Qi(k)、Vi2(k)<

如果收敛，迭代到此结束，转入计算各线路潮流和平衡节点的功率，并打印输出计算结果。不收敛则继续计算。

（3）利用式（11-49）和式（11-50）计算雅克比矩阵的各元素。 （4）解修正方程式（11-48）求节点电压的修正量ei(k)和fi(k)。 （5）修正各节点的电压

ei(k1)ei(k)ei(k),fi(k1)fi(k)fi(k)

（6）迭代计数加1，返回第一步继续迭代过程。

迭代结束后，还要计算出平衡节点的功率和网络中的功率分布。 输电线路功率的计算公式如下：

SijPijjQijViIijViyi0Vi(ViVj)yij

•2•••

3.3牛顿-拉夫逊法潮流计算程序框图

输入原始数据 --完整版学习资料分享---- 形成节点导纳矩阵 (0)(0)-----WORD格式--可编辑--专业资料-----

图2 牛顿-拉夫逊法潮流计算程序框图

4 C语言实现直角坐标表示的牛顿拉夫逊的计算

4.1计算节点导纳矩阵

该电力系统的等值电路如图所示

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

图3电力系统等值电路

根据节点方程和等值电路图，计算节点导纳矩阵的C语言子程序如下。

G[1][1]=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m);

B[1][1]=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r* z14r+z14m*z14m)+y140+y120;

G[2][2]=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);

B[2][2]=-z12m/(z12r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z24m)+y240+y120; G[3][3]=z13r/(z13r*z13r+z13m*z13m); B[3][3]=-z13m/(z13r*z13r+z13m*z13m);

G[4][4]=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);

B[4][4]=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140; G[1][2]=G[2][1]=-z12r/(z12r*z12r+z12m*z12m); B[1][2]=B[2][1]=z12m/(z12r*z12r+z12m*z12m); G[1][3]=G[3][1]=-k*z13r/(z13r*z13r+z13m*z13m); B[1][3]=B[3][1]=k*z13m/(z13r*z13r+z13m*z13m); G[1][4]=G[4][1]=-z14r/(z14r*z14r+z14m*z14m); B[1][4]=B[4][1]=z14m/(z14r*z14r+z14m*z14m); G[2][3]=G[3][2]=0.0; B[2][3]=B[3][2]=0.0;

G[2][4]=G[4][2]=-z24r/(z24r*z24r+z24m*z24m);

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

B[2][4]=B[4][2]=z24m/(z24r*z24r+z24m*z24m); G[3][4]=G[4][3]=0.0; B[3][4]=B[4][3]=0.0; for(i=1;i<5;i++) {for(j=1;j<5;j++)

{printf(\"%f+%fj\printf(\" \"); }

printf(\"\\n\");

4.2计算各节点不平衡量

按上一次迭代计算出的节点电压值e(k)和f(k)（当k=0时即为给定值），利用式（1-2和式（1-3）计算各类节点的不平衡量Pi(k)、Qi(k)和Vi2(k)。

4(0)4(0)(0)(0)(0)(0)PPPB1jej)1sP1se1(G1jejB1jfj)f1(G1jfj j1j14(0)4 (0)(0)(0)(0)(0)(0)(0)Q1Q1sQ1Q1sf1(G1jejB1jfj)e1(G1jfjB1jej)j1j1(0)1(0)1

同样的可以算出

P2(0)P2sP2(0 )(0)(0)Q2Q2sQ2

其部分C语言程序如下：

P3(0)P3sP3(0)V2(0)3V3sV3s2(0)2printf(\"迭代次数k1=%d\\n\for (i=1;i<5;i++) {float a=0,b=0; for(j=1;j<5;j++)

{a+=G[i][j]*e[j]-B[i][j]*f[j]; b+=G[i][j]*f[j]+B[i][j]*e[j]; }

P[i]=Ps[i]-(e[i]*a+f[i]*b);//计算有功功率的增量 Q[i]=Qs[i]-(f[i]*a-e[i]*b);//计算无功功率的增量

V32=V3s*V3s-e[3]*e[3]; }

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

printf(\"有功功率增量P[1]=%f\ printf(\" ,\");

printf(\"有功功率增量P[2]=%f\ printf(\" ,\");

printf(\"有功功率增量P[3]=%f\

printf(\"无功功率增量Q[1]=%f\ printf(\" ,\");

printf(\"无功功率增量Q[2]=%f\ printf(\" ,\");

printf(\"电压增量V32=%f\

printf(\"\\n\");

计算后校验收敛条件

(k)(k)2(k) maxPi、Qi、Vi< （2-0）

根据给定的容许误差105，按公式（2-0）校验是否收敛。根据C语言编程如下。

max=fabs(P[1])>fabs(P[2])?fabs(P[1]):fabs(P[2]); max=max>fabs(P[3])?max:fabs(P[3]); max=max>fabs(Q[1])?max:fabs(Q[1]); max=max>fabs(Q[2])?max:fabs(Q[2]); max=max>fabs(V32)?max:fabs(V32); printf(\"max=%f\\n\

//******************************************** while (max>0.00001)

经计算，各节点的不平衡量都未满足收敛条件，于是继续计算。

4.3计算雅克比矩阵

按公式（2-1）和（2-2）计算雅克比矩阵的各元素，方程中雅克比矩阵的各元素，可

以对式（1-2）和式（1-3）求偏导数获得。

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

当ij时， 当ji时，

计算雅可比矩阵的C程序如下：

for(j=1;j<4;j++)

PiQi(GijeiBijfi)ejfjPiQiBijeiGijfi)fjejV2iV2i0ejfj{if(1==j)

nPi(GikekBikfk)GiieiBiifieik1nPi(GikfkBikek)BiieiGiififik1nQi(GikfkBikek)BiieiGiifieik1nQi(GikekBikfk)GiieiBiififik12Vi2eiei2Vi2fifi{float c=0,d=0; int m;

for(m=1;m<5;m++)

{c+=G[1][m]*e[m]-B[1][m]*f[m]; d+=G[1][m]*f[m]+B[1][m]*e[m]; }

J[1*N-1][j*N-1]=-c-G[1][j]*e[1]-B[1][j]*f[1]; J[1*N-1][j*N]=-d+B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N-1]=d+B[1][j]*e[1]-G[1][j]*f[1];

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

J[1*N][j*N]=-c+G[1][j]*e[1]+B[1][j]*f[1]; } else

{J[1*N-1][j*N-1]=-G[1][j]*e[1]-B[1][j]*f[1]; J[1*N][j*N]=G[1][j]*e[1]-B[1][j]*f[1]; J[1*N-1][j*N]=B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N-1]=B[1][j]*e[1]-G[1][j]*f[1]; } }

4.4解修正方程

修正方程式

WJV

式中

22WPQ...PQPV...PV11mmm1m1n1n1

T 修正各节点电压后，迭代计数加1，然后返回继续迭代过程。 4.5计算平衡节点功率

迭代结束后，还要算出平衡点的功率和网络中的功率分布。输电线路功率的计算公式如下。

SijPijjQijViIijViyi0Vi(ViVj)yij

..2.....

计算平衡点功率和网络中功率分布的程序如下。

p12=-2*e[1]*f[1]*y120-(e[1]*(e[1]-e[2])-f[1]*(f[1]-f[2]))*G[1][2]+(e[1]*(f[1]-f[2])+f[1]*(e[1]-e[2]))*B[1][2];

q12=-(e[1]*e[1]-f[1]*f[1])*y120+(e[1]*(e[1]-e[2])+f[1]*(f[1]-f[2]))*B[1][2]+(e[1]*(f[1]-f[2]) +f[1]*(e[1]-e[2]))*G[1][2];

5程序调试结果

程序调试结果如图所示

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

图4程序调试结果图

总 结

本设计采用直角坐标形式的牛顿—拉夫逊法作常规潮流计算。P-Q分解法利用了电力系统的一些特有的运行特点，对牛顿—拉夫逊法作了简化，可提高计算速度，但

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

较难理解，牛顿—拉夫逊法的雅克比矩阵在每次迭代过程中都有变化，需要重新形成和求解，这占用了计算的大部分时间，成为牛顿—拉夫逊计算速度不能提高的主要因素，但收敛性好，物理概念也较为清晰。应用C语言进行编程层次清晰，便于按模块化方式组织，易于调试和维护，C语言的表现能力和处理能力极强，极大的提高了计算效率，同时也培养了我们快速的掌握应用程序的能力和解决实际问题的能力。

通过本次设计，我们更好的理解和掌握电力系统潮流计算的原理和方法，也更加系统的理解C程序设计语言的编程技巧。今后的研究方向仍围绕着如何改善传统潮流计算方法以及提高计算机潮流计算的可靠性、收敛性以及速度和灵活性。

参考文献

[1] 何仰赞等.电力系统分析上册[M]．武汉：华中科技大学出版社. [2] 何仰赞等.电力系统分析下册[M]．武汉：华中科技大学出版社. [3] 陈珩.电力系统稳态分析. 水利电力出版社,1994

[4] 熊静琪. 计算机控制技术. 电子工业出版社，ISBN 7- 5053-8706- 5. [5]张伯明,陈寿孙.高等电力网络分析[M].北京：清华大学出版社,1996.

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

主程序如下：

#include #include

#include #define N 2 #define M 7 main()

附录

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

{float z12r,z12m,y120,z13r,z13m,k,z14r,z14m,y140,z24r,z24m,y240,G[5][5],B[5][5],J[7][7]; float e[5]={0,1,1,1.1,1.05},f[5]={0},P[5],Q[5],Ps[5]={0,-0.30,-0.55,0.5},xe[4],xf[4]; float Qs[5]={0,-0.18,-0.13},V3s=1.1,V4S=1.05; float V32,max,P4,Q4; float a1=0,b1=0; int i,j,n,s,k1=0;

float L[M][M]={0},U[M][M]={0},sigma1,sigma2,b[M],y[M],x[M];

float p12,p13,p14,p21,p24,p31,p41,p42,q12,q13,q14,q21,q24,q31,q41,q42; printf(\"请输入z12的实部和虚部\\n\"); //输入电路中的阻抗 scanf(\"%f%f\

printf(\"请输入z13的实部和虚部\\n\"); scanf(\"%f%f\

printf(\"请输入z14的实部和虚部\\n\"); scanf(\"%f%f\

printf(\"请输入z24的实部和虚部\\n\"); scanf(\"%f%f\printf(\"请输入y120的值\\n\"); scanf(\"%f\

printf(\"请输入y140的值\\n\"); scanf(\"%f\

printf(\"请输入y240的值\\n\"); scanf(\"%f\

printf(\"请输入变比k的值\\n\"); scanf(\"%f\

//****************计算导纳矩阵*******************

G[1][1]=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m);

B[1][1]=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z14r+z14m*z14m)+y140+y120;

G[2][2]=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);

B[2][2]=-z12m/(z12r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z24m)+y240+y120; G[3][3]=z13r/(z13r*z13r+z13m*z13m); B[3][3]=-z13m/(z13r*z13r+z13m*z13m);

G[4][4]=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);

B[4][4]=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140; G[1][2]=G[2][1]=-z12r/(z12r*z12r+z12m*z12m); B[1][2]=B[2][1]=z12m/(z12r*z12r+z12m*z12m); G[1][3]=G[3][1]=-k*z13r/(z13r*z13r+z13m*z13m); B[1][3]=B[3][1]=k*z13m/(z13r*z13r+z13m*z13m); G[1][4]=G[4][1]=-z14r/(z14r*z14r+z14m*z14m); B[1][4]=B[4][1]=z14m/(z14r*z14r+z14m*z14m); G[2][3]=G[3][2]=0.0; B[2][3]=B[3][2]=0.0;

G[2][4]=G[4][2]=-z24r/(z24r*z24r+z24m*z24m); B[2][4]=B[4][2]=z24m/(z24r*z24r+z24m*z24m); G[3][4]=G[4][3]=0.0;

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

B[3][4]=B[4][3]=0.0; for(i=1;i<5;i++) {for(j=1;j<5;j++)

{printf(\"%f+%fj\printf(\" \"); }

printf(\"\\n\");//形成节点导纳矩阵

//******************************************* }

printf(\"\\n\");

//******************************************** //计算各节点不平衡量 loop1:

printf(\"迭代次数k1=%d\\n\for (i=1;i<5;i++) {float a=0,b=0; for(j=1;j<5;j++)

{a+=G[i][j]*e[j]-B[i][j]*f[j]; b+=G[i][j]*f[j]+B[i][j]*e[j]; }

P[i]=Ps[i]-(e[i]*a+f[i]*b);//计算有功功率的增量 Q[i]=Qs[i]-(f[i]*a-e[i]*b);//计算无功功率的增量

V32=V3s*V3s-e[3]*e[3]; }

printf(\"有功功率增量P[1]=%f\ printf(\" ,\");

printf(\"有功功率增量P[2]=%f\ printf(\" ,\");

printf(\"有功功率增量P[3]=%f\

printf(\"无功功率增量Q[1]=%f\ printf(\" ,\");

printf(\"无功功率增量Q[2]=%f\ printf(\" ,\");

printf(\"电压增量V32=%f\

printf(\"\\n\");

//************筛选出最大值***********************

max=fabs(P[1])>fabs(P[2])?fabs(P[1]):fabs(P[2]); max=max>fabs(P[3])?max:fabs(P[3]); max=max>fabs(Q[1])?max:fabs(Q[1]); max=max>fabs(Q[2])?max:fabs(Q[2]); max=max>fabs(V32)?max:fabs(V32); printf(\"max=%f\\n\

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

//******************************************** while (max>0.00001) {

//****形成雅克比矩阵********************** for(j=1;j<4;j++) {if(1==j)

{float c=0,d=0; int m;

for(m=1;m<5;m++)

{c+=G[1][m]*e[m]-B[1][m]*f[m]; d+=G[1][m]*f[m]+B[1][m]*e[m]; }

J[1*N-1][j*N-1]=-c-G[1][j]*e[1]-B[1][j]*f[1]; J[1*N-1][j*N]=-d+B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N-1]=d+B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N]=-c+G[1][j]*e[1]+B[1][j]*f[1]; } else

{J[1*N-1][j*N-1]=-G[1][j]*e[1]-B[1][j]*f[1]; J[1*N][j*N]=G[1][j]*e[1]-B[1][j]*f[1]; J[1*N-1][j*N]=B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N-1]=B[1][j]*e[1]-G[1][j]*f[1]; } }

for(j=1;j<4;j++) {if(2==j)

{float c=0,d=0; int m;

for(m=1;m<5;m++)

{c+=G[2][m]*e[m]-B[2][m]*f[m]; d+=G[2][m]*f[m]+B[2][m]*e[m]; }

J[2*N-1][j*N-1]=-c-G[2][j]*e[2]-B[2][j]*f[2]; J[2*N-1][j*N]=-d+B[2][j]*e[2]-G[2][j]*f[2]; J[2*N][j*N-1]=d+B[2][j]*e[2]-G[2][j]*f[2]; J[2*N][j*N]=-c+G[2][j]*e[2]+B[2][j]*f[2]; } else

{J[2*N-1][j*N-1]=-G[2][j]*e[2]-B[2][j]*f[2]; J[2*N][j*N]=G[2][j]*e[2]-B[2][j]*f[2]; J[2*N-1][j*N]=B[2][j]*e[2]-G[2][j]*f[2]; J[2*N][j*N-1]=B[2][j]*e[2]-G[2][j]*f[2]; } }

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

for(j=1;j<4;j++) {if(3==j)

{float c=0,d=0; int m;

for(m=1;m<5;m++)

{c+=G[3][m]*e[m]-B[3][m]*f[m]; d+=G[3][m]*f[m]+B[3][m]*e[m]; }

J[3*N-1][j*N-1]=-c-G[3][j]*e[3]-B[3][j]*f[3]; J[3*N-1][j*N]=-d+B[3][j]*e[3]-G[3][j]*f[3]; J[3*N][j*N-1]=-2*e[3]; J[3*N][j*N]=-2*f[3]; } else

{J[3*N-1][j*N-1]=-G[3][j]*e[3]-B[3][j]*f[3]; J[3*N-1][j*N]=B[3][j]*e[3]-G[3][j]*f[3]; J[3*N][j*N-1]=0; J[3*N][j*N]=0; } }

printf(\"雅克比矩阵是:\\n\"); for(i=1;i<7;i++) {for(j=1;j<7;j++) {printf(\"%f\printf(\" \"); }

printf(\"\\n\"); }

//********计算修正方程************* for(i=1;ifor(i=1;iU[1][i]=J[1][i];

L[i][1]=J[i][1]/U[1][1]; }

for(n=2;nfor(j=n;jsigma1=0;

for(s=0;s<=n-1;s++)

sigma1+=L[n][s]*U[s][j]; U[n][j]=J[n][j]-sigma1;

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

}

for(i=n;isigma2=0;

for(s=0;s<=n-1;s++)

sigma2+=L[i][s]*U[s][n];

L[i][n]=(J[i][n]-sigma2)/U[n][n]; } }

b[1]=P[1];b[2]=Q[1];b[3]=P[2]; b[4]=Q[2];b[5]=P[3];b[6]=V32;

for(i=1;isigma1=0;

for(n=1;n<=i-1;n++) sigma1+=L[i][n]*y[n]; y[i]=b[i]-sigma1; }

for(i=M-1;i>=1;i--) {

sigma2=0;

for(n=i+1;nxe[1]=-x[1];xe[2]=-x[3];xe[3]=-x[5]; xf[1]=-x[2];xf[2]=-x[4];xf[3]=-x[6]; printf(\"节点电压:\\n\"); for(i=1;i<4;i++) {e[i]+=xe[i]; f[i]+=xf[i]; }

for(i=1;i<4;i++) {printf(\"e[%d]=\ printf(\"%f\ printf(\" ,\"); }

for(i=1;i<4;i++) {printf(\"f[%d]=\ printf(\"%f\ printf(\" ,\"); }

printf(\"\\n\"); k1=k1+1;

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

goto loop1; }

for(j=1;j<5;j++)

{a1+=G[4][j]*e[j]-B[4][j]*f[j]; b1+=G[4][j]*f[j]+B[4][j]*e[j]; }

P4=e[4]*a1+f[4]*b1; Q4=f[4]*a1-e[4]*b1;

printf(\"P4+Q4=%f+j%f\printf(\"\\n\");

p12=-2*e[1]*f[1]*y120-(e[1]*(e[1]-e[2])-f[1]*(f[1]-f[2]))*G[1][2]+(e[1]*(f[1]-f[2])+f[1]*(e[1]-e[2]))*B[1][2];

q12=-(e[1]*e[1]-f[1]*f[1])*y120+(e[1]*(e[1]-e[2])+f[1]*(f[1]-f[2]))*B[1][2]+(e[1]*(f[1]-f[2])+f[1]*(e[1]-e[2]))*G[1][2];

p13=(e[1]*e[1]-f[1]*f[1])*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+2*e[1]*f[1]*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-(e[1]*(e[1]-e[3])-f[1]*(f[1]-f[3]))*G[1][3]+(e[1]*(f[1]-f[3])+f[1]*(e[1]-e[3]))*B[1][3]; q13=(e[1]*e[1]-f[1]*f[1])*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-2*e[1]*f[1]*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+(e[1]*(e[1]-e[3])-f[1]*(f[1]-f[3]))*B[1][3]-(e[1]*(f[1]-f[3])+f[1]*(e[1]-e[3]))*G[1][3]; p14=-2*e[1]*f[1]*y140-(e[1]*(e[1]-e[4])-f[1]*(f[1]-f[4]))*G[1][4]+(e[1]*(f[1]-f[4])+f[1]*(e[1]-e[4]))*B[1][4];

q14=-(e[1]*e[1]-f[1]*f[1])*y140+(e[1]*(e[1]-e[4])-f[1]*(f[1]-f[4]))*B[1][4]+(e[1]*(f[1]-f[4])+f[1]*(e[1]-e[4]))*G[1][4];

p21=-2*e[2]*f[2]*y120-(e[2]*(e[2]-e[1])-f[2]*(f[2]-f[1]))*G[2][1]+(e[2]*(f[2]-f[1])+f[2]*(e[2]-e[1]))*B[2][1];

q21=-(e[2]*e[2]-f[2]*f[2])*y120+(e[2]*(e[2]-e[1])+f[2]*(f[2]-f[1]))*B[2][1]+(e[2]*(f[2]-f[1])+f[2]*(e[2]-e[1]))*G[2][1];

p24=-2*e[2]*f[2]*y240-(e[2]*(e[2]-e[4])-f[2]*(f[2]-f[4]))*G[2][4]+(e[2]*(f[2]-f[4])+f[2]*(e[2]-e[4]))*B[2][4];

q24=-(e[2]*e[2]-f[2]*f[2])*y240+(e[2]*(e[2]-e[4])+f[2]*(f[2]-f[4]))*B[2][4]+(e[2]*(f[2]-f[4])+f[2]*(e[2]-e[4]))*G[2][4];

p31=(e[3]*e[3]-f[3]*f[3])*(1-k)*z13r/(z13r*z13r+z13m*z13m)+2*e[3]*f[3]*(1-k)*z13m/(z13r*z13r+z13m*z13m)-(e[3]*(e[3]-e[1])-f[3]*(f[3]-f[1]))*G[3][1]+(e[3]*(f[3]-f[1])+f[3]*(e[3]-e[1]))*B[3][1];

q31=-(e[3]*e[3]-f[3]*f[3])*(1-k)*z13m/(z13r*z13r+z13m*z13m)+2*e[3]*f[3]*(1-k)*z13r/(z13r*z13r+z13m*z13m)-(e[3]*(e[3]-e[1])-f[3]*(f[3]-f[1]))*B[3][1]+(e[3]*(f[3]-f[1])+f[3]*(e[3]-e[1]))*G[3][1];

p41=-2*e[4]*f[4]*y140-(e[4]*(e[4]-e[1])-f[4]*(f[4]-f[1]))*G[4][1]+(e[4]*(f[4]-f[1])+f[4]*(e[4]-e[1]))*B[4][1];

q41=-(e[4]*e[4]-f[4]*f[4])*y140+(e[4]*(e[4]-e[1])-f[4]*(f[4]-f[1]))*B[4][1]+(e[4]*(f[4]-f[1])+f[4]*(e[4]-e[1]))*G[4][1];

p42=-2*e[4]*f[4]*y240-(e[4]*(e[4]-e[2])-f[4]*(f[4]-f[2]))*G[4][2]+(e[4]*(f[4]-f[2])+f[4]*(e[4]-e[2]))*B[4][2];

q42=-(e[4]*e[4]-f[4]*f[4])*y240+(e[4]*(e[4]-e[2])+f[4]*(f[4]-f[2]))*B[4][2]+(e[4]*(f[4]-f[2])+f[4]*(e[4]-e[2]))*G[4][2];

printf(\"s12=%f+j%f\printf(\"s13=%f+j%f\

--完整版学习资料分享----

-----WORD格式--可编辑--专业资料-----

printf(\"s14=%f+j%f\printf(\"s21=%f+j%f\printf(\"s24=%f+j%f\printf(\"s31=%f+j%f\printf(\"s41=%f+j%f\printf(\"s42=%f+j%f\}

本科生课程设计成绩评定表

姓 名 专业、班级 性 别 课程设计题目：控制系统的超前校正 课程设计答辩或质疑记录： --完整版学习资料分享----

-----WORD格式--可编辑--专业资料----- 成绩评定依据： 评 定 项 目 1．选题合理、目的明确（10分） 2．设计方案正确、具有可行性、创新性（20分） 3．设计结果（20分） 4．态度认真、学习刻苦、遵守纪律（15分） 5．设计报告的规范化、参考文献充分（不少于5篇）（10分） 6．答辩（25分） 总 分 评分成绩 最终评定成绩（以优、良、中、及格、不及格评定） 指导教师签字：

年 月 日

--完整版学习资料分享----