当前位置：文档库 › 基于牛顿拉夫逊法潮流计算的matlab实验报告

基于牛顿拉夫逊法潮流计算的matlab实验报告

一、实验目的

应用MATLAB 语言编写具有一定通用性的牛顿-拉夫逊法潮流计算程序。要求：

（1）潮流计算方法为牛顿-拉夫逊法。（2）编程语言为MATLAB 。（3）程序具有较强通用性。二、程序流程图

所用公式

12222[()()]

[()()]()

j n

i i i ij j ij j i ij j ij j j j n i i i ij j ij j i ij j ij j j i i i i

P P e G e B f f G f B e Q Q f G e B f e G f B e U U e f ====?

?=--++??

?=---+??

??=-+???

∑∑

i j ≠

200

i ij ij i ij i

i i ij

ij i ij i ij i

i ij ij i ij i ij

i i ij i i ij i P H B e G f f P N G e B f e Q J G e B f N f Q L B e G f H e U R f U S e ??==-+?????==+???

??==--=-???

???==-+=?????==???

==???

i j

222i ii ij i ij i ii

i ii ij i ij i ii

i i ii

ij i ij i ii i

i ii ij i ij i ii i i ii i

i i ii i

i P H B e G f b f P N G e B f a e Q J G e B f a f Q L B e G f b e U R f

f U S e e ??

==-++?????==++???

??==--+????

??==-+-?????==???

==???

其中

()

j n

ii ii i ii i ij j ij j j j i j n

ii i ii i ij j ij j

j j i a G e B f G e B f b G f B e G f B e ==≠==≠?

=-+-?????=+++???

∑∑

三、求解问题及其结果

开始

形成节点导纳矩阵输入原始数据设节点电压(0)(0)i i e f ，i=1,2…,n,i ≠s 置迭代次数0k = 置节点号i=1 按式（3-3），（3-4）计算雅克比矩阵元素按式（3-2）计算PQ 节点的()k i P ?，()k i Q ?，PV 节点的()k i P ?，()2k i

U ?

求解修正方程式，得()k i e ?，()k i f ?

雅克比矩阵是否已全部形成？计算平衡节点及PV 节点功率求()max ||k e ?，()max ||k f ?

迭代次数 k=k+1 i=i+1 ()()max max ||,||k k e f ε??≤？潮流计算完成计算各节点电压的新值： (1)()()k k k i e e e +=+? (1)()()k k k i f f f +=+?

IEEE-美国新英格兰10机39节点测试系统

一、系统单线图

G G G

225

3729

5418

3624

2034

3310

632

二、系统参数

1）系统容量基准值为100MV A。

2) 负荷数据见表D-1

表D-1 负荷数据

节点号有功（MW）无功（Mvar）节点号有功（MW）无功（Mvar）

BUS- 3 BUS- 4 BUS- 7 BUS- 8 BUS- 12 BUS- 15 BUS- 16 BUS- 18 BUS- 20 322.0

500.0

233.8

522.0

8.5

320.0

329.0

158.0

680.0

2.4

184.0

84.0

176.0

88.0

153.0

32.3

30.0

103.0

BUS- 21

BUS- 23

BUS- 24

BUS- 25

BUS- 26

BUS- 27

BUS- 28

BUS- 29

BUS- 39

274.0

247.5

308.0

224.0

139.0

281.0

206.0

283.5

1104.0

115.0

84.6

-92.2

47.2

17.0

75.5

27.6

26.9

3）发电机数据见表D-2

表D-2 发电机数据

发电机节点号有功（MW）无功（Mvar）电压（p.u.）

PV PQ PQ PV PV PQ PV PV PV BUS-30

BUS-32

BUS-33

BUS-34

BUS-35

BUS-36

BUS-37

BUS-38

BUS-39

250.00

650.00

632.00

508.00

650.00

560.00

540.00

830.00

1000.00

175.90

103.35

96.88

1.04750

1.01230

1.04930

1.02780

1.02650

1.03000

平衡节点BUS-31 0.0(电压幅角) 1.0（幅值）

4）线路参数见表D-3

表D-3 线路参数

序号节点I 节点J R(p.u.) X(p.u.) B/2(p.u.)

LN1 LN2 LN3 LN4 LN5 LN6 LN7 LN8 LN9 BUS-2

BUS-39

BUS-3

BUS-25

BUS-4

BUS-18

BUS-5

BUS-14

BUS-6

BUS-1

BUS-2

BUS-3

BUS-4

BUS-5

0.00350

0.00100

0.00130

0.00700

0.00130

0.00110

0.00080

0.00020

0.04110

0.02500

0.01510

0.00860

0.02130

0.01330

0.01280

0.01290

0.00260

0.34935

0.3750

0.12860

0.07300

0.11070

0.10690

0.06710

0.06910

0.02170

LN10 LN11 LN12 LN13 LN14 LN15 LN16 LN17 LN18 LN19 LN20 LN21 LN22 LN23 LN24 LN25 LN26 LN27 LN28 LN29 LN30 LN31 LN32 LN33 LN34 BUS-8

BUS-7

BUS-11

BUS-8

BUS-9

BUS-39

BUS-11

BUS-13

BUS-14

BUS-15

BUS-16

BUS-17

BUS-19

BUS-21

BUS-24

BUS-18

BUS-27

BUS-22

BUS-23

BUS-24

BUS-26

BUS-27

BUS-28

BUS-29

BUS-5

BUS-6

BUS-7

BUS-8

BUS-9

BUS-10

BUS-13

BUS-14

BUS-15

BUS-16

BUS-17

BUS-21

BUS-22

BUS-23

BUS-25

BUS-26

BUS-28

0.00080

0.00060

0.00070

0.00040

0.00230

0.00100

0.00040

0.00090

0.00180

0.00090

0.00070

0.00160

0.00080

0.00030

0.00070

0.00130

0.00080

0.00060

0.00220

0.00320

0.00140

0.00430

0.00570

0.00140

0.01120

0.00920

0.00820

0.00460

0.03630

0.02500

0.00430

0.01010

0.02170

0.00940

0.00890

0.01950

0.01350

0.00590

0.00820

0.01730

0.01400

0.00960

0.03500

0.03230

