matlab电力系统潮流计算

1、华中科技大学信息工程学院课程设计报告书题目:电力系统潮流计算专 业：电气工程及其自动化班 级：学 号：学生：指导教师：学号»不无”,比尔、电气工程及其自学生专业班级）.、设计题目电力系统潮流计算设计技术参数1 .变压器的阻抗，励磁损耗。2 .线路的阻抗，无功损耗。3 .变压器及线路的功率损耗4 .变压器及线路的电压降落5 .各节点的功率和电压设 计 要 求1 .计算过程包括理论计算和MATLAB仿真的大部分。2 .报告里面的理论分析（手工计算）过程必须要后详细推导过程，再把仿真过程 放到相应的部分中去。工 作 量共24页；4000字工 作 计 划第15周：对系统模型分析和资料收集，找

2、到计算的思路，总结计算方法；第16周：对整个模型进行详细计算，并完成word文档。参 考 资 料1何仰赞，温增银电力系统分析（第三版）M.华中科技大学，20022同娟.MATLAB在电路分析中的应用.电气电子教学学报.20023王守相，玉田 电力系统潮流计算研究现状 -电力技术1996年05期指导教师签字2015年11月12日信息工程学院课程设计成绩评定表学生： 学号：专业（班级）：n课程设计题目：电力系统潮流计算指导教师评语:成绩：指导教师：摘要电力系统稳态分析包括潮流计算和静态安全分析。 本文主要运用的事潮流计 算，潮流计算是电力网络设计与运行中最基本的运算， 对电力网络的各种设计方 案及

3、各种运行方式进行潮流计算， 可以得到各种电网各节点的电压，并求得网络 的潮流及网络中的各元件的电力损耗， 进而求得电能损耗。本位就是运用潮流计 算具体分析，并有MATLAB仿真。关键词：电力系统潮流计算 MATLAB仿真AbstractElectric power system steady flow calculation and analysis of the static safety analysis. This paper, by means of the calculation, flow calculation is the trend of the power network

4、design and operation of the most basic operations of electric power network, various design scheme and the operation ways to tide computation, can get all kinds of each node of the power grid voltage and seek the trend of the network and the network of the components of the power loss, and getting e

5、lectric power. The standard is to use the power flow calculation and analysis, the specific have MATLAB simulation.Key words: Power system; Flow calculation; MATLAB simulation1任务提出与方案论证 22总体设计32.1 潮流计算等值电路32.2 建立电力系统模型32.3 模型的调11t与运行33详细设计43.1 计算前提43.2 手工计算 74设计图及源程序114.1 MATLAB 仿真114.2 潮流计算源程序115总结

6、19参考文献201任务提由与方案论证潮流计算是在给定电力系统网络结构、参数和决定系统运行状态的边界条件的情况 下确定系统稳态运行状态的一种基本方法，是电力系统规划和运营中不可缺少的一个重要 组成部分。可以说，它是电力系统分析中最基本、最重要的计算，是系统安全、经济分析 和实时控制与调度的基础。常规潮流计算的任务是根据给定的运行条件和网路结构确定 整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损 耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。在电力系统运行方式和 规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可 靠性和经济性。同

7、时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮 流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。 在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的 实时监控中，则采用在线潮流计算。是电力系统研究人员长期研究的一个课题。它既是 对电力系统规划设计和运行方式的合理性、可靠性及经济性进行定量分析的依据，又是电力系统静态和暂态稳定计算的基础。潮流计算经历了一个由手工到应用数字电子计算机的发展过程，现在的潮流算法都 以计算机的应用为前提用计算机进行潮流计算主要步骤在于编制计算机程序，这是一项 非常复杂的工作。对系统进行潮流分析，本文利

8、用MATLAB中的Simpower Systems工具箱设计电力系统，在 simulink环境下，不仅可以仿真系统的动态过程，还可以对系统 进行稳态潮流分析。2总体设计SimPowerSystems使用Simulink环境，可以将该系统中的发电机、变压器，线路等 模型联结起来，形成电力系统仿真模拟图。在加人测量模块，并对各元件的参数进行设 置后，用measurement和sink中的仪器可以观察各元件的电压、电流、功率的大小。2.1 潮流计算等值电路10MW12CMW+10kVP0 15.7kWYN,d11ps 73kW16MWVA' ： I。 0.5丫 Vs% 10.5130km32

