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

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个独立节点的电流平衡方程如下：

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y34Y24Y42y24；；

Y44y40y24y34Y34Y43y34。

Y12Y21y12Y23Y32y23一般的，对于有n个独立节点的网络，可以列写n个节点方程

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

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

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

或缩写为

YUI （2-5）

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

Yij称为节点i、j 间的互导纳，

它等于直接接于节点i、j间的支路导纳的负值。若节点i、j间不存在直接支路，则有

Yij0。由此可知节点导纳矩阵是一个稀疏的对称矩阵。

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

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

牛顿—拉夫逊法（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

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

'于是，3-9便简化为

f(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)

)

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

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

x迭代过程的收敛判据为

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

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轴的交点便确定了下一个近似值yxxx(k1)。由此可见，牛顿－拉夫逊法实质上就是切线法，是一种逐步线性化的方

法。

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

x(0)1，

x(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1122nnn

(3-10)

将上式中的n个多元函数在初始值附近分别展成泰勒级数,并略去含有

(0)x1(0)，x2，……，xn二次及以上阶次的各项，便得

(0)



f(0)f(0)f|x|xxxff(0)(0)(0)(0)f(x,x,...,x)|x|x(0)xx1(x(0),x(0),...,x(0))12n11011022...11212n011022f(0)|xxf...x(0)|x10n10nn00n.

......fn(x(0),x(0),...,x(0))12nf(0)f(0)|x|xxx11011022f...x(0)|x10nn0 (3-11)

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

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) ffff|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 的线性方程组,称为牛顿法的修正方程式.利用高斯消去法或三角分解法可以解出修正量x1，

(0)x2(0)，……，xn。然后对初始近似值进行修正

(0)(0)

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

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

f1|kx1((k),x(k),...,x(k))1x12nf2(x(k),x(k),...,x(k))|k212n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也可以缩写为

FxJ(k)(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个元素

f是第n个函数Jxiijif(x,x,...,x,)对第j个变量xi12ni,,,j的偏导数；上

角标(k)表示J阵的每一个元素都在点过程一直到满足收敛判据

maxf(x(k)1x(k)2...,x(k)n)处取值。迭代

n1f(x(k),x(k),...,x(k))i12 (3-18)或

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

i为止。

1和

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

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

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

节点电压用直角坐标表示是的牛顿-拉夫逊法潮流计算

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

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

YijGijjBij

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

nn

Qifi(GijejBijfj)ei(GijfjBijej)j1j1nYUijjinnj的右端，展开并分出实部和虚

Piei(GijejBijfj)fi(GijfjBijej)（11-45）

j1j1

假定系统中的第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 222222ViVisViVis(eifi)0 （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）

nn已经具备了方程组（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

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f11mmm1m1n1n1PPPPPPPPmmmmmmmme1f1emfmem1fm1en1fn1QQQQQQQQmmmmmmmmef1emfmem1fm1en1fn11 JPPPPPPPPm1m1m1m1m1m1m1m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时

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）  Qnif(GikekBikfk)GiieiBiifiik1V2ie2eiiV2if2fii

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

W1J11J12J1.n1 W2J21J22J2.n1Wn1Jn1.1Jn1.2Jn1.n1（11-51）

式中，Wi和Vi都是二维列向量；Jij是22介方阵。

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

i对于PQ节点

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

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

i

V1V2Vn111-52） （

Piej JijV2iejPifj （11-53） 2Vifj从表达式（1-7）~（1-11）可以看到，雅克比矩阵有以下特点：

（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•••

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

输入原始数据 形成节点导纳矩阵 给定节点电压初值ei(0),fi(0) k=0 (k)计算P,Qi(k)及Vii2(k) 是 是否收敛 否 计算雅可比矩阵各元素 解修正方程式，求ei(k),fi(k) k=k+1 ei(k1)ei(k)ei(k)fi(k1)fi(k)fi(k) 计算平衡节点功率及全部线路功率 输出

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

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

计算节点导纳矩阵

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

图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]=; B[2][3]=B[3][2]=;

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]=; B[3][4]=B[4][3]=; for(i=1;i<5;i++) {for(j=1;j<5;j++)

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

printf(\"\\n\");

计算各节点不平衡量

按上一次迭代计算出的节点电压值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

同样的可以算出

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

P3(0)P3sP3(0)V2(0)3V3sV3s2(0)2

其部分C语言程序如下：

printf(\"迭代次数k1=%d\\n\ for (i=1;i<5;i++)

PiQi(GijeiBijfi)ejfjPiQiBijeiGijfi)fjejV2iV2i0ejfj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 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);maxPi、Qi、Vi(k)(k)2(k)105ijjiWJV..2.22WP1Q1...PmQmPm1Vm1...Pn1Vn1....TSijPijjQijViIijViyi0Vi(ViVj)yij力系统分析上册[M]．武汉：华中科技

大学出版社.

[2] 何仰赞等.电力系统分析下册[M]．武汉：华中科技大学出版社.

[3] 陈珩.电力系统稳态分析. 水利电力出版社,1994

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

主程序如下：

#include <>

附录

#include <>

#include <> #define N 2 #define M 7 main() {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,,},f[5]={0},P[5],Q[5],Ps[5]={0,,,},xe[4],xf[4]; float Qs[5]={0,,},V3s=,V4S=; 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]=; B[2][3]=B[3][2]=;

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]=; B[3][4]=B[4][3]=; 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\

//******************************************** while (max> {

//****形成雅克比矩阵********************** 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]; } }

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;i{L[i][i]=1; }

for(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; }

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;nx[i]=(y[i]-sigma2)/U[i][i]; }

xe[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; 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\ 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\ }

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

姓 名 专业、班级 性 别 课程设计题目：控制系统的超前校正 课程设计答辩或质疑记录： 成绩评定依据： 评 定 项 目 1．选题合理、目的明确（10分） 2．设计方案正确、具有可行性、创新性（20分） 3．设计结果（20分） 4．态度认真、学习刻苦、遵守纪律（15分） 5．设计报告的规范化、参考文献充分（不少于5篇）（10分） 6．答辩（25分） 总 分 评分成绩 最终评定成绩（以优、良、中、及格、不及格评定）

指导教师签字：

年 月 日