0.01470

0.04740

0.06250

0.01510

0.07380

0.05650

0.06945

0.03900

0.19020

0.60000

0.03645

0.08615

0.18300

0.08550

0.06710

0.15200

0.12740

0.03400

0.06595

0.16080

0.12825

0.09230

0.18050

0.25650

0.11980

0.39010

0.51450

0.12450

LN35: BUS-4接有并联电容器，B4＝1.0000

LN36: BUS-5接有并联电容器，B4＝2.0000

5)变压器参数见表D-4

表D-4 变压器参数

序号节点I 节点J R(p.u.) X(p.u.) 变比(%)

TR37 TR38 TR39 TR40 TR41 TR42 TR43 TR44 TR45 TR46 TR47 TR48 BUS-11

BUS-13

BUS-30

BUS-31

BUS-32

BUS-34

BUS-33

BUS-35

BUS-36

BUS-37

BUS-38

BUS-20

BUS-12

BUS-2

BUS-6

BUS-10

BUS-20

BUS-19

BUS-22

BUS-23

BUS-25

BUS-29

BUS-19

0.00160

0.00000

0.00090

0.00070

0.00000

0.00050

0.00060

0.00080

0.00070

0.04350

0.01810

0.02500

0.02000

0.01800

0.01420

0.01430

0.02720

0.02320

0.01560

0.01380

100.6

102.5

107.0

100.9

107.0

102.5

100.0

102.5

106.0

%IEEE-美国新英格兰10机39节点测试系统% 1 2 3 4 5 6 % bus volt angle p q type bus=[ 1 1.0000 0.00 0.00 0.00 1

2 1.0000 0.00 0.00 0.00 1

3 1.0000 0.00 -3.22 -0.02

4 1

4 1.0000 0.00 -5.00 -1.84 1

5 1.0000 0.00 0.00 0.00 1

6 1.0000 0.00 0.00 0.00 1

7 1.0000 0.00 -2.338 -0.84 1

8 1.0000 0.00 -5.22 -1.76 1

9 1.0000 0.00 0.00 0.00 1

10 1.0000 0.00 0.00 0.00 1

11 1.0000 0.00 0.00 0.00 1

12 1.0000 0.00 -0.085 -0.88 1

13 1.0000 0.00 0.00 0.00 1

14 1.0000 0.00 0.00 0.00 1

15 1.0000 0.00 -3.20 -1.53 1

16 1.0000 0.00 -3.29 -0.323 1

17 1.0000 0.00 0.00 0.00 1

18 1.0000 0.00 -1.58 -0.30 1

19 1.0000 0.00 0.00 0.00 1

20 1.0000 0.00 -6.80 -1.03 1

21 1.0000 0.00 -2.74 -1.15 1

22 1.0000 0.00 0.00 0.00 1

23 1.0000 0.00 -2.475 -1.15 1

24 1.0000 0.00 -3.08 -0.922 1

25 1.0000 0.00 -2.24 -0.472 1

26 1.0000 0.00 -1.39 -0.17 1

27 1.0000 0.00 -2.81 -0.755 1

28 1.0000 0.00 -2.06 -0.276 1

29 1.0000 0.00 -2.835 -0.269 1

30 1.0475 0.00 2.50 0.00 2

31 1.0000 0.00 0.00 0.00 3

32 1.0000 0.00 6.50 1.759 1

33 1.0000 0.00 6.32 1.0335 1

34 1.0123 0.00 5.08 0.00 2

35 1.0493 0.00 6.50 0.00 2

36 1.0000 0.00 5.60 0.9688 1

37 1.0278 0.00 5.40 0.00 2

38 1.0265 0.00 8.30 0.00 2

39 1.0300 0.00 -1.04 0.00 2];

% 1 2 3 4 5 6 7 % line: from bus to bus R, X, G, B/2 K

line=[ 2 1 0.00350 0.04110 0 0.34935 0;

39 1 0.00100 0.02500 0 0.37500 0;

3 2 0.00130 0.01510 0 0.12860 0;

25 2 0.00700 0.00860 0 0.07300 0;

4 3 0.00130 0.02130 0 0.11070 0;

18 3 0.00110 0.01330 0 0.10690 0;

5 4 0.00080 0.01280 0 0.06710 0;

14 4 0.00080 0.01290 0 0.06910 0;

6 5 0.00020 0.00260 0 0.02170 0;

8 5 0.00080 0.01120 0 0.07380 0;

7 6 0.00060 0.00920 0 0.05650 0;

11 6 0.00070 0.00820 0 0.06945 0;

8 7 0.00040 0.00460 0 0.03900 0;

9 8 0.00230 0.03630 0 0.19020 0;

39 9 0.00100 0.02500 0 0.60000 0;

11 10 0.00040 0.00430 0 0.03645 0;

13 10 0.00040 0.00430 0 0.03645 0;

14 13 0.00090 0.01010 0 0.08615 0;

15 14 0.00180 0.02170 0 0.18300 0;

16 15 0.00090 0.00940 0 0.08550 0;

17 16 0.00070 0.00890 0 0.06710 0;

19 16 0.00160 0.01950 0 0.15200 0;

21 16 0.00080 0.01350 0 0.12740 0;

24 16 0.00030 0.00590 0 0.03400 0;

18 17 0.00070 0.00820 0 0.06595 0;

27 17 0.00130 0.01730 0 0.16080 0;

22 21 0.00080 0.01400 0 0.12825 0;

23 22 0.00060 0.00960 0 0.09230 0;

24 23 0.00220 0.03500 0 0.18050 0;

26 25 0.00320 0.03230 0 0.25650 0;

27 26 0.00140 0.01470 0 0.11980 0;

28 26 0.00430 0.04740 0 0.39010 0;

29 26 0.00570 0.06250 0 0.51450 0;

29 28 0.00140 0.01510 0 0.12450 0;

4 0 0 0 0 1.0000 0;

5 0 0 0 0 2.0000 0;

11 12 0.00160 0.04350 0 0 100.60000/100;

13 12 0.00160 0.04350 0 0 100.60000/100;

30 2 0.00000 0.01810 0 0 102.50000/100 ;

31 6 0.00000 0.02500 0 0 107.00000/100 ;

32 10 0.00000 0.02000 0 0 107.00000/100 ;

34 20 0.00090 0.01800 0 0 100.90000/100 ;

33 19 0.00070 0.01420 0 0 107.00000/100 ;

35 22 0.00000 0.01430 0 0 102.50000/100 ;

36 23 0.00050 0.02720 0 0 100.00000/100 ;

37 25 0.00060 0.02320 0 0 102.50000/100 ;

38 29 0.00080 0.01560 0 0 102.50000/100 ;

20 19 0.00070 0.01380 0 0 106.00000/100] ;