9、km/25MW /25MW110<VYN,d112 20MWVA20MWP0Ps10%Vs%4 15MW18.6 kW89kWO530MW18.710kVP0 15.7kWps 73<WI0% 0,5 gVs% 10.5YN, d114 63MWVA-w-<D>-po 44<WPs 12kWIo% 0.3570km110<V15MWXd 0.136X2 0.16 x 0.073 cos n 0.84 63MWXd 0.134X2 0.161X 0.06cos n 0.85P0Ps10%Vs%10%110VYN,d112 16MWVA63M10kVRRI0%V

10、s%G G：3 12MW12kW 0.351。5、G.1 50MWx 0.136X2 0.161X0 0.075 cos n 0.8Xd 0.128X2 0.154x 0.054cos n 0.8511kW 50k 0.55 105132kW63kW0.55Vs（1 2）% 10.580MW Vs（2 3）% 6.5 5MW+ Vs（1 3）% 17535M0 02 25MW10kVYN,d112 10MWVA35kV10kVXd 0.128X2 0.157X0 0.0591cos n 0.8YN,Y,d11 2 10MVA35kV2.2 建立电力系统模型在Simulink中按照电力系统原型选

11、择元件进行建模。所建立的模型和建立的方法在详细设计中详述。在电力系统模型的建立工程中主要涉及到的是：元器件的选择及其参数的设置；发电机选型；变 压器选择；线路的选择；负荷模型的选择；母线选择。2.3 模型的调试与运行Simulink环境下进行仿真运行。运行的具体结果和建立系统模型，并设置好参数以后，就可以在 分析也在详细设计中详述。3.1计算前提首先是发电机的参数计算，先对发电机G1 :PQi发电机G2:P2Q2发电机G3:P3Q3发电机G4:P4Q4发电机G5:P5Q53详细设计5个发电厂简化为5台发电机来计算。4 15 60MW60 tan(arccos0.8)4 63 252MW252

12、tan(arccos0.85)3 12 36MW36 tan(arccos0.8)1 50 50MW50 tan(arccos0.85)2 25 50MW50其次是变电站的参数计算，我们还是对变压器T1:Rt1ps VN变压器T2:（双并联）45MVar156MVar27MVar31MVartan(arccos0.8) 37.5MVar7个变电站简化为7台变压器来计算。103Xt12Vs%VN2Sn10. Io%p0j而RT2XT2Ps V；S2SNVs% vN103S02Sn10/. Io%(po j10073 11023 2 (16 10 )_210.5 110216 1033103.45

13、010 79.406SNSn)(0.0157 j0.0800)MVA_289 1102(20 103)210.5 110220 103103 1.34610 31.7625(0.0372j0.2000)MVA变压器T3:（四并联）变压器T4:（双并联）RT3Xt3S03Ps V；sNVs% Vn2Sn10310Sn)2 121 1102 (63 1 03)2_210.5 1102363 10(0.1760103 0.09210 5.042j0.8820)MVA_1 _RT4Rt1 1.725021XT4 Xt1 39.7030 2S04 2 S01 (0.0314 j0.1600)MVA变压器

14、T5:Rt5 4Rt3 0.3680XT5 4XT3 20.168c 1 c&513 (0.0440 j0.2205)MVA4变压器T6:（两个三绕组变压器并联）RT6 1RT6 2RT6 32163 352二 Z _ 3 22 (10 10 )1030.3861c,c,Vs1%- Vs(1 2) % Vs(1 3)% Vs(2 3)%10.75c,1c,c,c,Vs2%-Vs(12)%Vs(23)%Vs(13)%0.25、，1Vs3%- Vs(i 3)% Vs(2 3)% Vs(i 2)%6.75XT6 1Vsi%V；XT6 2SNVs2%VN210XT6 3SN_ ,2Vs3%Vn

