直角坐标系下牛顿法潮流计算
1电力系统潮流计算
潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角)、网络中的功率分布及功率损耗等。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性.可靠性和经济性。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。
2节点导纳矩阵的形成
在图1（a ）的简单电力系统中，若略去变压器的励磁功率和线路电容，负荷用阻抗表示，便可以得到一个有5个节点（包括零电位点）和7条支路的等值网络，如图1（b ）所示。将接于节点1和4的电势源和阻抗的串联组合变换成等值的电流源和导纳的并联组合，变得到图1（c ）的等值网络，其中1101I y E =和
4404I y E =分别称为节点1和4的注入电流源。
(a)
2
4?
?4
y (c)
图1 电力系统及其网络
以零电位点作为计算节点电压的参考点，根据基尔霍夫定律，可以写出4个独立节点的电流平衡方程如下：
10112121
12212022323242423323434244234434044()()()()0()()0()()y U y U U I y U U y U y U U y U U y U U y U U y U U y U U y U I ?+-=?
-++-+-=?
?
-+-=?
?-+-+=? （2-1） 上述方程组经过整理可以写成
1111221211222233244322333344422433444400Y U Y U I Y U Y U Y U Y U Y U Y U Y U Y U Y U Y U I ?
+ =?
+++=?
?
++=?? ++=? （2-2）
式中，
111012
Y y y =+；
2220232412
Y y y y y =+++；332334Y y y =+；44402434Y y y y =++；
122112Y Y y ==-；
233223
Y Y y ==-；
244224
Y Y y ==-；
344334
Y Y y ==-。
一般的，对于有n 个独立节点的网络，可以列写n 个节点方程
11112211211222221122n n n n n n nn n n Y U Y U Y U I Y U Y U Y U I Y U Y U Y U I ?
+++=?
+++=?
? ?
?++
+=?
（2-3）
也可以用矩阵写成
1111121212222212n n n n nn n n U I Y Y Y Y Y Y U I Y Y Y U I ????
???????? ??????=?????? ?????? ?????????? （2-4）
或缩写为
YU I = （2-5）
矩阵Y 称为节点导纳矩阵。它的对角线元素ii Y 称为节点i 的自导纳，其值等于接于节点i 的所有支路导纳之和。非对角线元素
ij
Y 称为节点i 、j 间的互导纳，
它等于直接接于节点i 、j 间的支路导纳的负值。若节点i 、j 间不存在直接支路，则有
ij Y =。由此可知节点导纳矩阵是一个稀疏的对称矩阵。
3牛顿-拉夫逊法潮流计算
牛顿-拉夫逊法的基本原理
牛顿—拉夫逊法（Newton —Raphson 法）是求解非线性方程代数方程组的有效迭代计算方法。在牛顿—拉夫逊法的每一次迭代过程中，对非线性方程通过线性化处理逐步近似。下面以单变量加以说明。
设有单变量非线性方程
()0f x = (3-1)
求解此方程时。先给出解的近似值
(0)
x
它与真解的误差为
(0)
x
?，则
(0)(0)
x x x
=+?将满足方程，即
(0)(0)
()0f x x +?= (3-2)
将(3-8)式左边的函数在(0)
x
附近展成泰勒级数，于是便得
2
'
''
()
(0)
(0)(0)
(0)
(0)
(0)
(0)
(0)
(0)
()()()()
......()
....
2!
!
()
()
n
n f f n x x f f f
x
x x x
x x
x +?=+
?+
++
+?? （3-3）
式中'
(0)
()f
x ,……
()
(0)()n f
x 分别为函数()f x 在(0)
x 处的一阶导数，….，n
阶导数。
如果差值(0)
x ?很小，3-9式右端(0)
x
?的二次及以上阶次的各项均可略去。
于是，3-9便简化为
'
(0)
(0)(0)
(0)
(0)
(
)()()f f f x
x x x
x
+?=+
?＝0 (3-4)
这是对于变量的修正量(0)
x
?的现行方程式，亦称修正方程式。解此方程可得
修正量
(0)
(0)
'
(0)
()
()
f x x
f x
?=-
(3-5)
用所求的(0)
x
?去修正近似解，变得
(0)
(1)
(0)
(0)
(0)
'
(0)
()
()
f x x
x x
x f x
=+?=-
(3-6)
由于3-10是略去高次项的简化式，因此所解出的修正量(0)
x
?也只是近似
值。修正后的近似解(1)
x 同真解仍然有误差。但是，这样的迭代计算可以反复进行下去，迭代计算的通式是
()
(1)
()
'
()
()
()
k k k k f x x
x
f
x +=-
(3-7)
迭代过程的收敛判据为
()
1
()k f x ε< (3-8)
或
()
2
k x
ε?< (3-9)
式中
1
ε，2
ε
为预先给定的小正数。
这种解法的几何意义可以从图3－1得到说明。函数y ＝f(x)为图中的曲线。
f(x)＝0的解相当于曲线与x 轴的交点。如果第k 次迭代中得到()
k x
，则过
()()(),()k k k f y x x ??=????
点作一切线，此切线同x 轴的交点便确定了下一个近似值(1)
k x
+。由此可见，牛顿－拉夫逊法实质上就是切线法，是一种逐步线性化的方
法。
应用牛顿法求解多变量非线性方程组3-1时，假定已给出各变量的初值
1
(0)
x ，
2
(0)
x ….
(0)
n
x ，令
1
(0)
x ?，
2
(0)
x ?，…..
(0)
n
x ?分别为各变量的修
正量，使其满足方程3-2即
11
12221
1221122(0)(0)(0)(0)(0)(0)(,,....,)0
(0)(0)(0)(0)(0)(0)(,,....,)0......(0)(0)(0)(0)(0)(0)(,,....,)0n n n n n
n n f x x x x x x f x x x x x x f x x x x x x ?+?+?+?=??+?+?+?=??
?
?
+?+?+?=??
(3-10)
将上式中的n 个多元函数在初始值附近分别展成泰勒级数,并略去含有
)0(1x ?，)
0(2
x ?，……，)
0(n x ?二次及以上阶次的各项，便得
11
10
01
1
2
12
1
2
11
10
021212
1
2
1
01
2
1
1
(0)(0)(0)(0)(0)
(0)
(,,...,)...0(0)(0)(0)(0)(0)(0)(,,...,) 0
......
(0)(0)(0)(0)(,,...,)|||||||n
n
n n n
n
n
n
f f f f x x x x x x x x
x f f f f x x x x x x x x
x
f f
x x x x x ?
??
+
?+?++
?=??????+?+?++?=????
?+
?+?1
1002
2
(0)
(0) 0
||n
n
f f x x x
x
???????????
??++?=????
.
(3-11)
方程式3-17也可以写成矩阵形式
11100
0121122
2221200
1
2
1200
01
2
...(0)(0)(0)(,,...,)(0)(0)(0)(,,...,)..................
...
(0)(0)(0)(,,...,)...|||||||||n
n n n
n n n n n n f f f x x x f x x x f f f f x x x x x
x f
x x x f f f x x
x ???
?
??????
????
?????????=-????????
??
?
?
??
?????1
2
(0)(0)...(0)n x x x ???????????
?
?????
??
????
?????????
?????
(3-12)
方程式3-18是对于修正量)
0(1x ?，)
0(2x ?，……，)
0(n x ? 的线性方程组,称为牛顿法的修正方程式.利用高斯消去法或三角分解法可以解出修正量)
0(1x ?，
)0(2x ?，……，)
0(n x ?。然后对初始近似值进行修正
(1)(0)(0)i
i
i
x x x =+? (i=1,2,….,n) (3-13)
如此反复迭代，在进行k ＋1次迭代时，从求解修正方程式
111121122
222121
2
121
2
...()()()(,,...,)()()()(,,...,)..................
...
()()()(,,...,)...|||||||||k
k
k n
n n k k
k
n
n n n n n k k
k n k k k k k k k k k f f f x x x f x x x f f f f x x x x x
x f
x x x f f f x x
x ???
?
??????
????
?????????=-????????
??
?
?
??
?????1
2
()()...()n k k k x x x ?
??????????
?
?????
??
????
?????????
?????
(3-14)
得到修正量1
()k x ?，2
()k x ?，()n
k x ?，并对各变量进行修正
(1)()()
i
i
i
k k k x x x +=+? (i=1,2,…,n) (3-15)
式3-20和3-21也可以缩写为
())()
()
(k k k x J
x F
?-= (3-16)
和 )()()
1(k k k x x x
?+=+ (3-17)
式中的X 和X ?分别是由n 个变量和修正量组成的n 维列向量；F(X)是由n 个多元函数组成的n 维列项量；J 是n 阶方阵，称为雅可比矩阵，它的第i 、j 个元素
i
ij
i
f J
x
=??是第n 个函数1
2
(,,...,,)n
i
f x x x 对第j 个变量j
x
的偏导数；上
角标(k)表示J 阵的每一个元素都在点,
,
,
()()()(...,)12i
k k k n f x x x 处取值。迭代
过程一直到满足收敛判据
{}1
1
2
()()()max
(,,...,)i
n
k k k f x x x ε
< (3-18)或
{}
2()max i
k x ε?< (3-19)
为止。
1
ε
和
2
ε
为预先给定的小正数。
将牛顿－拉夫逊法用于潮流计算，要求将潮流方程写成形如方程式3-1的
形式。由于节点电压可以采用不同的坐标系表示，牛顿－拉夫逊法潮流计算也将相应的采用不同的计算公式。
节点电压用直角坐标表示是的牛顿-拉夫逊法潮流计算
采用直角坐标时，节点电压可表示为
i
i i jf e V += 导纳矩阵元素则表示为
ij ij ij jB G Y +=
将上述表示式代入n
i j
i i i i i ij j i
S P jQ U I U Y U
*
**
==+==∑的右端，展开并分出实部和虚
部，便得
1
1
()()n
n
i i ij j ij j i ij j ij j j j P e G e B f f G f B e ===-++∑∑
（11-45）
假定系统中的第1，2，3···，m 号节点为PQ 节点，第i 个节点的给定功率设为is P 和is Q ，对对该节点可列写方程
0)()(0)()(1
1
1
1=+---=-=?=+---=-=?∑∑∑∑====n
j j ij j ij i n j j ij j ij i is i is i n
j j ij j ij i n j j ij j ij i is i is i e B f G f f B e G e P P P P e B f G f f B e G e P P P P
（i=1,2，···，m ） （11-46）
假定系统中的第m+1，m+2，···，n-1号节点为PV 节点，则对其中每一个节点可以列写方程
??
?
?
?=+-=-=?=+---=-=?∑∑==0)(0)()(2
2222211i i is i is i n j n
j j ij j ij i j ij j ij i is i is i f e V V V V e B f G f f B e G e P P P P （i=m+1，m+2，···，n-1） （11-47）
第n 号节点为平衡点，其电压n n n jf e V +=是给定的，故不参加迭代。 式（11-46）和式（11-47）总共包含了2（n-1）个方程，待求的变量有
1111,,...,,--n n f e f e 也是2（n-1）个。我们还可看到，方程（11-46）和式（11-47）
已经具备了方程组（11-34）的形式。因此，不难写出如下的修正方程式
V J W ?-=? （11-48）
??
???
11
()()n n i i ij j ij j i ij j ij j j j Q f G e B f e G f B e ===--+∑∑??
?
??
式中
[]
T
n n m m m m V P V P Q P Q P W 2
1121111......--++????????=?
[]T n n m m m m f e f e f e f e V 111111......--++????????=?
?????????
?
?
???????????
??
???
???????????????
??????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????=----+-+---------+-+------+-+++++++++-+-+++++++++--++--++--++--++112
112112112121211
211211111111111111112
112112112121211211211111111111111111111
111111111111111111111111111111111
1
n n n n m n m n m n m n n n n n n n m n m n m n m n n n n m n m m m m m m m m m m m n m n m m m m m m m m m m m n m n m m m m m m m m m m m n m n m m m m m m m m m m m n n m m m m n n m m m m f V e V f V e V f V e V f V e V f P e P f P e P f P e P f P e P f V e V f V e V f V e V f V e V f P e P f P e P f P e P f P e P f Q e Q f Q e Q f Q e Q f Q e Q f P e P f P e P f P e P f P e P f Q e Q f Q e Q f Q e Q f Q e Q f P e P f P e P f P e P f P e P J 上述方程中雅克比矩阵的各元素，可以对式（11-46）和式（11-47）求偏导数获得。当j i ≠时
?????
?
??
???=???=???-=???=???+-=???-=???0
))(22j i
j i i ij i ij j
i j i i ij i ij j i
j i f V e V f G e B e Q f P f B e G f Q e P （11-49） 当i j =时
???
?????
?
?
???
???
???
-=???-=???++--=???-++=???-++-=???----=???∑∑∑∑====i
i i
i
i
i n
k i ii i ii k ik k ik i i
n k i ii i ii k ik k ik i i
n k i ii i ii k ik k ik i i
n
k i ii i ii k ik k ik i i f f V e e V f B e G f B e G f Q f G e B e B f G e Q f G e B e B f G f P f B e G f B e G e P 22)()()()(2
2
1
111
（11-50）
修正方程式（11-48）还可以写成分块矩阵的形式
????
?
??????????????????????-=?
?????
?????????--------1211.12.11
.11.222211.11211121n n n n n n n n V V V J J J J J J J J J W W W （11-51）
式中，i W ?和i V ?都是二维列向量；ij J 是22?介方阵。
??
??????=?i i i f e V
对于PQ 节点
?
?
?
?????=?i i i Q P W ??
??
?
????
???????????????=j i j i j i j
i ij f Q e Q f P e
P J （11-52） 对于PV 节点
???
???????=?2i i i V P W
?????
?
?
????
???????????????=j i j i j i j
i ij f V e V f P e
P J 22 （11-53） 从表达式（1-7）~（1-11）可以看到，雅克比矩阵有以下特点：
（1）雅克比矩阵各元素都是节点电压的函数，它们的数值将在迭代过程中不断的改变。
（2）雅克比矩阵的子块ij J 中的元素的表达式只用到导纳矩阵中的对应元素ij Y 。若0=ij Y ，则必有0=ij J 。因此，式（1-9）式中分块形式的雅克比矩阵同节点导纳矩阵一样稀疏，修正方程的求解同样可以用稀疏矩阵的求解技巧。
（3）无论在式（1-6）或式（1-9）中雅克比矩阵的元素或子块都不具有对称性。
用牛顿-拉夫逊法计算潮流的流程框图示于图11-27。首先要输入网络的原始数据以及各节点的给定值并形成节点导纳矩阵。输入节点电压初值)0(i e 和)0(i f ，置迭代计数k=0。然后开始进入牛顿法的迭代过程。在进行第k+1次迭代时，其计算步骤如下：
（1）按上一次迭代计算出的节点电压值)(k e 和)(k f ，利用式（11-46）和式（11-47）计算各类节点的不平衡量)(k i P ?、)(k i Q ?和)(2k i V ?。
（2）按条件（11-43）校验收敛，即
{}
)(2)((k )i P max k i k i V Q ???、、<ε
如果收敛，迭代到此结束，转入计算各线路潮流和平衡节点的功率，并打印输出计算结果。不收敛则继续计算。
（3）利用式（11-49）和式（11-50）计算雅克比矩阵的各元素。 （4）解修正方程式（11-48）求节点电压的修正量)(k i e ?和)(k i f ?。 （5）修正各节点的电压
)()()1()()()1(,
k i k i k i k i k i k i f f f e e e ?+=?+=++
（6）迭代计数加1，返回第一步继续迭代过程。
迭代结束后，还要计算出平衡节点的功率和网络中的功率分布。 输电线路功率的计算公式如下：
ij j i i i i ij i ij ij ij y V V V y V I V jQ P S *
?
?
?
*
*
?
-+==+=)(02
牛顿-拉夫逊法潮流计算程序框图
图2 牛顿-拉夫逊法潮流计算程序框图
4 C语言实现直角坐标表示的牛顿拉夫逊的计算
计算节点导纳矩阵
该电力系统的等值电路如图所示
图3电力系统等值电路
根据节点方程和等值电路图，计算节点导纳矩阵的C语言子程序如下。
G[1][1]=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14 r*z14r+z14m*z14m);
B[1][1]=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z1 4r* 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",G[i][j],B[i][j]); printf(" "); }
printf("\n");
计算各节点不平衡量
按上一次迭代计算出的节点电压值)(k e 和)(k f （当k=0时即为给定值），利用式（1-2和式（1-3）计算各类节点的不平衡量)(k i P ?、)(k i Q ?和)(2k i V ?。
同样的可以算出
?
?????----=-=???
????-+--=-=?∑∑∑∑====4141)
0(1)0(1)0(1)0(1)0(1)0(11)0(11)0(14141)
0(1)0(1)0(1)0(1)0(1)0(11)0(1
1)
0(1
)()()()(j j j j j j j j j j s s j j j j j j j j j j s s e B f G e f B e G f Q Q Q Q e B f G f f B e G e P P P P 2
)0(32
3)
0(23
)0(33)0(3)0(2
2)0(2)
0(22)0(2s
s s s s V V V
P P P Q Q Q P P P -=?-=?-=?-=?
其部分C 语言程序如下：
printf("迭代次数k1=%d\n",k1); 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);{
})
(2)((k )i P max k i
k i V Q ???、、ε5
10
-=εj i ≠i
j =V
J W ?-=?[]
T
n n m m m m V P V P Q P Q P W 2
1121111......--++????????=?ij j i i i i ij i ij ij ij y V V V y V I V jQ P S .
.
.
.
0.
2.
.
)(-+==+=力系统分析上册[M]．武汉：华中科技
大学出版社.
[2] 何仰赞等.电力系统分析下册[M]．武汉：华中科技大学出版社.
?????
?
??
???=???=???-=???=???+-=???-=???0
))(22j i
j i i ij i ij j
i j i i ij i ij j i
j i f V e V f G e B e Q f P f B e G f Q e P ???
?????
?
?
???
???
???
-=???-=???++--=???-++=???-++-=???----=???∑∑∑∑====i
i i
i
i
i n
k i ii i ii k ik k ik i i
n k i ii i ii k ik k ik i i
n k i ii i ii k ik k ik i i
n
k i ii i ii k ik k ik i i f f V e e V f B e G f B e G f Q f G e B e B f G e Q f G e B e B f G f P f B e G f B e G e P 22)()()()(2
2
1
111
[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",&z12r,&z12m);
printf("请输入z13的实部和虚部\n");
scanf("%f%f",&z13r,&z13m);
printf("请输入z14的实部和虚部\n");
scanf("%f%f",&z14r,&z14m);
printf("请输入z24的实部和虚部\n");
scanf("%f%f",&z24r,&z24m);
printf("请输入y120的值\n");
scanf("%f",&y120);
printf("请输入y140的值\n");
scanf("%f",&y140);
printf("请输入y240的值\n");
scanf("%f",&y240);
printf("请输入变比k的值\n");
scanf("%f",&k);
//****************计算导纳矩阵*******************
G[1][1]=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z1 4r+z14m*z14m);
B[1][1]=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z 14r+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",G[i][j],B[i][j]);
printf(" ");
}
printf("\n");//形成节点导纳矩阵
//******************************************* }
printf("\n");
//******************************************** //计算各节点不平衡量
loop1:
printf("迭代次数k1=%d\n",k1);
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",P[1]);
printf(" ,");
printf("有功功率增量P[2]=%f",P[2]);
printf(" ,");
printf("有功功率增量P[3]=%f",P[3]);
printf("无功功率增量Q[1]=%f",Q[1]);
printf(" ,");
printf("无功功率增量Q[2]=%f",Q[2]);
printf(" ,");
printf("电压增量V32=%f",V32);
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",max);
//******************************************** 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",J[i][j]);
printf(" ");
}
printf("\n");
}
//********计算修正方程*************
for(i=1;i {L[i][i]=1; } for(i=1;i { U[1][i]=J[1][i]; L[i][1]=J[i][1]/U[1][1]; } for(n=2;n { for(j=n;j { sigma1=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;i { sigma2=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;i { sigma1=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;n sigma2+=U[i][n]*x[n]; x[i]=(y[i]-sigma2)/U[i][i]; } xe[1]=-x[1];xe[2]=-x[3];xe[3]=-x[5]; 前言 电力系统潮流计算是研究电力系统稳态运行情况的一种计算,它根据给定的运行件及系统接线情况确定整个电力系统各部分的运行状态。在电力系统规划设计和现有电力系统运行方式的研究中,都需要利用潮流计算来定量分析、比较供电方案或运行方式的合理性、可靠性和经济性。本次课程设计任务是闭式网络的潮流计算，用到的方法为牛顿拉夫逊极坐标法潮流计算。 牛顿法是数学中解决非线性方程式的典型方法，有较好的收敛性。解决电力系统潮流计算问题是以导纳距阵为基础的，因此，只要在迭代过程中尽可能保持方程式系数距阵的稀疏性，就可以大大提高牛顿法潮流程序的放率。自从20 世纪60 年代中期利用了最佳顺序消去法以后，牛顿法在收敛性、内存要求、速度方面都超过了阻抗法，成为直到目前仍在广泛采用的优秀方法。 目录 1任务书 (2) 2.模型简介及等值电路 (3) 3.设计原理 (4) 4.修正方程的建立 (7) 5.程序流程图及MATLAB程序编写 (9) 6.结果分析 (15) 7.设计总结 (25) 8.参考文献 (26) 《电力系统分析》 课程设计任务书 题目极坐标表示的牛顿拉夫逊法潮流计算程序设计学生姓名学号专业班级 设计内容与要求1. 设计要求 掌握MATLAB语言编程方法；理解和掌握运用计算机进行潮流计算的基本算法原理；针对某一具体电网，进行潮流计算程序设计。 其目的在于加深学生对电力系统稳态分析中课程中基本概念和计算方法的理解，培养学生运用所学知识分析和解决问题的能力。 2. 内容 1）学习并掌握MATLAB语言。 2）掌握变压器非标准变比概念及非标准变比变压器的等值电路。掌握节点导纳矩阵的概念及导纳矩阵的形成和修改方法。 3）掌握电力系统功率方程、变量和节点分类。 4）掌握利用极坐标表示的牛-拉法进行潮流计算的方法和步骤。 5）选择一个某一具体电网，编制程序流程框图。 6）利用MATLAB语言编写该模型的潮流计算程序，并上机调试程序，对计算结果进行分析。 7）整理课程设计论文。 起止时间2013 年7 月 4 日至2013 年7月10日指导教师签名年月日 系（教研室）主任签 名 年月日学生签名年月日 2 模型简介及等值电路 2.1课程设计模型：模型3 摘要 本文，首先简单介绍了基于在MALAB中行潮流计算的原理、意义，然后用具体的实例，简单介绍了如何利用MALAB去进行电力系统中的潮流计算。 众所周知，电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各线的电压、各元件中流过的功率、系统的功率损耗等等。在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性和经济性。 此外，在进行电力系统静态及暂态稳定计算时，要利用潮流计算的结果作为其计算的基础；一些故障分析以及优化计算也需要有相应的潮流计算作配合；潮流计算往往成为上述计算程序的一个重要组成部分。以上这些，主要是在系统规划设计及运行方式安排中的应用，属于离线计算范畴。 牛顿－拉夫逊法在电力系统潮流计算的常用算法之一，它收敛性好，迭代次数少。本文介绍了电力系统潮流计算机辅助分析的基本知识及潮流计算牛顿－拉夫逊法，最后介绍了利用MTALAB程序运行的结果。 关键词：电力系统潮流计算，牛顿－拉夫逊法，MATLAB ABSTRACT This article first introduces the flow calculation based on the principle of MALAB Bank of China, meaning, and then use specific examples, a brief introduction, how to use MALAB to the flow calculation in power systems. As we all know, is the study of power flow calculation of power system steady-state operation of a calculation, which according to the given operating conditions and system wiring the entire power system to determine the operational status of each part: the bus voltage flowing through the components power, system power loss and so on. In power system planning power system design and operation mode of the current study, are required to quantitatively calculated using the trend analysis and comparison of the program or run mode power supply reasonable, reliability and economy. In addition, during the power system static and transient stability calculation, the results of calculation to take advantage of the trend as its basis of calculation; number of fault analysis and optimization also requires a corresponding flow calculation for cooperation; power flow calculation program often become the an important part. These, mainly in the way of system design and operation arrangements in the application areas are off-line calculation. Newton - Raphson power flow calculation in power system is one commonly used method, it is good convergence of the iteration number of small, introduce the trend of computer-aided power system analysis of the basic knowledge and power flow Newton - Raphson method, introduced by the last matlab run results. Keywords：power system flow calculation, Newton – Raphson method, matlab >> %本程序的功能是用牛顿拉夫逊法进行潮流计算n=input('请输入节点数:n='); nl=input('请输入支路数:nl='); isb=input('请输入平衡母线节点号:isb='); pr=input('请输入误差精度:pr='); B1=input('请输入由各支路参数形成的矩阵:B1='); B2=input('请输入各节点参数形成的矩阵:B2='); Y=zeros(n); e=zeros(1,n);f=zeros(1,n);V=zeros(1,n); O=zeros(1,n);S1=zeros(nl); for i=1:nl if B1(i,6)==0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end Y(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; Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; 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); end for i=1:n S(i)=B2(i,1)-B2(i,2); B(i,i)=B(i,i)+B2(i,5); end P=real(S);Q=imag(S); ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; while IT2~=0 IT2=0;a=a+1; for i=1:n if i~=isb C(i)=0; D(i)=0; for j1=1:n C(i)= C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1); D(i)= D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1); end 牛顿插值法 插值法是利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。当插值节点增减时全部插值基函数均要随之变化，这在实际计算中很不方便。为了克服这一缺点，提出了牛顿插值。牛顿插值通过求各阶差商，递推得到的一个公式： f(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0 )...(x-xn-1)+Rn(x)。 插值函数 插值函数的概念及相关性质[1] 定义：设连续函数y-f(x) 在区间[a,b]上有定义，已知在n+1个互异的点 x0,x1,…xn上取值分别为y0,y1,…yn （设a≤ x1≤x2……≤xn≤b)。若在函数类中存在以简单函数P(x) ，使得P(xi)=yi,则称P(x) 为f(x)的插值函数. 称x1,x2,…xn 为插值节点，称[a,b]为插值区间。 定理：n次代数插值问题的解存在且唯一。 牛顿插值法C程序 程序框图#include xxxx实验报告 学生姓名：学号：专业班级： 实验类型：□验证□综合■设计□创新实验日期：实验成绩： 一、实验项目名称 P-Q分解法潮流计算实验 二、实验目的与要求： 目的：电力系统分析的潮流计算是电力系统分析的一个重要的部分。通过对电力系统潮流分布的分析和计算，可进一步对系统运行的安全性，经济性进行分析、评估，提出改进措施。电力系统潮流的计算和分析是电力系统运行和规划工作的基础。 潮流计算是指对电力系统正常运行状况的分析和计算。通常需要已知系统参数和条件，给定一些初始条件，从而计算出系统运行的电压和功率等；潮流计算方法很多：高斯-塞德尔法、牛顿-拉夫逊法、P-Q分解法、直流潮流法，以及由高斯-塞德尔法、牛顿-拉夫逊法演变的各种潮流计算方法。 本实验采用P-Q分解法进行电力系统分析的潮流计算程序的编制与调试，获得电力系统中各节点电压，为进一步进行电力系统分析作准备。通过实验教学加深学生对电力系统潮流计算原理的理解和计算，初步学会运用计算机知识解决电力系统的问题，掌握潮流计算的过程及其特点。熟悉各种常用应用软件，熟悉硬件设备的使用方法，加强编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。 要求：编制调试电力系统潮流计算的计算机程序。程序要求根据已知的电力网的数学模型（节点导纳矩阵）及各节点参数，完成该电力系统的潮流计算，要求计算出节点电压、功率等参数。 三、主要仪器设备及耗材 每组计算机1台、相关计算软件1套 四、实验内容： 1.理论分析： P-Q分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思想是：把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，把有功功率和无功功率迭代分开来进行。 牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系统时， 课程设计说明书 题目电力系统分析系(部) 专业(班级) 学号 指导教师 起止日期 电力系统分析课程设计任务书系（部）：专业：指导教师： 目录 一、潮流计算基本原理 1.1潮流方程的基本模型 1.2潮流方程的讨论和节点类型的划分 1.3、潮流计算的意义 二、牛顿－拉夫逊法 2.1牛顿-拉夫逊法基本原理 2.2节点功率方程 2.3修正方程 2.4牛顿法潮流计算主要流程 三、收敛性分析 四、算例分析 总结 参考文献 电力系统分析潮流计算 一、潮流计算基本原理 1.1潮流方程的基本模型 电力系统是由发电机、变压器、输电线路及负荷等组成，其中发电机及负荷是非线性元件，但在进行潮流计算时，一般可以用接在相应节点上的一个电流注入量来代表。因此潮流计算所用的电力网络系由变压器、输电线路、电容器、电抗器等静止线性元件所构成，并用集中参数表示的串联或并联等值支路来模拟。结合电力系统的特点，对这样的线性网络进行分析，普通采用的是节点法，节点电压与节点电流之间的关系 V Y I = （1－1） 其展开式为 j n j ij i V Y I ∑==1 ),,3,2,1 (n i = （1－2） 在工程实际中，已经的节点注入量往往不是节点电流而是节点功率，为此必须应用联 系节点电流和节点功率的关系式 i i i i V jQ P I * -= ),,3,2,1(n i = （1－3） 将式（1－3）代入式（1－2）得到 j n j ij i i i V Y V jQ P ∑=* =-1 ),,3,2,1(n i = （1－4） 交流电力系统中的复数电压变量可以用两种极坐标来表示 i j i i e V V θ= （1－5） 或 i i i jf e V += （1－6） 而复数导纳为 %本程序的功能是用牛顿——拉夫逊法进行潮流计算 % B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳 % 5、支路的变比；6、支路首端处于K侧为1，1侧为0 % B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；% 3为PV节点； clear; n=10;%input('请输入节点数：n='); nl=10;%input('请输入支路数：nl='); isb=1;%input('请输入平衡母线节点号：isb='); pr=0.00001;%input('请输入误差精度：pr='); B1=[1 2 0.03512+0.08306i 0.13455i 1 0; 2 3 0.0068+0.18375i 0 1.02381 1; 1 4 0.05620+0.13289i 0.05382i 1 0; 4 5 0.00811+0.24549i 0 1.02381 1; 1 6 0.05620+0.13289i 0.05382i 1 0; 4 6 0.04215+0.09967i 0.04037i 1 0; 6 7 0.0068+0.18375i 0 1.02381 1; 6 8 0.02810+0.06645i 0.10764i 1 0; 8 10 0.00811+0.24549i 0 1 1; 8 9 0.03512+0.08306i 0.13455i 1 0] B2=[0 0 1.1 1.1 0 1; 0 0 1 0 0 2; 0 0.343+0.21256i 1 0 0 2; 0 0 1 0 0 2; 0 0.204+0.12638i 1 0 0 2; 0 0 1 0 0 2; 0 0.306+0.18962i 1 0 0 2; 0 0 1 0 0 2; 0.5 0 1.1 1.1 0 3; 0 0.343+0.21256i 1 0 0 2] ;%input('请输入各节点参数形成的矩阵：B2='); Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl); % % %--------------------------------------------------- for i=1:nl %支路数 if B1(i,6)==0 %左节点处于1侧 p=B1(i,1);q=B1(i,2); else %左节点处于K侧 p=B1(i,2);q=B1(i,1); end 潮流例题：根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择Matlab进行设计）。 2.在给定的电力网络上画出等值电路图。 3.运用计算机进行潮流计算。 4.编写设计说明书。 一、设计原理 1．牛顿-拉夫逊原理 牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不 平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。 牛顿—拉夫逊迭代法的一般步骤： （1）形成各节点导纳矩阵Y。 （2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。 （3）计算各个节点的功率不平衡量。 （4）根据收敛条件判断是否满足，若不满足则向下进行。 （5）计算雅可比矩阵中的各元素。 （6）修正方程式个节点电压 （7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。 （8）计算平衡节点输出功率和各线路功率 2．网络节点的优化 １）静态地按最少出线支路数编号 这种方法由称为静态优化法。在编号以前。首先统计电力网络个节点的出线支路数，然后，按出线支路数有少到多的节点顺序编号。当由n 个节点的出线支路相同时，则可以按任意次序对这n 个节点进行编号。这种编号方法的根据是导纳矩阵中，出线支路数最少的节点所对应的行中非零元素也２）动态地按增加出线支路数最少编号在上述的方法中，各节点的出线支路数是按原始网络统计出来的，在编号过程中认为固定不变的，事实上，在节点消去过程中，每消去一个节点以后，与该节点相连的各节点的出线支路数将发生变化(增加，减少或保持不变)。因此，如果每消去一个节点后，立即修正尚未编号节点的出线支路数，然后选其中支路数最少的一个节点进行编号，就可以预期得到更好的效果，动态按最少出线支路数编号方法的特点就是按出线最少原则编号时考虑了消去过程中各节点出线支路数目的变动情况。 3．MATLAB编程应用 Matlab 是“Matrix Laboratory”的缩写，主要包括：一般数值分析，矩阵运算、数字信号处理、建模、系统控制、优化和图形显示等应用程序。由于使用Matlab 编程运算与人进行科学计算的思路和表达方式完全一致，所以不像学习高级语言那样难于掌握，而且编程效率和计算效率极高，还可在计算机上直接输出结果和精美的图形拷贝，所以它的确为一高效的科研助手。 二、设计内容 1.设计流程图 实验名称：插值计算 1引言 在生产和科研中出现的函数是多种多样的。常常会遇到这样的情况：在某个实际问题中，虽然可以断定所考虑的函数f（x）在区间[a,b]上存在且连续，但却难以找到它的解析表达式，只能通过实验和观测得到在有限个点上的函数值。用这张函数表来直接求出其他点的函数值是非常困难的，在有些情况下，虽然可以写出f(x)的解析表达式，但由于结构十分复杂，使用起来很不方便。面对这些情况，构造函数P(x)作为f(x)的近似，插值法是解决此类问题比较古老却目前常用的方法，不仅直接广泛地应用与生产实际和科学研究中，而且是进一步学习数值计算方法的基础。 设函数y=f(x)在区间[a,b]上连续，且在n+1个不同的点a≤x0,x1……,xn≤b上分别取值y0,y1……,yn. 插值的目的就是要在一个性质优良、便于计算的函数φ中，求一简单函数P(x),使P(xi)=yi(i=0,1…,n)而在其他点x≠xi上，作为f(x)的近似。 通常，称区间[a,b]为插值区间，称点x0,x1,…,xn为插值节点，上式为插值条件，称函数类φ为插值函数类，称P(x)为函数f(x)在节点x0,x1,…,xn处的插值函数，求插值函数P(x)的方法称为插值法。 2实验目的和要求 用matlab定义分段线性插值函数、分段二次插值函数、拉格朗日插值函数，输入所给函 数表，并利用计算机选择在插值计算中所需的节点，计算f(0.15),f(0.31),f(0.47)的近似值。 3算法描述 1.分段线性插值流程图 2.分段二次插值流程图 3.拉格朗日插值流程图 4程序代码及注释 1.分段线性插值 %该程序仅针对《电力系统分析下》何仰赞P61 的4节点算例。 %节点电压用直角坐标表示时的牛顿-拉夫逊法潮流计算（matlab程序，可能有小错误，仅供学习之用，如果想要通用的程序，可自己动手改或再到pudn、csdn等网站搜索更好的）%南昌大学电力061 李圣涛2009年5月编写，2012年5月上传 clc %清空command windows clear all %清空workspace %为了提高可移植性、可读性、通用性，设置以下变量 N=4; %独立节点数 NPQ=2; %PQ节点数 NPV=1; %PV节点数 K=0; %迭代次数 %请输入最大迭代次数Kmax。可从0开始，以观察第Kmax次迭代的结果 Kmax=input('\n\n 请输入最大迭代次数后回车(可从零开始) Kmax=\n'); small=10^(-5); %ε不能太小 %i为节点标号，其中1号……NPQ号为PQ节点，（NPQ+1） %号……（N-1）号为PV节点，N号节点为平衡节点 %节点导纳矩阵Y的实部 G=[1.042093 -0.588235 0 -0.453858; -0.588235 1.069005 0 -0.480769; 0 0 0 0 ; -0.453858 -0.480769 0 0.934627 ]; %节点导纳矩阵Y的虚部 B=[ -8.242876 2.352941 3.666667 1.891074 ; 2.352941 -4.727377 0 2.403846 ; 3.666667 0 -3.3333333 0 ; 1.891074 2.403846 0 -4.261590 ]; %Y矩阵 Y=complex(G,B); %给定PQ节点的Pnode、Qnode，PV节点的Pnode、Vnode。（Vnode为节点电压的幅值）Pnode=[ -0.3 -0.55 0.5 ];%PQ、PV节点的初值P Qnode=[ -0.18 -0.13 0 ];%PQ节点的初值Q Vnode=[ 0 0 1.10 ];%PV节点的初值V %迭代初值 1电力系统潮流计算 潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角)、网络中的功率分布及功率损耗等。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性.可靠性和经济性。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。 2节点导纳矩阵的形成 在图1（a ）的简单电力系统中，若略去变压器的励磁功率和线路电容，负荷用阻抗表示，便可以得到一个有5个节点（包括零电位点）和7条支路的等值网络，如图1（b ）所示。将接于节点1和4的电势源和阻抗的串联组合变换成等值的电流源和导纳的并联组合，变得到图1（c ）的等值网络，其中1101I y E =和 4404I y E =分别称为节点1和4的注入电流源。 (a) 2 4? ?4 y (c) 图1 电力系统及其网络 以零电位点作为计算节点电压的参考点，根据基尔霍夫定律，可以写出4个独立节点的电流平衡方程如下： 10112121 12212022323242423323434244234434044()()()()0()()0()()y U y U U I y U U y U y U U y U U y U U y U U y U U y U U y U I ?+-=? -++-+-=? ? -+-=? ?-+-+=? （2-1） 上述方程组经过整理可以写成 1111221211222233244322333344422433444400Y U Y U I Y U Y U Y U Y U Y U Y U Y U Y U Y U Y U I ? + =? +++=? ? ++=?? ++=? （2-2） 式中， 111012 Y y y =+； 2220232412 Y y y y y =+++；332334Y y y =+；44402434Y y y y =++； 122112Y Y y ==-； 233223 Y Y y ==-； 244224 Y Y y ==-； 344334 Y Y y ==-。 一般的，对于有n 个独立节点的网络，可以列写n 个节点方程 11112211211222221122n n n n n n nn n n Y U Y U Y U I Y U Y U Y U I Y U Y U Y U I ? +++=? +++=? ? ? ?++ +=? （2-3） 也可以用矩阵写成 1111121212222212n n n n nn n n U I Y Y Y Y Y Y U I Y Y Y U I ???? ???????? ??????=?????? ?????? ?????????? （2-4） 或缩写为 YU I = （2-5） 矩阵Y 称为节点导纳矩阵。它的对角线元素ii Y 称为节点i 的自导纳，其值等于接于节点i 的所有支路导纳之和。非对角线元素 ij Y 称为节点i 、j 间的互导纳， 它等于直接接于节点i 、j 间的支路导纳的负值。若节点i 、j 间不存在直接支路，则有 ij Y =。由此可知节点导纳矩阵是一个稀疏的对称矩阵。 /*利用牛顿-拉夫逊迭代法(极坐标形式)，计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。所有参数应归算至标幺值下。 /*可计算最大节点数为100，可计算PQ,PV,平衡节点*/ /*可计算非标准变比和平行支路*/ #include （ 二 〇 一 四 年 十 二 月 课 程 论 文 学校代码： 10128 学 号: 20141100304 题 目：五节点系统计算机潮流计算编程 学生姓名：张佳羽 学 院：电力学院 系 别：电力系 专 业：电力系统及其自动化 指导教师：郭力萍 程序设计 % 本程序的功能是用牛顿拉夫逊法进行潮流计算 n=input('请输入节点数:n='); nl=input('请输入支路数:nl='); isb=input('请输入平衡母线节点号:isb='); pr=input('请输入误差精度:pr='); B1=input('请输入由各支路参数形成的矩阵:B1='); B2=input('请输入各节点参数形成的矩阵:B2='); X=input('请输入由节点号及其对地阻抗形成的矩阵：X='); Y=zeros(n); e=zeros(1,n);f=zeros(1,n);V=zeros(1,n); O=zeros(1,n);S1=zeros(nl); for i=1:nl if B1(i,6)==0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end Y(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; Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; 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); end for i=1:n S(i)=B2(i,1)-B2(i,2); B(i,i)=B(i,i)+B2(i,5); end P=real(S);Q=imag(S); ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; while IT2~=0 IT2=0;a=a+1; for i=1:n if i~=isb C(i)=0; D(i)=0; for j1=1:n C(i)= C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1); 2 模型简介及等值电路 2.1课程设计模型：模型3 电力网络接线如下图所示，各支路阻抗标幺值参数如下：Z12=0.02+j0.06，Z13=0.08+j0.24，Z23=0.06+j0.18，Z24=0.06+j0.12，Z25=0.04+j0.12，Z34=0.01+j0.03， Z45=0.08+j0.24，k=1.1。该系统中，节点1为平衡节点，保持 11.060 V j =+ &为定值；节点2、3、4都是PQ节点，节点5为PV节点，给定的注入功率分别为: 20.200.20 S j =+， 3-0.45-0.15 S j =， 40.400.05 S j =--， 50.500.00 S j =-+， 51.10 V= &。各节点电压（初值）标幺值参数如下： 节点 1 2 3 4 5 Ui（0）=ei（0） +jfi（0） 1.06+j0.0 1.0+j0.0 1.0+j0.0 1.0+j0.0 1.1+j0.0 计算该系统的潮流分布。计算精度要求各节点电压修正量不大于10-5。 图2-1 2.2模型分析 节点类型介绍 按变量的不同，一般将节点分为三种类型。 1 PQ节点 这类节点的有功功率和无功功率是给定的，节点(,) Vδ是待求量。通常变电所都是这一类型节点。由于没有发电设备，故其发电功率为零。有些情况下，系统中某些发电厂输出的功率在一段时间内是固定时，该发电厂母线也作为PQ节点。因此，电力系统中绝大多数节点属于这一类型。 2 PV节点 这类节点有功功率P和电压幅值V是给定的，节点的无功功率Q和电压的相位δ是待求量。这类节点必须有足够的可调无功容量，用以维持给定的电压幅值，因此又称为电压控制节点。一般选择有一定无功储备的发电厂和具有可调无功电源设备的变电所作为PV节点。 3 平衡节点 在潮流分布算出以前，网络中的功率损耗是未知的，因此，网络中至少有一个节点 基于MATLAB的直角坐标下牛顿-拉夫逊法潮流计算 摘要 潮流计算，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。它是基于配电网络特有的层次结构特性，论文提出了一种新颖的分层前推回代算法。该算法将网络支路按层次进行分类，并分层并行计算各层次的支路功率损耗和电压损耗，因而可大幅度提高配电网潮流的计算速度。论文在MATLAB环境下，利用其快速的复数矩阵运算功能，实现了文中所提的分层前推回代算法，并取得了非常明显的速度效益。另外，论文还讨论发现，当变压器支路阻抗过小时，利用Π型模型会产生数值巨大的对地导纳，由此会导致潮流不收敛。为此，论文根据理想变压器对功率和电压的变换原理，提出了一种有效的电压变换模型来处理变压器支路，从而改善了潮流算法的收敛特性。 关键词：电力系统；潮流分析；MATLAB Abstract Flow calculation is an important analysis function of power system and is the necessary facility of fault analysis, relay protection setting and security analysis. In addition, the traditional design method is a structured program design method based on functional decomposition, the entire software engineering as a combination of objects, as the domain of a particular issue, the composition of the object will remain basically unchanged Therefore, this decomposition method based on object design software structure relatively stable, easy to maintain and expand. . Combine the characteristics of power systems, software running on the use of MATLAB language WINDOWS OS graphical flow calculation software. The main features of the system are simple and intuitive graphical interface and stable operation. Calculated accurately Calculations, the algorithm has done a number of improvements to enhance the computing speed, the various types of effective package makes the procedure has good modularity maintainability and reusability. The MATLAB language is used to calculate flow distribution of power system in this paper. The typical examples explain that the method has the characteristics of simple programming high calculation efficiency and matching people habit the calculation result can satisfy the engineering calculation needs and at the same time verify the usefulness of the method. Key words: Electric power system; flow calculation; MATLAB 英文文献 Power Flow Calculation by Combination of Newton-Raphson Method and Newton’s Method in Optimization. Andrey Pazderin, Sergey Yuferev URAL STATE TECHNICAL UNIVERSITY – UPI E-mail: pav@https://www.sodocs.net/doc/4214641420.html,tu.ru, usv@https://www.sodocs.net/doc/4214641420.html,tu.ru Abstract--In this paper, the applicati on of the Newton’s method in optimization for power flow calculation is considered. Convergence conditions of the suggested method using an example of a three-machine system are investigated. It is shown, that the method allows to calculate non-existent state points and automatically pulls them onto the boundary of power flow existence domain. A combined method which is composed of Newton-Raphson method and Newton’s method in optimization is presented in the paper. Index Terms—Newton method, Hessian matrix, convergence of numerical methods, steady state stability Ⅰ.INTRODUCTION The solution of the power flow problem is the basis on which other problems of managing the operation and development of electrical power systems (EPS) are solved. The complexity of the problem of power flow calculation is attributed to nonlinearity of steady-state equations system and its high dimensionality, which involves iterative methods. The basic problem of the power flow calculation is that of the solution feasibility and iterative process convergence [1]. The desire to find a solution which would be on the boundary of the existence domain when the given nodal capacities are outside the existence domain of the solution, and it is required to pull the state point back onto the feasibility boundary, motivates to develop methods and algorithms for power flow calculation, providing基于极坐标的牛顿拉夫逊潮流计算
牛顿拉夫逊法潮流计算
基于MATLAB牛顿拉夫逊法进行潮流计算
牛顿插值法原理及应用
PQ分解法潮流计算实验
电力系统课程设计-牛顿拉夫逊法潮流计算
牛顿——拉夫逊法进行潮流计算
牛顿法潮流计算综述
计算方法实验报告 插值
直角坐标牛顿-拉夫逊法潮流计算matlab程序(仅供参考)
直角坐标系下牛顿法潮流计算
牛顿-拉夫逊迭代法极坐标潮流计算C语言程序
5节点电力系统牛顿-拉夫逊法潮流计算
基于极坐标的牛顿拉夫逊潮流计算修正版
基于MATLAB的直角坐标下牛顿拉夫逊法潮流计算
外文翻译--基于优化的牛顿——拉夫逊法和牛顿法的潮流计算
相关文档
- 电力系统三种潮流计算方法的比较
- 牛顿-拉夫逊算法(极坐标)潮流计算算例
- 牛顿-拉夫逊迭代法极坐标潮流计算C语言程序
- 牛顿拉夫逊法潮流计算matlab程序
- 牛顿法潮流计算流程图
- 牛顿法最优潮流汇编
- 基于牛顿拉夫逊法潮流计算的matlab实验报告
- 牛顿法潮流计算..
- 直角坐标系下牛顿法潮流计算
- PQ分解法潮流计算实验
- 直角坐标系下牛顿法潮流计算
- 牛顿法潮流计算
- 潮流计算的基本算法与使用方法
- 牛顿法潮流计算
- 牛顿法潮流计算Matlab程序
- 外文翻译--基于优化的牛顿——拉夫逊法和牛顿法的潮流计算
- 直角坐标系下牛顿法潮流计算
- 第三节牛顿 拉夫逊法潮流计算
- 牛顿拉夫逊法潮流计算
- P-Q分解法潮流计算