计算结果牛顿－拉夫逊法潮流计算结果节点计算结果：

n节点节点电压节点相角（角度）节点注入功率

1 1.049185 -8.874991 0.000000 + j 0.000000

2 1.053167 -6.367180 0.000000 + j 0.000000

3 1.041493 -9.207297 -3.220000 + j -0.024000

4 1.036574 -10.04258

5 -5.000000 + j -1.840000

5 1.044652 -8.959237 0.000000 + j 0.000000

6 1.043883 -8.293104 0.000000 + j 0.000000

7 1.032645 -10.342431 -2.338000 + j -0.840000

8 1.031177 -10.811816 -5.220000 + j -1.760000

9 1.042715 -10.595648 0.000000 + j 0.000000

10 1.046426 -6.010476 0.000000 + j 0.000000

11 1.044322 -6.792462 0.000000 + j 0.000000

12 1.030736 -6.795388 -0.085000 + j -0.880000

13 1.042351 -6.675491 0.000000 + j 0.000000

14 1.036310 -8.232337 0.000000 + j 0.000000

15 1.018517 -8.519794 -3.200000 + j -1.530000

16 1.025492 -7.051856 -3.290000 + j -0.323000

17 1.032750 -8.077118 0.000000 + j 0.000000

18 1.034779 -8.936485 -1.580000 + j -0.300000

19 1.044862 -2.382169 0.000000 + j 0.000000

20 0.988148 -3.811032 -6.800000 + j -1.030000

21 1.024926 -4.596980 -2.740000 + j -1.150000

22 1.042650 -0.070512 0.000000 + j 0.000000

23 1.032952 -0.245457 -2.475000 + j -1.150000

24 1.021125 -6.906503 -3.080000 + j -0.922000

25 1.060163 -4.952002 -2.240000 + j -0.472000

26 1.052697 -6.205207 -1.390000 + j -0.170000

27 1.037683 -8.217337 -2.810000 + j -0.755000

28 1.050444 -2.695196 -2.060000 + j -0.276000

29 1.050163 0.063077 -2.835000 + j -0.269000

30 1.004392 1.594781 6.500000 + j 1.759000

31 0.991632 2.892572 6.320000 + j 1.033500

32 1.050539 7.797786 5.600000 + j 0.968800

33 1.047500 -3.957598 2.500000 + j 1.211174

34 1.012300 1.385774 5.080000 + j 1.826359

35 1.049300 4.925324 6.500000 + j 2.637566

36 1.027800 1.819476 5.400000 + j -0.108224

37 1.026500 7.125579 8.300000 + j 0.214225

38 1.030000 -10.390696 -1.040000 + j -2.291639

39 1.000000 0.000000 5.628660 + j 1.384403

线路计算结果：

n节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)

2 1 1.178698 + j -0.360055 -1.174311 + j -0.360481 0.004386 + j -0.720536 39 1 6.405845 + j -2.096152 -6.361848 + j 2.408287 0.043997 + j 0.312135

3 2 -3.633961 + j -0.542613 3.649983 + j 0.446577 0.016021 + j -0.096036 25 2 2.370242 + j -1.109311 -2.328681 + j 0.997356 0.041562 + j -0.111955

4 3 -0.750370 + j -0.307172 0.751094 + j 0.080014 0.000724 + j -0.227159 18 3 0.337560 + j -0.66385

5 -0.337133 + j 0.438599 0.000427 + j -0.22525

6 5 4 1.635254 + j 0.499000 -1.633054 + j -0.609119 0.002200 + j -0.110119 14 4 2.621711 + j -0.216428 -2.616576 + j 0.15077

7 0.005135 + j -0.065651 6 5 4.826035 + j -0.675350 -4.821682 + j 0.684607 0.004353 + j 0.009257

8 5 -3.178130 + j -1.041836 3.186428 + j 0.998989 0.008297 + j -0.042847 7 6 -4.249274 + j -0.969559 4.259899 + j 1.010657 0.010625 + j 0.041098

11 6 3.465003 + j -0.270003 -3.457273 + j 0.209136 0.007730 + j -0.060866

8 7 -1.909893 + j -0.196732 1.911274 + j 0.129559 0.001381 + j -0.067173

9 8 0.132235 + j 0.116464 -0.131977 + j -0.521432 0.000258 + j -0.404968 39 9 7.617154 + j -1.902126 -7.557438 + j 2.142687 0.059717 + j 0.240561

11 10 -3.483660 + j -0.203064 3.488121 + j 0.171352 0.004461 + j -0.031712

13 10 -3.008372 + j -0.730489 3.011879 + j 0.688680 0.003508 + j -0.041809

14 13 -2.934129 + j -0.411463 2.941429 + j 0.307264 0.007300 + j -0.104199

15 14 -0.311115 + j -0.998556 0.312417 + j 0.627891 0.001303 + j -0.370665

16 15 2.896296 + j 0.430232 -2.888885 + j -0.531444 0.007411 + j -0.101212

17 16 -2.048841 + j 0.950740 2.052282 + j -1.049122 0.003441 + j -0.098383 19 16 4.542969 + j 0.681545 -4.511670 + j -0.625873 0.031300 + j 0.055672

21 16 3.324778 + j -0.302389 -3.316338 + j 0.177006 0.008440 + j -0.125383 24 16 0.410793 + j -0.811601 -0.410571 + j 0.744757 0.000222 + j -0.066844 18 17 -1.917560 + j 0.363855 1.920087 + j -0.475208 0.002527 + j -0.111353 27 17 -0.128621 + j 0.132648 0.128754 + j -0.475531 0.000133 + j -0.342884 22 21 6.093176 + j 1.070437 -6.064778 + j -0.847611 0.028398 + j 0.222826