15、10SN106.5840.1534.134Io%S06 2 ( P06 j 10) (0.0264 j0.1100)MVA 100变压器T7:(双并联)线路线路线路线路RT7XT7S37再次是传输线参数计算,根据教材查得L1:L2:L3:（双回路）L4:线路L5:（双回路）r。0.21ps VN 1q3 * 102Vs% VNSN10/ c ； I。p。woSN)50 352 (10 103)2 10.5 35210 103(0.02205条传输线的具体计算如下。/km xq0.4/ kmbg2.8RL1X L1Bl1Ql1RL2Xl2Bl2Xl3BL3r。110.21408.4Xqber。X

16、qb。11110.42.812 BL1VN2l2l21240160.210.42.8QL2QL3r。10 640 1.120.6776MVar13013010 627.352103 0.30610 6.431j0.1100)MVA106S/km10 4S130 3.64 10 4S2 BL2VN22.2022MVarl3Xqbol3l30.21 700.4 702.8 1067.351470 3.92 10 4S2.3716MVarRL4r。l40.216012.6Xl4XqI40.46024BL4bQl42.810 6一一一 4 一60 1.68 10 S12一Bl3Vn2QL41.0164

17、MVar1 Bl4V。211RL5-r015- 0.21 202.122、，1.1 c ccXL5x0 150.4 20 422BL5 2 b0 15 2 2.8 10 6 20 1.12 10 4S12QL5 BL3VN20.0686MVar23.2手工计算FLR1:ST1P2LRt1jX)102(3.4501102j74.406) (0.0285 j0.6562)MVASa 10MWST1P2 Q2Sl1 (Rl1VNP2 Q2ST 272( RT2SbScSG1 20ST2 60 j45 20 0.4032 j9.5156 (39.5968 j35.4844)MVA FLR2Sb Sa

18、25 jQL1SL1 (4.4826 j35.9144)MVA01 j Qli (10.0442 j0.1142)MVA10.04422 0.11422jX li)(8.4 j16) (0.070 j0.1334)MVA402 452jXT?)1102 (1.346 j31.7625) (0.4032 j9.5156)MVAST3P2 Q2(FT3 jXT3)Vn2522 15622-(0.092 j5.042) (0.6679 j36.6024)MVA1103：,，Sc (4.4931 j34.1048) MVA2_ 2P2 Q2SL2- (RL2Vn2Sd Sg2 Sc 1202_ 2P2

19、 Q2ST 42( Rt4VnjXL2)ST3jXT4)一，2 一 一 一4.493134.104811022-(27.3 j52) (2.67 j5.0854)MVAFLRS0362jQ L2SL227211022-(1.725(132.9792 j149.229)MVAj39.703) (0.1091 j2.5101)MVA,，Sd (133.5955 j149.9956) MVAc P2 Q2133.59552 149.99562SL3 k(R3 'J -"""- a?5 j14) (24.51 j46.682)MVA_ '_ _ _ _ _

20、 _ _ _ _ _ _ _SeSG3Sd30258T4S04jQL3SL3(89.945j130.0151)MVAFLR4:S P2 Q2Sr5 ySe (92.7481502 312(Rt5 jXT5)1102 (0.368 j20.168)j133.9937)MVASL4P2 Q2Vn2(RL4jXL4)92.74872 133.99372_21102(12.6(0.1052 j5.7687)MVAj24) (27.654 j52.674)MVASfSg4Se 80SrsS05 jQL4SL4(34.9449j107.3469)MVAFLR5:ST7152r (0.306352j6.43

21、1) (0.0562j1.1812)MVASh 15St7j QL5 (15.0782 j0.3422)MVASL515.07822 0.34222SShSt 6352Sl5S06152 37.52352(2.1 j4) (0.3899 j0.743)MVAj QL5(0.3865 (20.4945 j1.1266)MVAj4.34)(0.514j5.7793)MVAST 6ST6 1SgSf计算每一个FLR1:FLR2:功率分布:SL220.65052 0.5451235226.3362 98.73692352ST6 1 SG5ST6(0.386(0.3862Sr 6FLR的功率分布和电压分

22、布计算如下:*ZT3*ZL2j0.153)j6.584)(0.1345 j0.0533)MVA(3.2905 j56.1256)MVA(25.5114 j194.12)MVAVT2PR QXVnVb 115 Vt2PR QXVbVli*ZT3Sd40 1.346 45 31.7625 12.8970kV115102.1030kV10.0442 8.4 0.1442 160.8489kV102.1030101.2541kV(Vb Vn)V*NZL2 ZT3(0.092 j5.042) (132.9792 j149.229) 1418.6727.392 j57.042(4.8812 j13.809

23、7)MVAST3*ZL2 s* j(Vb Vn)* *ZL2 ZT3 ZL2 ZT3 (108.687 j122.62)MVA(27.3 j52) (132.9792 j149.229) 1418.6727.392 j57.042电压分布:Sci Sl2Sl2 (4.8812 j13.8097) (2.67 j5.0854) (7.5512 j8.7243)MVA7.5512 27.3 8.7243 52Vl2102.1030Vd VbVl2 102.103 ( 2.424)2.424kV104.527kVFLR3:功率分布:ZT4 s(VG3 Vd)-3* J *ZL3 ZT4 ZL3 ZT

24、4 (59.444 j16.846)MVA *ZL3 s(Vb Vn)-4*De*ZL3 ZT4ZL3 ZT4(31.811 j60.1256)MVA(1.725 j39.703) (89.945 j130.0151) 1037.9279.075 j53.73(7.35 j14) (89.945 j130.0151) 1037.9279.075 j53.73电压分布：01SL3SL3 (59.444 j19.846) (24.51 j46.682) (83.954 j26.836)MVAVl383.954 7.35 26.836 14105.56439.404kVFLR4:Ve VdVL3 9

25、6.16kV功率分布:ZT5 S(VG3 Vd).4*f *ZL4 ZT5 ZL4 ZT5 (20.843 j19.689)MVA *ZL4 S(VG3 Vd)-4*f *ZL4 ZT5 ZL4 ZT5 (1.398 j44.389)MVA(0.368 j20.168) (34.9449 j107.3469) 1037.92712.968 j44.168电压分布：01SL3VL4VN 二(12.6 j24) (34.9449 j107.3469) 1037.92712.968 j 44.168SL3 (59.444 j16.846) (24.51 j46.682) (83.954 j63.52

26、8)MVA83.954 12.6 63.528 24105.564324.464kVVe VdVL3 81.10kVFLR5:这里我们先将f点和发电机G5当做电源，经过ZT61和Zt63构成两端供电网络以 g点作为运算负荷进行计算。ST6(0.386 j4.134) (20.2656 j70.9293) (22.0938 37) 35(3.900 j25i75)MVAST4电压分布：(0.386 j6.584) (20.2656 j70.9293) (22.0938 37) 350.772 j10.718(16.5061 j91.7905)MVAST631ST63ST63 (16.6421 j

27、97.5698)MVAVT6316.6421 0.386 97.5698 4.134 10.9186kV37Vg 37VT63 26.0814V20.2656 0.386 70.9293 ( 0.153)VT62"26.08140.1162kVVi VgVT62 26.197620.4945 2.1 1.1266 4Vl5 26.1976Vh ViVL5 24.38261.815kV4设计图及源程序4.1MATLAB 仿真相关的原始数据输入格式如下：1、B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写。2、对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧

28、节点编号，将变压器的 串联阻抗置于低压侧处理，第三列为支路的串列阻抗参数，第四列为支路的对地导纳参数，第五烈为 含变压器支路的变压器的变比，第六列为变压器是否是否含有变压器的参数，其中1”为含有变压器，0”为不含有变压器。3、B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中1”为平衡节点，2”为PQ节点，3”为PV节点参数。4、X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地参数。4.2潮流计算源程序close allclear alln=input('请输入节点数：n=');n1=

29、input('请输入支路数：n1=');isb=input('请输入平衡节点号：isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数旧1=');B2=input('请输入节点参数旧2=');X=input('节点号和对地参数：X=');Y=zeros(n);Times=1;%置迭代次数为初始值I建节点导纳矩阵for i=1:n1if B1(i,6)=0%含变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y

30、(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else喻有变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);Y(q,q户丫(q,q)+1/(B1(i,5)A2*B1(i,3);endendYOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%|各 OrgS、DetaS 初始化施|J建OrgS ,用于

31、存储初始功率参数h=0;j=0;for i=1:n%寸PQ节点的处理if i=isb&B2(i,6)=2h=h+1;forj=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*i mag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j

32、)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B 2(j,3)+imag(Y(i,j)*real(B2(j,3);endendendfor i=1:n%寸PV节点的处理，注意这时不可再将h初始化为0if i=isb&B2(i,6)=3h=h+1;forj=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*i mag(B2(j,3)+imag(Y(i,j

33、)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B 2(j,3)+imag(Y(i,j)*real(B2(j,3);endendendOrgS施建PVU用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)=3t=t+1;PVU(t,1)=B2(i,3);endendPVU施ij建DetaS,用于存储有功功率、无功功率和电

34、压幅值的不平衡量h=0;for i=1:n%寸PQ节点的处理if i=isb&B2(i,6)=2h=h+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1);endendt=0;for i=1:n%寸PV节点的处理，注意这时不可再将h初始化为0if i=isb&B2(i,6)=3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1)A2+imag(PVU(t,1)A

35、2-real(B2(i,3)A2-imag(B2(i,3)A2; endendDetaS施建I,用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:nif i=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3);endendI%J建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n%寸PQ节点的处理if B2(i,6)=2h=h+1;forj=1:nif j=isbk=k+1;if i=j%对角元素的处理Jacbi(2*h-1,2*k-1)=-ima

36、g(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1);Jacbi(2*h-1,2*k户real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1);else%曰对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i

37、,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k=(n-i) %将用于循环的指针置于初始值，以确保雅可比矩阵换行 k=0;endendendendendk=0;for i=1:n%寸PV节点的处理if B2(i,6)=3h=h+1;forj=1:nif j=isbk=k+1;if i=j%寸角元素的处理Jacbi(2*h-1,

38、2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=2*imag(B2(i,3);Jacbi(2*h,2*k)=2*real(B2(i,3);else%曰对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3);Jacbi(

39、2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k=(n-1)%各用于循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendJacbi%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%慎正节点电压j=0;for i=1:n%寸PQ节点处理if B2(i,6)=2j=j+1；B2(i,3)=B2(i,3)+DetaU(2*j,1

40、)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n%对PV节点的处理if B2(i,6)=3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB24始循环 *while abs(max(DetaU)>prOrgS=zeros(2*n-2,1);%!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0; j=0; for i=1:nif i=isb&B2(i,6)=2h=h+1;forj=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+

41、real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*i mag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B 2(j,3)+imag(Y(i,j)*real(B2(j,3);endendendfor i=1:n

42、if i=isb&B2(i,6)=3h=h+1;forj=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*i mag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(

43、Y(i,j)*imag(B20,3)+imag(Y(i,j)*real(B2(j,3);endendendOrgS%建 DetaSh=0;for i=1:nif i=isb&B2(i,6)=2h=h+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1);endendt=0;for i=1:nif i=isb&B2(i,6)=3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=real(P

44、VU(t,1)A2+imag(PVU(t,1)A2-real(B2(i,3)A2-imag(B2(i,3)A2; endendDetaS屈建Ii=zeros(n-1,1);h=0;for i=1:nif i=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3);endendI%J建 JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)=2h=h+1;forj=1:nif j=isbk=k+1;if i=jJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*r

45、eal(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1);Jacbi(2*h-1,2*k户real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1);Jacbi(2*h,2*k户Jacbi(2*h-1,2*k-1)-2*imag(I(h,1);elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k=(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)=3h=h+1;forj=1:nif j=isbk=k+1;if i=jJacbi(2*h-1,2*k-