|
软件程序
软件/程序名称: |
牛顿拉夫讯 |
软件/程序大小: |
1k |
软件/程序语言: |
简体中文 |
运行平台: |
Windows XP/2003/Vista |
功能简介: |
牛顿拉夫讯 |
授权方式: |
免费版 |
研究/处理: |
汉化 |
马上加入,结交更多好友,共享更多资料,让你轻松玩转电力研学社区!
您需要 登录 才可以下载或查看,没有账号?立即加入
×
谁能帮忙看看这个牛顿拉夫讯怎么改进的,能用于配电网潮流计算?
" q& q6 P7 O& l# i牛顿拉夫讯
- u# e) ]& p8 O' Ydouble PowerFlow(int iNode, int iBalc, ROAD *pRoad, NodePQ *pNodePQ, bool vOutputVoltage = false)
l ?) N& i& F+ f{
- ?9 v, }3 }: K# g; c1 E int M=2*(iNode-1);
' k Q) K2 R0 b' w1 J7 e int flag1[N]; // PV节点标志
. f; R) c5 f' E- ^ 2 K: T& B$ v& W' |! i5 ]8 @
double g[N][N]; // 电导' q4 p- W+ @- x* j
double b[N][N]; // 电纳
4 z( H; o3 p9 T+ d% I) c: Q double dvv[N]; // 电压微分- h* ]7 `7 P" x& @" z
double angle[N]; // 相角
0 `8 U) Y( [3 k ] double p[N],dp[N]; // 有功及其微分
4 W/ T9 N9 X; L$ L: m$ I& J double q[N],dq[N]; // 无功及其微分
, o+ a- Q5 e* D) x. Y1 j double e[N]; // 节点电压实部
4 U1 M+ J! Q! O* k# Q3 o double f[N]; // 节点电压虚部+ A: G: g2 V. {. B
double *jcbb = new double[M];
; J$ P; ~: T4 G( R5 ^ double *jcb = new double [M][M];
# e8 j# S; q$ R5 s8 i v
# ^! V T$ {8 B8 `6 V0 V8 p8 Q& I int i,j,i1,j1,n,m,t,k1,balc,node/*,flag,npv*/;
3 l8 a+ c" b9 K4 X, _- u double rr,xx,bb,k,com1,com2,sump=0.0,sumq=0.0;- t, i9 ]$ H) ~0 ~. L2 w, F4 U
double lossp=0.0, lossq=0.0;4 c$ M+ o0 d% t! ~; c& l
node = iNode; // 节点数8 n8 k' Z# {6 k
balc = iBalc - 1; // 平衡节点% o1 [# }# _: g. I$ ^) {
for(i=0;i<node;i++){
' N/ S1 I- Q; Q7 F; e flag1[i]=0;/ z; @0 K- d" W; a
dp[i]=0.0;
2 k3 W; f8 v' b+ m dq[i]=0.0;
0 A3 `/ F- } P dvv[i]=0.0;
3 E1 g& V# S2 \; c; n4 t }' s9 u; g- A* u* k. ^
- L' U& l* ]4 V! f2 s& H
for(i=0;i<node;i++) {
, y; X7 Q' z6 A& X+ S6 }7 k for(j=0;j<node;j++){
, O" S# B4 f' V& j; B0 ` g[i][j]=0.0;
& V! M p5 k! I/ D9 w8 Q b[i][j]=0.0;
- s$ ~4 B, r, o {& N8 ^ }
# ^ I/ K+ ~6 ?3 R }- K3 M+ q0 |/ V1 q5 P$ I
for(n = 0; n < node-1; n++){
/ K4 T V$ S9 L, g9 J i = pRoad[n].head;
' ?! i! m2 o( W9 ]; i j = pRoad[n].tail; |' R, @4 z* F- N4 i% x
rr = pRoad[n].r;# y$ _+ B6 z3 g1 `* }5 `9 [( {
xx = pRoad[n].x;
{7 A u l* B, e! @0 a- Q: G// printf("%d\t%d\t%.3f\t%.3f\n", pRoad[n].head, pRoad[n].tail, pRoad[n].r, pRoad[n].x);
" u1 ?6 a2 @+ V* f7 Q bb = 0;8 G" h8 G- _' n5 N {7 Y0 p! s
k = 1;
. F% d% \* w! r ) L6 w: Z# f2 U* S6 G+ |
g[i-1][j-1]+=-(rr/(rr*rr+xx*xx)/k);
5 Q$ U6 k9 d: i W b[i-1][j-1]+=-(-xx/(rr*rr+xx*xx)/k);
2 w5 j* f4 j" k0 w8 Q g[j-1][i-1]+=-(rr/(rr*rr+xx*xx)/k);7 G6 u" e. z9 l$ B1 g
b[j-1][i-1]+=-(-xx/(rr*rr+xx*xx)/k);
$ v% o6 E& x& h2 I- A" Y g[i-1][i-1]+=rr/((rr*rr+xx*xx)*k*k);/ J/ C3 w, T) f7 F$ z1 h
g[j-1][j-1]+=rr/(rr*rr+xx*xx);
3 Z/ _& P1 P0 N: r% } b[i-1][i-1]+=bb+(-xx/((rr*rr+xx*xx)*k*k));
1 q" A/ x. {& ~; n& o- A0 T b[j-1][j-1]+=bb+(-xx/(rr*rr+xx*xx));7 k5 o6 _0 q& C* l. z
}
4 L' w8 q7 }# o' R4 j for(i=0;i<node;i++){
7 q. h5 e& Y9 q$ | angle[i] = 0;
! t1 O9 D4 Z' i5 i. { v[i] = 1; i; | u: j2 j( w0 l
p[i] = pNodePQ[i].p;
7 F9 e: }: }3 m q[i] = pNodePQ[i].q;( n( V) n4 A( F; O5 {) u' W
// printf("P:%.3f\t Q:%.3f\n", pNodePQ[i].p, pNodePQ[i].q);3 F; v2 a; P# _7 U0 Q& @6 w' r
e[i]=v[i]*cos(angle[i]*PI/180);9 m; E1 W& J5 k$ C$ S/ @
f[i]=v[i]*sin(angle[i]*PI/180);
) t1 V. A0 i0 Q# @ }3 m' t0 Z1 X9 q. x$ x
// flag1[1]=1;& @* q5 R1 ~" d" |; |
. S! a! [) X: ^
// printf("共有n=%d条支路\n",n);( W7 D x \5 f3 g' Q: y3 _9 ?! _
// printf("节点导纳矩阵如下:\n");
- }6 D4 _- l7 {% j [/ i + r( B5 A( o, x" i) K! P9 O1 r" L" z$ |
// for(i=0;i<node;i++){* e0 A2 ~ A: l& r
// for(j=0;j<node;j++)
0 L5 `2 M2 x, b8 b; A// printf("%lf+j%lf ", g[i][j],b[i][j]);
9 u" v3 q+ J9 M; m, W* F// printf("\n\n");}5 p. _' ]- j. X# a6 W: Y" }
//
* O9 b6 F6 ]5 o' A% {/ G, p$ N// printf("初始值如下:\n");+ n! v: @7 @8 t. c$ {
// printf("节点号 电压 相角 有功 无功 e f\n");
, a9 X% x' P% O: c// for(i=0;i<node;i++)
3 T w6 S- l$ M8 G" i; g: F8 } p// printf("%d %lf %lf %lf %lf %lf %lf\n", (i+1), v[i], angle[i], p[i], q[i], e[i], f[i]);7 Q4 s% V$ i: ?2 o( E7 k1 x* ]
6 M' M1 i6 [6 e# q) T5 n for(t=0;t<50;t++){( ~( g" X+ K! f& ^! }' m; o- ^$ E! \
// printf("\n第%d次迭代开始:\n", t);5 l3 M5 v1 d) x3 ?! Y: k
for(i=0;i<node;i++){
2 t' F5 q9 e* e$ t1 T! T if(i==balc) /*如果是平衡节点*/) F2 ], j/ ]3 M' ?
continue;9 @" H+ G" e' U9 T
else{
- A, H. I$ j' Z( O4 X com1=0.0, com2=0.0;* O9 K# y& v% b$ q; U2 a; A, C
for(j=0;j<node;j++){* }( Q! S' m. {# ~
com1+=g[i][j]*e[j]-b[i][j]*f[j];
q3 K2 O/ }) q8 A com2+=g[i][j]*f[j]+b[i][j]*e[j];
/ j* o ^" y. y; b2 H$ E( u }( a% M5 v: K. G R( Q F; J
if(flag1[i]!=1){ /*如果是PQ节点*/) s, G* N3 y& ^
dp[i]=p[i]-e[i]*com1-f[i]*com2;
/ m4 b* ]) o' o3 Q% F dq[i]=q[i]-f[i]*com1+e[i]*com2;2 x, i# G% W; [8 l4 G
}8 L( q" c0 g$ X( B
else{$ _. E8 c9 p' V: {
dp[i]=p[i]-e[i]*com1-f[i]*com2;' F, }+ d) Q2 l, R0 m3 {$ Q
dvv[i]=v[i]*v[i]-(e[i]*e[i]+f[i]*f[i]);, b' J9 [/ K! k
}: r" d( f# K" ^7 t3 A/ e; h
}! j, t. h! @# X0 f
}
6 g) D: y. B1 j! z
# O6 S/ X. C2 K4 T; S' x! _ for(i=0,j=0;i<node;i++){
/ D9 N$ S8 q( z% \% k/ l; W if(i==balc) continue;
3 Z5 s( D0 j7 A0 L/ k/ i else if(flag1[i]!=1){
: I/ K7 L1 p5 I/ D* C/ ~' T jcbb[j++]=dp[i];' J* O8 l# X1 F9 w& H1 M
jcbb[j++]=dq[i];; z$ u# @1 _; P
// printf("\ndp[%d]=%lf\ndq[%d]=%lf",(i+1),dp[i],(i+1),dq[i]);6 b5 d3 \2 ]6 y% U8 Q8 s
}% Q2 ^! Q- A* w2 H. ^8 U2 r
else{
. S; B; }! u% k5 H. ?0 Y, q( E jcbb[j++]=dp[i];
# a5 m+ q1 d, C/ U* A7 d jcbb[j++]=dvv[i];0 H2 R8 I* W+ C* R: g. n, H+ K
// printf("\ndp[%d]=%lf\ndvv[%d]=%lf",(i+1),dp[i],(i+1),dvv[i]);
% K& o F5 H: n+ o* p# R }
& D% H3 p9 o$ _- E- g7 g }; u/ N3 ~3 ]. I' i3 R1 p6 L
// printf("\n");
8 g1 y( f% \) q# X) I j - J+ f/ g$ F) J) F3 P6 z
for(i=0;i<node;i++) {
4 B7 r: x4 {# H* I if(fabs(dp[i])>=E||fabs(dq[i])>=E||fabs(dvv[i])>=E)
t* ^3 p5 ^3 E7 S+ d6 ]/ q$ b break;; c* @# m7 g* i7 C% U7 O
}
0 P1 b# _* @. C3 ]4 w9 |. ]4 a if(i>=node)
% u4 E. g1 M; t break; 3 h( q | N8 Y6 `/ J
: i0 U( G9 r M6 W! q
5 n4 A: Z( X4 Q4 C7 d) _% O# x for(i=0,i1=0;i<node;i++){
$ d$ N, l2 ? y" ] if(i==balc)
6 V+ W$ \9 c, b continue;
4 i# H2 E! O$ X i1++;
* i! C) H- b; {4 v( x n=2*i1;% B4 f5 P6 L; R3 p( B3 L0 @
for(j=0,j1=0;j<node;j++){" Q' O0 t+ J2 X8 A6 q# k8 T
if(j==balc) continue;
: _; @' u* ~9 @" X' S/ u8 g! P j1++;
3 W- o) [; D& Q: @2 W m=j1*2;
& {3 r2 E* L% O! t2 M2 T if(i!=j){ $ e! D6 E( d$ i9 @( U- _/ l$ `
jcb[n-2][m-2]=-(g[i][j]*e[i]+b[i][j]*f[i]);
, f2 F5 K. ^3 I& u jcb[n-2][m-1]=b[i][j]*e[i]-g[i][j]*f[i];) s1 C& q0 H8 ?
if(flag1[i]!=1)% n" B3 q" P6 ~# m% L+ P# ]
{# K0 o% L1 b9 Y# C4 X g
jcb[n-1][m-2]=jcb[n-2][m-1];
, O; R. N4 L/ |& b- p jcb[n-1][m-1]=-jcb[n-2][m-2];
( K0 X9 W2 i8 u. ^4 [/ ~5 @ }
; l1 J$ J5 w5 n% S! ^" S else1 M0 g D5 L0 A# t+ V
jcb[n-1][m-2]=jcb[n-1][m-1]=0.0;
/ j% u5 }$ J& U" W }
8 z# s3 c; I8 T else{; B# r/ x' O* M6 L+ v, ^: u; e
com1=0.0, com2=0.0;7 C6 V$ w! F" [3 x; w
for(k1=0;k1<node;k1++)
$ f2 h( ~# \" E! k# U+ q! E {5 V2 B: _. m7 ~9 I1 c6 Q
com1+=g[i][k1]*e[k1]-b[i][k1]*f[k1];
' F3 } Y* W- p- H) g4 z com2+=g[i][k1]*f[k1]+b[i][k1]*e[k1];
2 m+ [8 O2 [# _+ I/ ~+ i }
2 ]" ?; N3 T, N6 m. X3 G jcb[n-2][m-2]=-com1-g[i][i]*e[i]-b[i][i]*f[i];5 e" G9 p. ^& {
jcb[n-2][m-1]=-com2+b[i][i]*e[i]-g[i][i]*f[i];& K+ C* l9 o, @) h' J; z. k
if(flag1[i]!=1)
0 A3 u( M3 R- b2 [! H- I {9 U4 P1 \9 @9 b9 A2 ]/ `
jcb[n-1][m-2]=com2+b[i][i]*e[i]-g[i][i]*f[i];! ^2 t3 r- h. o' D' C4 v, B3 c6 ?
jcb[n-1][m-1]=-com1+g[i][i]*e[i]+b[i][i]*f[i];
% v4 Y( r: w c' L. O. ? }
9 f, M- H K, c0 R0 {2 m else( v8 o: O0 K6 s2 {; b
{8 a& _$ p6 c8 u2 o, S0 l+ N
jcb[n-1][m-2]=-2*e[i];
; M2 _$ X; o" L& o6 k3 O9 G1 ?% n jcb[n-1][m-1]=-2*f[i]; g% C `! L6 n2 m, {. l
}
6 V4 J" }/ i9 D5 o }
8 v1 p+ M. e; ]( a3 J7 n; ]+ F+ A }
2 O- S7 N3 o% S }
0 @1 ?6 T$ ]5 T0 X' h% i l% B) V
2 t3 C4 h9 O( w% C- F @// printf("\n第%d次输出-J矩阵如下:\n", t);
! Z" r. _. X2 ^ for(i=0;i<2*(node-1);i++){
* `* i3 w, c% q2 p- ?" v for(j=0;j<2*(node-1);j++){* q0 L) ], m2 l8 \# I/ a" Y
jcb[i][j]*=-1.0;2 {! ]* K/ { D3 N- Y' |/ v
// printf("%lf ", jcb[i*M+j]);
" G, N+ C& W# f/ Z2 w! B }
8 E" m. s/ d0 P& Y: \// printf("\n\n");' B# m; _& ?) N' H4 h% X+ S4 ~9 N0 [
}0 U# t7 C6 ]& v5 a6 |
4 ~7 [' i9 F2 y$ {7 J7 p' r /*算线形方程组*/
6 c/ X. s; M( G g9 |) ]+ [ if (agaus(&jcb[0],&jcbb[0],2*(node-1))!=0) {" U: y8 N! ~) i5 F5 S6 n+ I4 w
for (i=0,j=0;i<2*(node-1);i++) {
V/ ~7 g& M- X if(j==balc) 1 s: k( O) ] z5 D- C+ \+ s
j++;: R5 m7 ? M1 I! T! J
# o( {; B0 {& t0 u& {7 s# u9 p1 m if(i%2==0){2 ]/ N7 `6 u4 H) I" p' x0 }
e[j]+=jcbb[i];4 h, N: A2 I+ {' t2 Q" q
// printf("de[%d]=%lf e[%d]=%lf\n", (j+1),jcbb[i],(j+1),e[j]);
+ H, {) o. I9 Q0 @3 b }
s$ _3 n' t% } else{
a9 X0 d3 u8 p _) P" d f[j]+=jcbb[i];+ F4 N8 w9 M5 j+ ^
// printf("df[%d]=%lf f[%d]=%lf\n",(j+1),jcbb[i],(j+1),f[j]);
7 f! E. S/ ] ?- s j++;
# c" I* K: o5 K3 o1 W& u }0 t! K' Z7 K5 y" U' G% I
}' M" {# x% p5 i9 m9 r+ }- ~3 }
}
) K2 g/ G* b0 p }
; V/ u$ {4 D" p & ^5 n! I5 L: H! K' T4 i2 L7 s
// printf("共经过了%d次迭代\n", t);" g. f/ |1 _" H! O
* g& \; }7 u0 R' j- Q /*求平衡节点P,Q*/! j) I$ r) Z% t4 Y" J- G- Y0 J# r
i=balc;2 _& O& H8 I& d5 c
com1=0.0, com2=0.0;. d6 D& m- v7 d I+ v
for(j=0;j<node;j++){
! o) z5 s. h$ |& \: B5 H5 |" F, |- a2 K com1+=g[i][j]*e[j]-b[i][j]*f[j];; D# V0 s% V. t
com2+=g[i][j]*f[j]+b[i][j]*e[j];' T9 l" d7 H5 t( c8 w( I8 O1 H
}
* s! V; f1 u9 ^: h% ^, `: A; z p[i]=e[i]*com1+f[i]*com2;' A. k. q) A# z6 _/ f1 A
q[i]=f[i]*com1-e[i]*com2;1 s J$ W$ K( ^$ |
W) w j# q$ J( L6 O
/*求其他各节点P,Q*/
1 Q$ {/ w6 V/ o" i$ r for(i=0;i<node;i++){
7 s4 t4 S! P4 i& U) \$ m# W if(i==balc) continue;
4 l, @% V0 f; K6 ~8 { else
1 }2 [7 K+ ?% r& \ {
2 m0 \: O0 C) S com1=0.0, com2=0.0;
6 E; ?6 ^& O: f# g" {: R( L' k for(j=0;j<node;j++)
4 j. B4 A0 A* X {
0 l0 j# x( a: t7 L com1+=g[i][j]*e[j]-b[i][j]*f[j]; t0 l3 r/ V8 o4 a$ }
com2+=g[i][j]*f[j]+b[i][j]*e[j];
" h9 D$ a8 k: `1 i2 w) C }
, h! e& l! F* }. a: x& S' j8 K p[i]=e[i]*com1+f[i]*com2;
* R9 o6 E9 D4 F- }+ G: w# O q[i]=f[i]*com1-e[i]*com2;
, z% f- N9 c4 p5 y' v* h6 C& E }
+ A. F9 d* q$ ]" B }
( @/ C) i$ S5 ?2 o- d! V # F: s5 a4 D; j- e' }! M8 f
// printf("最终结果为:\n");& i! e: T6 `! g1 W0 V
// printf("节点号 电压 相角 有功 无功 e f\n");, {( w3 y' j0 q) S" h% m/ x; P, [
for(i=0;i<node;i++)
6 w1 B6 Y/ {2 j) X2 J9 F/ c {
: {! {9 i3 `! X: G" w7 c- _& | v[i]=sqrt(e[i]*e[i]+f[i]*f[i]);
. v9 z7 z6 p8 q5 w: j. M5 l angle[i]=atan(f[i]/e[i])*180/PI;( A+ t$ Q: L* s& L7 T0 C. e0 ?
sump+=p[i];
) b; h3 q e: j) `$ b sumq+=q[i];
, P9 V1 C6 J1 t. I/ V" o+ \// printf("%d %lf %lf %lf %lf %lf %lf\n", (i+1), v[i], angle[i], p[i], q[i], e[i], f[i]);5 l! O) g, Q, S9 ?
} |
|