23 22 -0.406149 + j -1.116063 0.406824 + j 0.928040 0.000675 + j -0.188024

24 23 -3.490793 + j -0.110399 3.516516 + j 0.138837 0.025723 + j 0.028438

26 25 -0.771398 + j -0.442881 0.773189 + j -0.111580 0.001791 + j -0.554460

27 26 -2.681379 + j -0.887648 2.691475 + j 0.731900 0.010096 + j -0.155748

28 26 1.416063 + j -0.565082 -1.408178 + j -0.210747 0.007885 + j -0.775830

29 26 1.921038 + j -0.679443 -1.901899 + j -0.248272 0.019138 + j -0.927715 29 28 3.491624 + j -0.395924 -3.476063 + j 0.289082 0.015561 + j -0.106842 4 0 0.000000 + j -1.074485 0.000000 + j 0.000000 0.000000 + j -1.074485

5 0 0.000000 + j -2.18259

6 0.000000 + j 0.000000 0.000000 + j -2.182596

11 12 0.018656 + j 0.473066 -0.018327 + j -0.464126 0.000329 + j 0.008940

13 12 0.066943 + j 0.423225 -0.066673 + j -0.415874 0.000270 + j 0.007351

30 2 7.897633 + j -0.731582 -7.897633 + j 1.860277 0.000000 + j 1.128695

31 6 7.506817 + j 1.371343 -7.506817 + j 0.109153 0.000000 + j 1.480496

32 10 12.260592 + j 5.296517 -12.260592 + j -2.064007 0.000000 + j 3.232509

34 20 5.080000 + j 1.826359 -5.054406 + j -1.314473 0.025594 + j 0.511886

33 19 -1.716763 + j 5.348910 1.736896 + j -4.940504 0.020133 + j 0.408405

35 22 6.500000 + j 2.637566 -6.500000 + j -1.998477 0.000000 + j 0.639089

36 23 1.402814 + j -0.195113 -1.401865 + j 0.246763 0.000949 + j 0.051650

37 25 9.586236 + j 0.419689 -9.533808 + j 1.607517 0.052428 + j 2.027206

38 29 -12.165903 + j 2.106593 12.280860 + j 0.135062 0.114957 + j 2.241655

20 19 -1.745594 + j 0.284473 1.747837 + j -0.240265 0.002242 + j 0.044208

结果分析：此程序的运行结果和试验程序给出的结果是一致的。说明程序无误，但在精确度上有微小差异，这主要是和导纳矩阵的精确度以及显示精度有关。

心得：本程序分模块进行，先是排序，再是求导纳阵，然后求雅阁比，再进行迭代运算，程序本身很简洁明了，运行的时候只需要在matlab里输入main就行了，然后打开BUS和line所在的.m文件，结果就会自动存在result文件中了，通过编写牛顿拉夫逊法matlab潮流计算程序复习了潮流计算的知识，也实现了计算机算法

附录：

实验源程序：

Main函数：

clear

[dfile,pathname]=uigetfile('*.m','Select Data File');

if pathname == 0

error(' you must select a valid data file')

else

lfile =length(dfile);

% strip off .m

eval(dfile(1:lfile-2));

end

[nb,mb]=size(bus);

[nl,ml]=size(line); % 计算bus和line矩阵的行数和列数[bus,line,nPQ,nPV,nodenum] = Num(bus,line); % 对节点重新排序的子程序

Y = y(bus,line) % 计算节点导纳矩阵的子程序

myf = fopen('Result.m','w');

fprintf(myf,'计算结果');

fclose(myf); % 在当前目录下生成“Result.m”文件，写入节点导纳矩阵

format long

EPS = 1.0e-10; % 设定误差精度

for t = 1:100 % 开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出

[dP,dQ] = dPQ(Y,bus,nPQ,nPV); % 计算功率偏差dP和dQ的子程序

J = Jac(bus,Y,nPQ); % 计算雅克比矩阵的子程序

UD = zeros(nPQ,nPQ);

for i = 1:nPQ

UD(i,i) = bus(i,2); % 生成电压对角矩阵 end

end

dAngU = J\[dP;dQ];

dAng = dAngU(1:nb-1,1); % 计算相角修正量

dU = UD*(dAngU(nb:nb+nPQ-1,1)); % 计算电压修正量

bus(1:nPQ,2) = bus(1:nPQ,2) - dU; % 修正电压

bus(1:nb-1,3) = bus(1:nb-1,3) - dAng; % 修正相角

if (max(abs(dU))

break

end% 判断是否满足精度误差，如满足则跳出，否则返回继续迭代

end

bus = PQ(bus,Y,nPQ,nPV); % 计算每个节点的有功和无功注入的子程序

[bus,line] = ReNum(bus,line,nodenum); % 对节点恢复编号的子程序

YtYm = YtYm(line); % 计算线路的等效Yt和Ym的子程序，以计算线路潮流

bus_res = bus_res(bus); % 计算节点数据结果的子程序

S_res = S_res(bus,line,YtYm); % 计算线路潮流的子程序

myf = fopen('Result.m','a');

fprintf(myf,'牛顿－拉夫逊法潮流计算结果节点计算结果：n节点节点电压节点相角（角度）节点注入功率\n');

for i = 1:nb

fprintf(myf,'%2.0f ',bus_res(i,1));

fprintf(myf,'%10.6f ',bus_res(i,2));

fprintf(myf,'%10.6f ',bus_res(i,3));

fprintf(myf,'%10.6f + j %10.6f\n',real(bus_res(i,4)),imag(bus_res(i,4))); end

fprintf(myf,'n线路计算结果：n节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)\n');

for i = 1:nl

fprintf(myf,'%2.0f ',S_res(i,1));

fprintf(myf,'%2.0f ',S_res(i,2));

fprintf(myf,'%10.6f + j %10.6f ',real(S_res(i,3)),imag(S_res(i,3)));

fprintf(myf,'%10.6f + j %10.6f ',real(S_res(i,4)),imag(S_res(i,4)));

fprintf(myf,'%10.6f + j%10.6f\n',real(S_res(i,5)),imag(S_res(i,5)));

end

fclose(myf); % 迭代结束后继续在“Result.m”写入节点计算结果和线路计算结果程序结束

"Num.m" 作用为对节点重排序，并修改相应的线路数据

function [bus,line,nPQ,nPV,nodenum] = Num(bus,line)

[nb,mb]=size(bus);

[nl,ml]=size(line);

nSW = 0; % number of swing bus counter

nPV = 0; % number of PV bus counter

nPQ = 0; % number of PQ bus counter

for i = 1:nb, % nb为总节点数

type= bus(i,6);

if type == 3,

nSW = nSW + 1; % increment swing bus counter

SW(nSW,:)=bus(i,:);

elseif type == 2,

nPV = nPV +1; % increment PV bus counter

PV(nPV,:)=bus(i,:);

else

nPQ = nPQ + 1; % increment PQ bus counter

PQ(nPQ,:)=bus(i,:);

end

end;

bus=[PQ;PV;SW];

newbus=[1:nb]';

nodenum=[newbus bus(:,1)];

bus(:,1)=newbus;

for i=1:nl

for j=1:2

for k=1:nb

if line(i,j)==nodenum(k,2)

line(i,j)=nodenum(k,1);

break

end

"y.m" 作用为计算节点导纳矩阵

function Y = y(bus,line)

[nb,mb]=size(bus);

[nl,ml]=size(line);

Y=zeros(nb,nb);

for k=1:nl

I=line(k,1); %读入线路参数

J=line(k,2);

Zt=line(k,3)+j*line(k,4);

Yt=1/Zt;

Ym=line(k,5)+j*line(k,6);

K=line(k,7);

if (K==0)&(J~=0) % 普通线路: K=0;

Y(I,I)=Y(I,I)+Yt+Ym;

Y(J,J)=Y(J,J)+Yt+Ym;

Y(I,J)=Y(I,J)-Yt;

Y(J,I)=Y(I,J);

end

if (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0; Y(I,I)=Y(I,I)+Ym;

end

if K>0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧

Y(I,I)=Y(I,I)+Yt+Ym;

Y(J,J)=Y(J,J)+Yt/K/K;

Y(I,J)=Y(I,J)-Yt/K;

Y(J,I)=Y(I,J);

end

if K<0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧

Y(I,I)=Y(I,I)+Yt+Ym;

Y(J,J)=Y(J,J)+K*K*Yt;

Y(I,J)=Y(I,J)+K*Yt;

Y(J,I)=Y(I,J);

end

"dPQ.m" 作用为计算功率偏差

function [dP,dQ] =dPQ(Y,bus,nPQ,nPV) % nPQ、nPV为相应节点个数

n = nPQ + nPV +1; % 总节点个数

dP = bus(1:n-1,4);

dQ = bus(1:nPQ,5); % 对dP和dQ赋初值 PV节点不需计算dQ 平衡节点不参与计算

for i = 1:n-1

for j = 1:n

dP(i,1) = dP(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))* sin(bus(i,3)-bus(j,3)));

if i

dQ(i,1) = dQ(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))* cos(bus(i,3)-bus(j,3)));

end

end% 利用循环计算求取dP和dQ

"Jac.m" 作用为计算雅克比矩阵

function J = Jac(bus,Y,nPQ)

[nb,mb]=size(bus);

H = zeros(nb-1,nb-1);

N = zeros(nb-1,nPQ);

K = zeros(nPQ,nb-1);

L = zeros(nPQ,nPQ); % 将雅克比矩阵分块，即：J = [H N; K L]，并初始化

Qi = zeros(nb-1,1);

Pi = zeros(nb-1,1);

for i = 1:nb-1

for j = 1:nb

Pi(i,1)=Pi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag( Y(i,j))*sin(bus(i,3)-bus(j,3)));

Qi(i,1)=Qi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag( Y(i,j))*cos(bus(i,3)-bus(j,3)));

end

end% 初始化并计算Qi和Pi

for i = 1:nb-1

for j = 1:nb-1

if i~=j

H(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))* cos(bus(i,3)-bus(j,3)));

else

H(i,j)=Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);

end% 分别计算H矩阵的对角及非对角元素

if j < nPQ+1

if i~=j

N(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))* sin(bus(i,3)-bus(j,3)));

else

N(i,j)=-Pi(i,1)-real(Y(i,j))*((bus(i,2))^2);

end

end% 分别计算N矩阵的对角及非对角元素

if i < nPQ+1

if i~=j

K(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*s in(bus(i,3)-bus(j,3)));

else

K(i,j)=-Pi(i,1)+real(Y(i,j))*((bus(i,2))^2);

end% 分别计算K矩阵的对角及非对角元素

if j < nPQ+1

if i~=j

L(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))* cos(bus(i,3)-bus(j,3)));

else

L(i,j)=-Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);

end

end% 分别计算L矩阵的对角及非对角元素

end

J = [H N; K L]; % 生成雅克比矩阵

"PQ.m" 作用为计算每个节点的功率注入

function bus = PQ(bus,Y,nPQ,nPV)

n = nPQ+nPV+1; % n为总节点数

for i = nPQ+1:n-1

for j = 1:n

bus(i,5)=bus(i,5)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-ima g(Y(i,j))*cos(bus(i,3)-bus(j,3)));

end

end% 利用公式计算PV节点的无功注入

for j =1:n

bus(n,4)=bus(n,4)+bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+ima g(Y(n,j))*sin(bus(n,3)-bus(j,3)));

bus(n,5)=bus(n,5)+bus(n,2)*bus(j,2)*(real(Y(n,j))*sin(bus(n,3)-bus(j,3))-ima g(Y(n,j))*cos(bus(n,3)-bus(j,3)));

end% 计算平衡节点的无功及有功注入

"ReNum.m" 作用为对节点和线路数据恢复编号

function [bus,line] = ReNum(bus,line,nodenum)

[nb,mb]=size(bus);

[nl,ml]=size(line);

bus_temp = zeros(nb,mb); % bus_temp矩阵用于临时存放bus矩阵的数据

k = 1;

for i = 1 :nb

for j = 1 : nb

if bus(j,1) == k

bus_temp(k,:) = bus(j,:);

k = k + 1;

end

end% 利用bus矩阵的首列编号重新对bus矩阵排序并存入bus_temp矩阵中

bus = bus_temp; % 重新赋值回bus，完成bus矩阵的重新编号

for i=1:nl

for j=1:2

for k=1:nb

if line(i,j)==nodenum(k,1)

line(i,j)=nodenum(k,2);

break

end

end% 恢复line的编号

"YtYm.m" 作用为计算线路的等效Yt和Ym，以计算线路潮流

function YtYm = YtYm(line)

[nl,ml]=size(line);

YtYm = zeros(nl,5); % 对YtYm矩阵赋初值0

YtYm(:,1:2) = line(:,1:2); % 矩阵前两列为线路两段节点编号，后三列分别为线路等效Yt，i侧的等效Ym，j侧的等效Ym

j = sqrt(-1);

for k=1:nl

I=line(k,1);

J=line(k,2);

Zt=line(k,3)+j*line(k,4);

if Zt~=0

Yt=1/Zt;

end

Ym=line(k,5)+j*line(k,6);

K=line(k,7);

if (K==0)&(J~=0) % 普通线路: K=0

YtYm(k,3) = Yt;

YtYm(k,4) = Ym;

YtYm(k,5) = Ym;

end

if (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0

YtYm(k,4) = Ym;

end

if K>0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧

YtYm(k,3) = Yt/K;

YtYm(k,4) = Ym+Yt*(K-1)/K;

YtYm(k,5) = Yt*(1-K)/K/K;

end

if K<0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧

YtYm(k,3) = -Yt*K;

YtYm(k,4) = Ym+Yt*(1+K);

YtYm(k,5) = Yt*(K^2+K);

end

"bus_res.m" 计算并返回节点数据结果

function bus_res = bus_res(bus)

[nb,mb]=size(bus);

bus_res = zeros(nb,4); % bus_res矩阵储存着节点计算结果

bus_res(:,1:2) = bus(:,1:2);

bus_res(:,3) = bus(:,3) *180 / pi; % 相角采用角度制

bus_res(:,4) = bus(:,4) + (sqrt(-1))*bus(:,5); % 注入功率

end

"S_res.m" 计算并返回线路潮流

function S_res = S_res(bus,line,YtYm)

[nl,ml]=size(line);

S_res = zeros(nl,5); % S_res矩阵储存着线路潮流计算结果

S_res(:,1:2) = line(:,1:2); % 前两列为节点编号

for k=1:nl

I = S_res(k,1);

J = S_res(k,2);

if (J~=0)&(I~=0)

S_res(k,3)=bus(I,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4)))-bus(I,2)*bus(J,2)*(c os(bus(I,3))+j*sin(bus(I,3)))*(conj(cos(bus(J,3))+j*sin(bus(J,3))))*conj(YtY m(k,3));

S_res(k,4)=bus(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))-bus(I,2)*bus(J,2)*(c onj(cos(bus(I,3))+j*sin(bus(I,3))))*(cos(bus(J,3))+j*sin(bus(J,3)))*conj(YtY m(k,3));

S_res(k,5)=S_res(k,3) + S_res(k,4); % 利用公式计算非接地支路的潮流else if(J==0 )

S_res(k,3)=bus(I,2)^2*conj(YtYm(k,4));

S_res(k,5)=S_res(k,3)+S_res(k,4);

else

S_res(k,4)=bus(J,2)^2*conj(YtYm(k,5));

S_res(k,5)=S_res(k,3)+S_res(k,4); % 利用公式计算接地支路的潮流end

end

作业

y=pifun(10000)

y =

3.1415

>> clear all >> tic;

>> y=pifun(100000); >> toc;

Elapsed time is 0.014727 seconds.

function y= pifun(n) y=0;

for i=1:n

y=y+(-1)^(i-1)*1/(2*i-1); end

y=y*4; end

>> a=0:pi/12:2*pi; >> x=5*cos(a); >> y=3*sin(a); >> plot(x,y)

-5

-4-3-2-1012345

-3-2

-1

牛顿插值法原理及应用

牛顿插值法插值法是利用函数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 void main() { float x[11],y[11][11],xx,temp,newton; int i,j,n; printf("Newton插值:\n请输入要运算的值:x="); scanf("%f",&xx); printf("请输入插值的次数(n<11):n="); scanf("%d",&n); printf("请输入%d组值:\n",n+1); for(i=0;i

牛顿拉夫逊法潮流计算

摘要本文，首先简单介绍了基于在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

基于极坐标的牛顿拉夫逊潮流计算

前言电力系统潮流计算是研究电力系统稳态运行情况的一种计算,它根据给定的运行件及系统接线情况确定整个电力系统各部分的运行状态。在电力系统规划设计和现有电力系统运行方式的研究中,都需要利用潮流计算来定量分析、比较供电方案或运行方式的合理性、可靠性和经济性。本次课程设计任务是闭式网络的潮流计算，用到的方法为牛顿拉夫逊极坐标法潮流计算。牛顿法是数学中解决非线性方程式的典型方法，有较好的收敛性。解决电力系统潮流计算问题是以导纳距阵为基础的，因此，只要在迭代过程中尽可能保持方程式系数距阵的稀疏性，就可以大大提高牛顿法潮流程序的放率。自从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

matlab实现数值分析报告插值及积分

Matlab实现数值分析插值及积分摘要：数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。其中Aitken插值计算的结果图如下：对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下：

问题重述问题一：已知列表函数表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。问题解决问题一：插值方法对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。一、拉格朗日插值法：拉格朗日插值多项式如下：首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得：

PQ分解法潮流计算实验

xxxx实验报告学生姓名：学号：专业班级：实验类型：□验证□综合■设计□创新实验日期：实验成绩：一、实验项目名称 P-Q分解法潮流计算实验二、实验目的与要求：目的：电力系统分析的潮流计算是电力系统分析的一个重要的部分。通过对电力系统潮流分布的分析和计算，可进一步对系统运行的安全性，经济性进行分析、评估，提出改进措施。电力系统潮流的计算和分析是电力系统运行和规划工作的基础。潮流计算是指对电力系统正常运行状况的分析和计算。通常需要已知系统参数和条件，给定一些初始条件，从而计算出系统运行的电压和功率等；潮流计算方法很多：高斯-塞德尔法、牛顿-拉夫逊法、P-Q分解法、直流潮流法，以及由高斯-塞德尔法、牛顿-拉夫逊法演变的各种潮流计算方法。本实验采用P-Q分解法进行电力系统分析的潮流计算程序的编制与调试，获得电力系统中各节点电压，为进一步进行电力系统分析作准备。通过实验教学加深学生对电力系统潮流计算原理的理解和计算，初步学会运用计算机知识解决电力系统的问题，掌握潮流计算的过程及其特点。熟悉各种常用应用软件，熟悉硬件设备的使用方法，加强编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。要求：编制调试电力系统潮流计算的计算机程序。程序要求根据已知的电力网的数学模型（节点导纳矩阵）及各节点参数，完成该电力系统的潮流计算，要求计算出节点电压、功率等参数。三、主要仪器设备及耗材每组计算机1台、相关计算软件1套四、实验内容： 1.理论分析： P-Q分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思想是：把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，把有功功率和无功功率迭代分开来进行。牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系统时，

基于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

电力系统课程设计-牛顿拉夫逊法潮流计算

课程设计说明书题目电力系统分析系(部) 专业(班级) 姓名学号指导教师起止日期

电力系统分析课程设计任务书系（部）：专业：指导教师：

目录一、潮流计算基本原理 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牛顿插值法例题与程序

题目一：多项式插值某气象观测站在8：00（AM ）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton ）逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度（x=10）。二、数学原理假设有n+1个不同的节点及函数在节点上的值（x 0，y 0），……（x n ，y n ），插值多项式有如下形式：）（））（（）（）（）（n 10n 102010n x -x )(x -x x -x x P x x x x x x -??-+??+-++=αααα （1）其中系数i α（i=0,1,2……n ）为特定系数，可由插值样条i i n y x P =）（（i=0,1,2……n ）确定。根据均差的定义，把x 看成[a,b]上的一点,可得 f(x)= f （0x ）+f[10x x ，]（0x -x ） f[x, 0x ]= f[10x x ，]+f[x,10x x ，] （1x -x ） …… f[x, 0x ，…x 1-n ]= f[x, 0x ，…x n ]+ f[x, 0x ，…x n ](x-x n ) 综合以上式子，把后一式代入前一式，可得到： f(x)= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n )+ f[x, 0x ，…x n ,x ]）（x 1n +ω= N n （x ）+）（x n R 其中

N n （x ）= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n ) （2））（x n R = f(x)- N n （x ）= f[x, 0x ， (x) n ,x ]）（x 1n +ω （3））（x 1n +ω=（0x -x ）…(x-x n ) Newton 插值的系数i α（i=0,1,2……n ）可以用差商表示。一般有 f k =α[k 10x x x ??，] （k=0,1,2，……，n ）（4）把（4）代入（1）得到满足插值条件N ）（）（i i n x f x =（i=0,1,2，……n ）的n 次Newton 插值多项式 N n （x ）=f （0x ）+f[10x x ，]（1x -x ）+f[210x x x ，，]（1x -x ）（2x -x ）+……+f[n 10x x x ??，]（1x -x ）（2x -x ）…（1-n x -x ）. 其中插值余项为：）（）！（）（）（）（）（x 1n f x N -x f x R 1n 1 n n +++==ωξ ξ介于k 10x x x ??，之间。三、程序设计 function [y,A,C,L]=newdscg(X,Y,x,M) % y 为对应x 的值，A 为差商表，C 为多项式系数，L 为多项式 % X 为给定节点，Y 为节点值，x 为待求节点 n=length(X); m=length(x); % n 为X 的长度 for t=1:m

电力系统稳态分析牛顿拉夫逊法

0引言潮流是配电网络分析的基础，用于电网调度、运行分析、操作模拟和设计规划，同时也是电压优化和网络接线变化所要参考的内容。潮流计算通过数值仿真的方法把电力系统的详细运行情况呈现给工作人员，从而便于研究系统在给定条件下的稳态运行特点。随着市场经济的发展，经济利益是企业十分看重的，而线损却是现阶段阻碍企业提高效益的一大因素。及时、准确的潮流计算结果，可以给出配电网的潮流分布、理论线损及其在网络中的分布，从而为配电网的安全经济运行提供参考。从数学的角度来看，牛顿-拉夫逊法能有效进行非线性代数方程组的计算且具有二次收敛的特点，具有收敛快、精度高的特点，在输电网中得到广泛应用。随着现代计算机技术的发展，利用编程和相关软件，可以更好、更快地实现配电网功能，本文就是结合牛顿-拉夫逊法的基本原理，利用C++程序进行潮流计算，计算结果表明该方法具有良好的收敛性、可靠性及正确性。 1牛顿-拉夫逊法基本介绍 1.1潮流方程对于N 个节点的电力网络（地作为参考节点不包括在内），如果网络结构和元件参数已知，则网络方程可表示为： =&& YV I （1-1）式中，Y 为N*N 阶节点导纳矩阵；&V 为N*1维节点电压列向量；&I 为N*1维节点注入电流列向量。如果不计网络元件的非线性，也不考虑移相变压器，则Y 为对称矩阵。电力系统计算中，给定的运行变量是节点注入功率，而不是节点注入电流，这两者之间有如下关系： ??=&&&EI S （1-2）式中，&S 为节点的注入复功率，是N*1维列矢量；? &S 为&S 的共轭；

??i diag ??=???? &&E V 是由节点电压的共轭组成的N*N 阶对角线矩阵。由（1-1）和（1-2），可得： ??=&&&S EYV 上式就是潮流方程的复数形式，是N 维的非线性复数代数方程组。将其展开，有： ?i i i ij j j i P jQ V Y V ∈-=∑&& j=1,2，….，N （1-3）式中，j i ∈表示所有和i 相连的节点j ，包括j i =。将节点电压用极坐标表示，即令i i i V V θ=∠&，代入式（1-3）中则有： ()i i i i ij ij j j j i P jQ V G jB V θθ∈-=∠-+∠∑ ()()cos sin i j ij ij ij ij j i V V G jB j θθ∈=+-∑ 故有： () ()cos sin sin cos i i j ij ij ij ij j i i i j ij ij ij ij j i P V V G B Q V V G B θθθθ∈∈?=+?? =-?? ∑∑ i=1,2,…,N （1-4）式（1-4）是用极坐标表示的潮流方程。而节点功率误差： (cos sin )θθ∈?=-+∑SP i i i j ij ij ij ij j i P P V V G B (1-5) (cos sin )θθ∈?=--∑SP i i i j ij ij ij ij j i Q Q V V G B (1-6) 式中：SP i P ，SP i Q 为节点i 给定的有功功率及无功功率。 1.2牛顿-拉夫逊法基本原理 1.2.1牛拉法的一般描述牛拉法是把非线性方程式的求解过程变成反复对相应的线性方程式的求解过程，即非线性问题通过线性化逐步近似，这就是牛拉法的核心。下面以非线性方程式的求解过程来进行说明。设电力网络的节点功率方程一般形式如下：

牛顿法潮流计算综述

潮流例题：根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择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.设计流程图

牛顿插值MATLAB算法

MATLAB程序设计期中作业 ——编程实现牛顿插值成员:刘川（P091712797）签名＿＿＿＿＿汤意（P091712817）签名＿＿＿＿＿王功贺（P091712799）签名＿＿＿＿＿班级:2009信息与计算科学学院:数学与计算机科学学院日期:2012年05月02日

牛顿插值的算法描述及程序实现一：问题说明在我们的实际应用中，通常需要解决这样的问题，通过一些已知的点及其对应的值，去估算另外一些点的值，这些数据之间近似服从一定的规律，于是，这就引入了插值法的思想。插值法是利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。利用插值基函数很容易得到拉格朗日插值多项式，公式结构紧凑，在理论分析中甚为方便，但当插值节点增减时全部插值基函数均要随之变化，整个公式也将发生变化，这在实际计算中是很不方便的，为了克服这一缺点，提出了牛顿插值。二：算法分析 newton 插值多项式的表达式如下： 010011()()()()()n n n N x c c x x c x x x x x x -=+-+???+--???- 其中每一项的系数c i 的表达式如下： 12011010 [,,,][,,,] [,,,]i i i i i f x x x f x x x c f x x x x x -???-???=???= - 即为f (x)在点01,,,i x x x ???处的i 阶差商，（[]()i i f x f x =，1,2,,i n = ），由差商01[,,,]i f x x x ???的性质可知： () 010 1 [,,,]()i i i j j k j k k j f x x x f x x x ==≠???=-∑∏ 牛顿插值的程序实现方法：第一步：计算[][][][]001012012,,,,,,,n f x f x x f x x x f x x x x 、、、、。第二步：计算牛顿插值多项式中01[,,,]i f x x x ???011()()()i x x x x x x ---???-，1,2,,i n = ，得到n 个多项式。

c语言编写的牛顿拉夫逊法解潮流程序

c语言编写的牛顿拉夫逊法解潮流程序闲来无事，最近把牛拉法用c语言重写一遍，和matlab相比，c语言编写潮流程序最大的难点在于矩阵求逆，我使用的求逆方法是初等行变换法，程序段如下： #include #define N 3 void main() { int i,j,k; float t; float Jacob[N][N]={{1,2,2},{1,3,4},{2,3,4}};//欲进行求逆的矩阵 float inv_J[N][N];//逆矩阵存储于此 //初始化inv_J[N][N] for(i=0;i

} //输出逆矩阵 for(i=0;i #include #define N 4 //节点数 #define n_PQ 2 //PQ节点数 #define n_PV 1 //PV节点数 #define n_br 5 //串联支路数 void main() { void disp_matrix(float *disp_p,int disp_m,int disp_n); //矩阵显示函数 float Us[2*N]={1.0,0,1.0,0,1.05,0,1.05,0}; //电压初值 float Ps[N]={0,-0.5,0.2}; //有功初值 float Qs[N]={0,-0.3}; //无功初值float G[N][N],B[N][N]; //各几点电导电纳 struct //阻抗参数 { int nl; //左节点 int nr; //右节点 float R; //串联电阻值 float X; //串联电抗值 float Bl; //左节点并联电导 float Br; //右节点并联电纳 }ydata[n_br]={ {1,2,0,0.1880,-0.6815,0.6040}, {1,3,0.1302,0.2479,0.0129,0.0129}, {1,4,0.1736,0.3306,0.0172,0.0172}, {3,4,0.2603,0.4959,0.0259,0.0259}, {2,2,0,0.05,0,0} };

牛顿插值法的MATLAB综合程序

6.3.5 牛顿插值法的MATLAB 综合程序求牛顿插值多项式、差商、插值及其误差估计的MATLAB 主程序 function [y,R,A,C,L]=newdscg(X,Y,x,M) n=length(X); m=length(x); for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y'; s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end q1=abs(q1*(z-X(j-1)));c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n))); for k=(n-1):-1:1 C=conv(C,poly(X(k))); d=length(C);C(d)=C(d)+A(k,k); end y(k)= polyval(C, z); end R=M*q1/c1;L(k,:)=poly2sym(C); 例6.3.6 给出节点数据00.27)00.4(=-f ，00.1)00.0(=f ，00.2)00.1(=f ，00.17)00.2(=f ，作三阶牛顿插值多项式，计算)345.2(-f ，并估计其误差. 解首先将名为newdscg.m 的程序保存为M 文件，然后在MATLAB 工作窗口输入程序 >> syms M,X=[-4,0,1,2]; Y =[27,1,2,17]; x=-2.345; [y,R,A,C,P]=newdscg(X,Y,x,M) 运行后输出插值y )345.2(-≈f 及其误差限公式R ，三阶牛顿插值多项式P 及其系数向量C ，差商的矩阵A 如下 y = 22.3211 R = 65133/562949953421312*M （即R =2.3503*M ） A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167 C = 0.9167 4.2500 -4.1667 1.0000 P = 11/12*x^3+17/4*x^2-25/6*x+1

直角坐标牛顿-拉夫逊法潮流计算matlab程序(仅供参考)

%该程序仅针对《电力系统分析下》何仰赞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 =。由此可知节点导纳矩阵是一个稀疏的对称矩阵。