电力系统分析课程设计文档

电力系统的潮流计算和短路故障的计算机算法程序设计

信息工程学院课程设计任务书学生姓名学号成绩设计题目电力系统的潮流计算和短路故障的计算机算法程序设计1.选择合适的计算机编程语言，在此选用matlab的m语言；2.理解牛顿拉夫逊算法计算方法和具体计算流程，并在计算机上编程和调试，和教材上的结果进行比较，最终得到正确结果；3.建立电力系统计算的相关数学模型，就是考虑影响问题的主要因素，设计内容而忽略一些次要因素，使数学模型既能正确地反映实际问题，又使计算不过于复杂，就是建立用于描述电力系统相应计算的有关参数间的相互关系的数学方程式；4.会做出三相对称短路的等效电路图以及不对称短路的正序，负序，零序等效图；5.理解正序定则，并且会用正序定则进行各种短路故障的计算；6.用simulink搭建出短路故障的模型，利用仿真结果与编程计算的结

果进行比较，并验证编程计算结果的正确性。潮流计算：1.在读懂程序的基础上画出潮流计算基本流程图；2.通过输入数据，对电力系统分析教材上的例11-5进行潮流计算并输出结果；设3.会用AutoCAD和multisim等软件画电力系统图和等效电路图。计短路计算：要1．在对称短路计算、简单不对称短路计算中都进行计算；求2．计算机语言自选；3．设计、编制、调试出相关的通用计算程序；4．输入／输出数据一律以文件格式形成；5．要求计算的题目：采用所编制的程序进行《电力系统分析上册》例6-3题，例8-1题。

2013年3月10日--3月12日查阅电力系统短路故障相关资料；2013年3月13日--3月15日网上查询电力系统的潮流计算和短路时故障的计算机算法资料；间2013年3月16--3月19日对相应的题目分析详细的解法，确定编安程语言，学习matlab编程和simulink建模；排2013年3月20--3月30日掌握matlab编写的算法，并上机调试；2013年4月6日--4月14日完成设计论文，并检查上交。[1]《电力系统分析》华中科技大学出版社何仰赞，温增银；[2]《电路原理》清华大学出版社汪建；参考[3]《Matlab/Simulink电力系统建模与仿真》机械工业出版社于群，曹娜；资[4]《电力系统分析学习指导书》中国电力出版社王葵；料[5]《matlab从入门到精通》人民邮电出版社胡晓东，董辰辉；[6]《电力系统故障的计算机辅助分析》重庆大学出版社米麟书等目录4

1电力系统图及初步分析1.1电力系统图及设计任务.1.2初步分析2牛顿-拉夫逊法简介.2.1概述2．2一般概念2.3潮流计算的修正方程2.4直角坐标表示的修正方程3程序设计3.1程序流程图5

输入原始数据形成节点导纳矩阵设电压初值e、f(0)(0)设迭代次数k=0计算误差向量∆P(k)、∆Q(k)、∆U2(k)是收敛否否求PU节点无功功率，求雅可比矩阵元素求平衡节点功率解修正方程，求解∆e(k)、∆f(k)求支路功率分布和损耗停机修正节点电压e(k+1)=e(k)-∆e(k)f(k+1)=f(k)-∆f(k)K=K+1K>Kmax不收敛停机3.2潮流计算程序运行结果如下：.二．三相短路计算2.1计算原理：利用节点阻抗矩阵计算短路电流注意：上述计算方法以及公式来源于电力系统分析上册P136-P1372.2三相短路计算流程图：6

2.3习题实例【例6-3】在如图2-3所示的电力系统中分别在节点1和节点5接入发电机支路，其标幺值参数为：.2.4三相短路计算程序及结果如下：三．不对称短路计算3.1不对称短路课程设计的题目3.2课程设计的设计任务及设计大纲3.3电力系统不对称故障时元件的序参数和等值电路错误!签。3.3.1电力系统不对称故障时用标幺值表示的各序等值电路未指定书签。3.4电力系统不对称故障时各序等值电路的化简与计算书签。3.4.1正序等值电3.4.2负序等值电3.4.3零序等值电3.5电力系统3.5.1理论分析路的化简计算路的化简计算路的化简计算不对称故障时元件参数的计算7

3.5.2各元件各序等值电路电抗标幺值的计算3.6电力系统不对称故障分析与计算3.6.1单相接地短路3.6.2两相直接接地短路3.6.3两相短路注释：以上程序中的计算公式都是根据正序等效定则得到的。3.7正序等3.8短路计算的matlab/simulink模型如下：3.9.1变压器和线路3.9.2短路模块和负载模块的参数设置3.9.3故障相短路设计总结.效定则的内容参数设置：相电流和相电压波形参考文献.附录一．潮流计算8

1电力系统图及初步分析1.1电力系统图及设计任务此电力系统图有AutoCAD2012软件画出网络各元件参数的标幺值如下：Z12=0.1+j0.4;y120=y210=j0.01538;z13=j0.13;k=1.1;z14=0.12+j0.5;y140=y410=j0.01920;z24=0.08+j0.4;y240=y420=j0.01413系统中节点1，2为PQ节点，节点3为P节点，节点4为平衡节点，已给定P1s+jQ1s=-0.3-j0.18,P2s+jQ2s=-0.55-j0.13P3s=0.5,V3s=1.10,V4s=1.050o容许误差为。试用牛顿法计算潮流分布1051.2初步分析潮流计算在数学上可归结为求解非线性方程组，其数学模型简写如下：2牛顿-拉夫逊法简介2.1概述9

牛顿－拉夫逊法是目前求解非线性方程最好的一种方法。这种方法的特点就是把对非线性方程的求解过程变成反复对相应的线性方程求解的过程，通常称为逐次线性化过程，就是牛顿－拉夫逊法的核心。牛顿-拉夫逊法的基本原理是在解的某一邻域内的某一初始点出发，沿着该点的一阶偏导数——雅可比矩阵J，朝减小方程的误差的方向前进一步，在新的点上再计算误差和雅可比矩阵，重复这一过程直到误差达到收敛标准，即得到了非线性方程组的解。因为越靠近解，偏导数的方向越准，收敛速度也越快，所以牛顿法具有二阶收敛特性。2．2一般概念对于非线性代数方程组(2－1)即fx,x,,x0i1,2,ni12n在待求量x的某一个初始计算值x0附件，将上式展开泰勒级数并略去二阶及以上的高阶项，得到如下的线性化的方程组x00(2－2)fxfx00上式称之为牛顿法的修正方程式。由此可以求得第一次迭代的修正量1fx(2－3)0x0fx010

将和相加，得到变量的第一次改进值。接着再从出xxx1x001发，重复上述计算过程。因此从一定的初值出发，应用牛顿法求x0解的迭代格式为(2－4)fxkfxxkk(2－5)xkxkxk1上两式中：是函数对于变量的一阶偏导数矩阵，即雅可比fxfxx矩阵J；k为迭代次数。由式（2－4）和式子（2－5）可见，牛顿法的核心便是反复形成求解修正方程式。牛顿法当初始估计值和方程的精确解足够接x0近时，收敛速度非常快，具有平方收敛特性。2.3潮流计算的修正方程运用牛顿－拉夫逊法计算潮流分布时，首先要找出描述电力系统的非线性方程。这里仍从节点电压方程入手，设电力系统导纳矩阵已知，则系统中某节点（i节点）电压方程为从而得进而有S•UYU•niiijjj1（2－6）•nPjQUYU0iiiijjj1式（2－6）中，左边第一项为给定的节点注入功率，第二项为由11

节点电压求得的节点注入功率。他们二者之差就是节点功率的不平衡量。现在有待解决的问题就是各节点功率的不平衡量都趋近于零时，各节点电压应具有的价值。由此可见，如将式（2－6）作为牛顿－拉夫逊中的非线性函数FX0，其中节点电压就相当于变量。建立了这种对应关系，就X可列出修正方程式，并迭代求解。但由于节点电压可有两种表示方式——以直角做表或者极坐标表示，因而列出的迭代方程相应地也有两种，下面分别讨论。2.4直角坐标表示的修正方程节点电压以直角坐标表示时，令•、，且将导•UejfjjjUejfiii纳矩阵中元素表示为YGjB，则式（2－7）改变为ijijij（2－7）PjQejfnGjBejf0iiiiijijjjj1再将实部和虚部分开，可得eGeBffGfBe0iijjijjiijjijjfGeBfeGfBe0iijjijjiijjijjPin（2－8）j1Qinj1这就是直角坐标下的功率方程。可见，一个节点列出了有功和无功两个方程。12

对于节点（），给定量为节点注入功率，记为、PQi1，2,,m1PiQ，则由式（2－8）可得功率的不平衡量，作为非线性方程iPPneGeBffGfBei（2－9）iiijjijjiijjijjj1QQnfGeBfeGfBeiiiijjijjiijjijjj1式中、PQ——分别表示第i节点的有功功率的不平衡量和无功功ii率的不平衡量。对于节点（），给定量为节点注入有功功率及im1,m2,,nPV电压数值，记为、，因此，可以利用有功功率的不平衡量和电PUii压的不平衡量表示出非线性方程，即有PPneGeBffGfBej（2－10）iiiijjijjiijjijj1U2U2efi22iii式中为电压的不平衡量。Ui对于平衡节点（），因为电压数值及相位角给定，所以im•也确定，不需要参加迭代求节点电压。UejfSSs个节点的系统只能列出个方程，其中有功功率方21n因此，对于nnm程个，无功功率方程个，电压方程个。将式（2－9）、n1m1式（2－10）非线性方程联立，称为n个节点系统的非线性方程组，且按泰勒级数在、（）展开，并略去高次项，得0e0i1,2,,n,imfii13

到以矩阵形式表示的修正方程如下。PHf1NLHJNLHJ1pHNL1pNHJ1nHNL1n1111112121p1p1nQeJ1111112121n1PHNNHNf221212222222p2p2p2n2n2n2QeJL21JL22JL2pJL2n2212（2－11）PfpHNSHRNSHRNHRNp1p1p1p2p2p2p2pppppppppnpnpnpU2eRSSpnpp1pPnfHNSHRNSHRNSHRNnnn1n1n2n2npnpnnnUR2Sennnn1n1n2n2npnpnnn上式中雅可比矩阵的各个元素则分别为将（2－11）写成缩写形式PHN212ff（－）QJLJeeRSU2对雅可比矩阵各元素可做如下讨论：当时，对于特定的，只有该特定点的和是变量，于是雅可jijfeii比矩阵中各非对角元素表示为当时，雅可比矩阵中各对角元素的表示式为ji由上述表达式可知，直角坐标的雅可比矩阵有以下特点：1)雅可比矩阵是阶方阵，由于2n1、等等，所以jiHHNNijijji它是一个不对称的方阵。14

2)雅可比矩阵中诸元素是节点电压的函数，在迭代过程中随电压的变化而不断地改变。3)雅可比矩阵的非对角元素与节点导纳矩阵中对应的非对角元YB素有关，当中的为零时，雅可比矩阵中相应的H、N、J、YYBijijijijL也都为零，因此，雅可比矩阵也是一个稀疏矩阵。ij3程序设计3.1程序流程图3.2潮流计算程序运行结果如下：请输入节点数：n=4请输入支路数:n1=4请输入平衡母线节点号isb=4请输入误差精度pr=0.00001请输入由之路参数形成的矩阵B1=[120.1+0.4i0.3056i10;130+0.3i01.10;140.12+0.5i0.0382i10;240.08+0.4i0.02826i10]请输入各节点参数形成的矩阵B2=[0-0.3-0.18i1002;0-0.55-0.13i1002;00.5+0i11.103;0015

11.0501]节点号和对地参数:X=[10;20;30;40]导纳矩阵Y=1.0421-7.4054i-0.5882+2.3529i0+3.3333i-0.4539+1.8911i-0.5882+2.3529i1.0690-4.5899i0-0.4808+2.4038i0+3.3333i00-3.3333i0-0.4539+1.8911i-0.4808+2.4038i00.9346-4.2617i初始功率参数OrgS=0.0000-0.1719-0.0000-0.166900功率和电压的不平衡量DetaS=16

-0.3000-0.0081-0.55000.03690.50000第一次迭代的雅克比矩阵Jacbi=7.57731.0421-2.3529-0.5882-3.33330-1.04217.23350.5882-2.35290-3.3333-2.3529-0.58824.75681.0690000.5882-2.3529-1.06904.422900-3.33330003.33330000002.0000第一次迭代的修正方程DetaU=-0.023617

-0.0005-0.1232-0.01860.12640节点1的电压是0.9995-0.0236i节点2的电压是0.9814-0.1232i节点3的电压是1.0000+0.1264i节点4的电压是1I=-0.2959+0.1801i18

-0.5294+0.1331i0.5000-0.0018i雅克比矩阵Jacbi=7.55690.9205-2.3378-0.6435-3.3315-0.0787-1.51247.19670.6435-2.33780.0787-3.3315-2.2368-0.86724.50601.0852000.8672-2.2368-2.14414.239800-3.33330.4213003.33150.078700000.25282.0000修正方程DetaU=0.0006-0.01170.0001-0.02150.002319

-0.0083I=-0.2993+0.1891i-0.5464+0.2039i0.5058-0.0132i雅克比矩阵Jacbi=7.47990.9009-2.3106-0.6353-3.2925-0.0769-1.49947.10160.6353-2.31060.0769-3.2925-2.1863-0.85434.47831.0448000.8543-2.1863-2.13774.070500-3.30570.4289003.29250.076900000.25741.9834修正方程DetaU=1.0e-003*-0.011920

-0.1922-0.0131-0.48090.0367-0.0420I=-0.2994+0.1893i-0.5468+0.2057i0.5060-0.0137i雅克比矩阵Jacbi=7.47860.9007-2.3101-0.6352-3.2919-0.0769-1.49947.10010.6352-2.31010.0769-3.2919-2.1851-0.85414.47791.0440000.8541-2.1851-2.13764.066500-3.30560.4291003.29190.076921

00000.25741.9834修正方程DetaU=1.0e-006*-0.0103-0.0909-0.0039-0.25680.0092-0.0028迭代次数为4节点1的电压是0.9876-0.0231i节点2的电压是0.9595-0.1231i22

节点3的电压是0.9917+0.1287i节点4的电压是1可见：上述计算结果，与电力系统分析教材上的结果基本一致。我们也可以用matlab/simulink中提供的图形用户分析界面powergui模块以及SimPowerSystem模块搭建模型，进行潮流计算分析，同样可以验证上述结果。另外，也可以运用中国电力科学院开发的电力系统分析综合程序软件PSASP进行潮流计算。由于时间有限，在此不再赘述。二．三相短路计算2.1计算原理：利用节点阻抗矩阵计算短路电流如图3-1所示假定系统中的节点f经过渡阻抗zf发生短路。这个过渡阻抗zf不参与形成网络的节点导纳矩阵，如果保持故障处的边界条件不变，把网络的原有部分同故障支路分开图3-1因此，对于正常的网络状态而言，发生短路相当于在故障节点f增加了一个注入电流-If，因此，网络中任一节点i的电压可以表示为（3-1）式中，G为网络内有源i节点的集合。VZIZIijjiffjG由上式可见，任一节点i的电压都由两项叠加而成，第一项表示当注入电流If=0时由网络内所有电源在节点i产生的电压，也就是短路前瞬间正常运行状态下的节点电压，这是节点电压的正常分量，记作Vi（0）。第二项是当网络中所有电流源都断开，电压源都短接时，仅仅由短路电流If在节点i产生的电压，这就是节点电压的故障分量。23

由此可知，式（3-1）又可表示为（3-2）VVZI(0)式（3-2）也适用于故障点f，于是有iiiff（3-3）VVZIff(0)式中，是故障点f的自阻抗，也称为输入阻抗f。ffZff根据边界条件（3-4）VzI0fff由式（3-3）和（3-4）可以得出（3-5）V(0)IfZzf即可求出短路电流。fff注意：上述计算方法以及公式来源于电力系统分析上册P136-P1372.2三相短路计算流程图：用公式（6-10）计算短路电流IfIf=1/（Zff+zf）2.3习题实例【例6-3】在如图2-3所示的电力系统中分别在节点1和节点5接入用公式（6-11）计算各点电压Vi=1-Zif/(Zff+zf)发电机支路，其标幺值参数为：。在节点3发生三相短路，计算短路电EE1.0,zj0.15,zj0.22用公式（6-9）计算指定支路的电路1515流及网络中的电流分输布。线路的电阻入数据计算节点和电容阻抗矩略去不计，阵If列元变压器的标幺变比等于1。各元件参数的标幺值如下：输出结果图2-3电力系统等值网络图图2-4三相短路时的等值网络图(用multisim软件可画出)24

Y=-j13.87160j9.5238000-j8.33330j4.76190j9.52380-j15.2329j2.2960j3.44400j4.7619j2.2960-j10.9646j3.93602.4三相短路计算程序及结果如下：n=input('请输入短路节点号f=');Y=[0-16.905j,9.5238j,0,0,0;0+9.5238j,37.4084j,15.3846j,12.5000j,0;0,15.3846j,-35.3846j,20.000j,0;0,12.5000j,20.000j,-37.9348j,5.4348j;0,0,0,5.4348j,-9.9802j];disp('导纳矩阵Y='),disp(Y)Z=inv(Y);%求逆矩阵，得到阻抗矩阵disp('阻抗矩阵Z='),disp(Z)disp('短路电流If为')If=1/0.1860i25

disp('故障后，各节点电压为')V1=1-0.0902i*IfV2=1-0.1533i*IfV3=0V4=1-0.1611i*IfV5=1-0.0877i*Ifdisp('故障后，各支路电流为')I54=(V5-V4)/0.184iI43=(V4-V3)/0.05iI23=(V2-V3)/0.065iI12=(V1-V2)/0.105iI24=(V2-V4)/0.08i运行结果如下：请输入短路节点号f=3导纳矩阵Y=26

0-16.9050i0+9.5238i0000+9.5238i0+37.4084i0+15.3846i0+12.5000i000+15.3846i0-35.3846i0+20.0000i000+12.5000i0+20.0000i0-37.9348i0+5.4348i0000+5.4348i0-9.9802i阻抗矩阵Z=0+0.0545i0-0.0082i0-0.0077i0-0.0074i0-0.0040i0-0.0082i0-0.0146i0-0.0137i0-0.0131i0-0.0071i0-0.0077i0-0.0137i0+0.0288i0+0.0116i0+0.0063i0-0.0074i0-0.0131i0+0.0116i0+0.0305i0+0.0166i0-0.0040i0-0.0071i0+0.0063i0+0.0166i0+0.1093i短路电流If为If=0-5.3763i故障后，各节点电压为27

V1=0.5151V2=0.1758V3=0V4=0.1339V5=0.5285故障后，各支路电流为I54=0-2.1447iI43=0-2.6774i28

I23=0-2.7047iI12=0-3.2309iI24=0-0.5242i可见：此计算结果与电力系统分析教材上的结果一样。三．不对称短路计算3.1不对称短路课程设计的题目电力系统简单结构图如图3.1所示。图3.1电力系统在*)K点发生不对称短路，系统=120MVA，U=10.5kV，次暂态电简单结构图(为简洁，不加下标动势标幺值1.67，各元件参数如下：发电机G1：Snn次暂态电抗标幺值X为0.9，负序电抗标幺值为0.45；=60MdVA，U%=10.5变压器T1：SnK=60MVA，U%=10.5变压器T2：Sn线路L=105km，单位长度电抗=60MVA，X=1.2，X=0.35=40MVA，X=1.2，X=0.35Kxx=0.4Ω/km，=31,x10负荷L1：Sn负荷L2：S12n12=120MVA和U取SBB为所在级平均额定电压Vav。29

3.2课程设计的设计任务及设计大纲⑴选择110kV为电压基本级，画出用标幺值表示的各序等值电路。并求出各序元件的参数（要求列出基本公式，并加说明）。⑵化简各序等值电路并求出各序总等值电抗。⑶K处发生单相直接接地短路，列出边界条件并画出复合相序图。求出短路电流。⑷设在K处发生两相直接接地短路，列出边界条件并画出复合相序图。求出短路电流。⑸讨论正序定则及其应用。并用正序定则直接求在K处发生两相直接短路时的短路电流。⑹思考提高：用Matlab仿真并比较结果。⑺附录：要画出完整各序等值电路图以及给出参数计算的程序。3.3电力系统不对称故障时元件的序参数和等值电路要求：选择110kV为电压基本级，画出用标幺值表示的各序等值电路。并求出各序元件的参数（要求列出基本公式，并加说明）。3.3.1电力系统不对称故障时用标幺值表示的各序等值电路图3.2电力系统不对称故障时用标幺值表示的正序等值电路图3.3电力系统不对称故障时用标幺值表示的负序等值电路图3.4电力系统不对称故障时用标幺值表示的零序等值电路3.4电力系统不对称故障时各序等值电路的化简与计算要求：化简各序等值电路并求出各序总等值电抗（戴维南等效电路）。3.4.1正序等值电路的化简计算图3.5正序等值电路首先求整个网络对短路点的正序等值电动势和正序等值电抗。在图3.5中，将支路1和支路5并联得支路7，它的电抗和电动势分别为：将支路7、2、4串联，得支路9，它的电抗为：将支路3、6串联得支路8，其电抗为：30

将支路8、9并联得：图3.6正序等值网络化简后的电路图3.4.2负序等值电路的化简计算图3.7负序等值电路首先求整个网络对短路点的负序等值电抗。在图3.7中，将支路1和支路5并联得支路7，它的电抗分别为：将支路7、2、4串联，得支路9，它的电抗为：将支路3、6串联得支路8，其电抗为：将支路8、9并联得：图3.8负序等值网络化简后的电路图3.4.3零序等值电路的化简计算图3.9零序等值电路将支路1和支路4串联得：图3.10负序等值网络化简后的电路图3.5电力系统不对称故障时元件参数的计算3.5.1理论分析进行电力系统计算时，采用有单位的阻抗、导纳、电压、电流、功率等的相对值进行运算、称为有名制。在作整个电力系统的等值网络图时，必须将其不同电压级的各元件参数阻抗、导纳以及相应的电压、电流归算至同一电压等级—基本级。而基本级一般电力系统中取最高电压级。式中，K、K、…K为变压器的变比；R’、X’、G’、B’、分别为归算12n前的有名值；R、X、G、B、分别为归算后的有名值。进行电力系统计算时，采用没有单位的阻抗、导纳、电压、电流、功31

率等的相对值进行运算、称为标幺制。标幺值的定义为：本设计中MVA，和所在级平均额定电压相等。在电力S120UU系统计算中，用平均额定电压之比代替变压器的实际变比时，元件参avbBB数和变量的标幺值的计算可大为简化。所以将元件参数和变量归算至基本级为：而求取电力系统各元件（发电机G、变压器T、电力线路l、电抗器L）电抗的标么值的计算公式如下：3.5.2各元件各序等值电路电抗标幺值的计算选取110kV为电压基本级，在电力系统暂态分析中，等值电路中的电阻可以忽略不计，所以有以下结论。⑴发电机G1的各序等值电路电抗标幺值：S发电机的正序电抗标幺值B。XXd1SG(1)*NS发电机的负序电抗标幺值B。XXd2SG(2)*N由于变压器的连接方式为连接，所以零序网络与发电机是断开Δ/Y的，无零序电流流过，其零序电抗为0。MATLAB程序如下：%求发电机参数的标幺值,计算公式：X=Xd1*(SB/SGN)clearSn=120;SB=120;Xdc1=0.9;Xdc2=0.45;XG1b=Xdc1*(SB/Sn);disp('一.发电机1的电抗值XG1b='),disp(XG1b)32

XG2b=Xdc2*(SB/Sn);disp('发电机2的电抗值XG2b='),disp(XG2b)程序运行结果为：一.发电机1的电抗值XG1b=0.9000发电机2的电抗值XG2b=0.4500即有发电机的正序电抗标幺值0.9，负序电抗标幺值X。X0.45G(1)*G(2)*⑵变压器T1和T2的各序等值电路电抗标幺值：变压器T1的正序电抗标幺值U%SB。100SXkT(1)*N变压器T1的负序电抗标幺值U%S100SB。XkT(2)*N变压器T1的零序电抗标幺值U%SB。100SXkT(0)*N由于变压器T1和变压器T2的参数一所以变压器T2的正序电样，、零序电抗的标幺值与变压器T1的正序电抗抗、负序电抗、负序电抗、零序电抗相等。MATLAB程序如下：33

%求变压器T的各序等值电路电抗的参数，计算公式：XT=Uk%/100*(SB/STN)clearST1=60;ST2=60;SB=120;Uk1=10.5;Uk2=10.5;XT1=(Uk1/100)*(SB/ST1);disp('二.变压器T的各序电抗XT1='),disp(XT1)XT2=(Uk2/100)*(SB/ST2);disp('XT2='),disp(XT2)XT0=XT1;%由于变压器是静止电器，所以各序参数相等disp('XT0='),disp(XT0)程序运行结果为：二.变压器T的各序电抗XT1=0.2100XT2=0.2100XT0=34

0.2100即有变压器T1(T2)的正序电抗标幺值，变压器T1(T2)的X0.21，变压器T1(T2)零序电抗标幺值T(1)*负序电抗标幺值XX0.21。T(2)*0.21T(0)*l⑶电力线路的各序等值电路电抗标幺值：lS电力线路的正序电抗标幺值l(1)*。。Xxl1BU2avnl电力线路的负序电抗标幺值SXxlBU2l(2)*1avnlSS电力线路的零序电抗标幺值L(0)*。Xxl3xlBBUU2012avnavnMATLAB程序如下：%输电线1的各序等值电路的电抗标幺值计算，计算公式：Xl1=Xl2=x0*(SB/Uav^2),Xl0=3*Xl1clearSB=120;x0=0.4;L=105;Uav=115;Xl1b=x0*L*(SB/Uav^2);disp('三.输电Xl1b='),disp(Xl1b)线的各序电抗Xl2b=x0*L*(SB/Uav^2);disp('Xl2b='),disp(Xl2b)35

Xl0b=3*x0*L*(SB/Uav^2);disp('Xl0b='),disp(Xl0b)程序运行结果为：三.输电线的各序电抗Xl1b=0.3811Xl2b=0.3811Xl0b=1.1433即有电力线路Xl的正序电抗标幺值，电力线路。1.1433l的负序电X0.3811l(1)*抗标幺值，电力线路l的零序电抗标幺值0.3811Xl(2)*L(0)*⑷负荷L1的各序等值电路电抗标幺值：负荷L1的正序电抗标幺值SB。SXX111L(1)*N负荷L1的负序电抗标幺值SB。SXX121L(2)*N负荷L2的各序等值电路电抗标幺值：负荷L1的正序电抗标幺值SB。SXX212L(1)*N36

负荷L1的负序电抗标幺值XX2L(2)*22SB。SNMATLAB程序如下：%负荷1的各序。计算公式：X1L=X1*(SB/Sn)等值电路的电抗标幺值clearSB=120;Sn=60;X11=1.2;X12=0.35;X1L1b=X11*(SB/Sn);disp('负荷1各序参数X1L1b='),disp(X1L1b)X1L2b=X12*(SB/Sn);disp('XlL2b='),disp(X1L2b)%负荷2的各序。计算公式：X2L=X1*(SB/Sn)等值电路的电抗标幺值SB=120;Sn=40;X21=1.2;X22=0.35;X2L1b=X21*(SB/Sn);disp('负荷2各序参数X2L1b='),disp(X2L1b)X2L2b=X22*(SB/Sn);disp('X2L2b='),disp(X2L2b)37

程序运行结果为：负荷1各序参数X1L1b=2.4000XlL2b=0.7000负荷2各序参数X2L1b=3.6000X2L2b=1.0500即负荷L1的正序电抗标幺值，负荷L1的负序电抗标幺值L2的正序电X抗标幺2.值4，负荷L2的负序3.62L(1)*。负荷XX0.71L(1)*X电抗标幺值。1.051L(2)*2L(2)*由于变压器的连接方式为连接，所以零序网络与负荷是断开的，Δ/Y无零序电流流过，其零序电抗为0。3.6电力系统不对称故障分析与计算要求：若K处发生单相直接接地短路，列出边界条件并画出复合相序图，求出短路电流；若在K处发生两相直接接地短路，列出边界条件并画出复合相序图，求出短路电流。电力系统中发生不对称短路时，无论是单相接地短路、两相短路还是两相接地短路，只是在短路点出现系统结构的不对称，而其它部分三相仍旧是对称的。根据对称分量法列a相各序电压方程式为上述方程式包含了六个未知量，必须根据不对称短路的具体边界条件38

列出另外三个方程才能求解。3.6.1单相接地短路图3.6.1单相接地短路⑴边界条件当电力系统中的K点发生单相（A相）直接短路接地故障时，其短路点的边界条件为A相在短路点K的对地电压为零，B相和C相从短路点流出的电流为零，即：⑵复合相序图将边界条件用对称分量法表示为：由上式可以作出单相接地短路的复合序网络图如图3.6.2所示。图3.6.2单相接地短路的复合序网络(Zf=0)所以有：MATLAB程序如下：%单相接地短路时的短路电流计算Uf0=0.95;Zkk1=0.83;Zkk2=0.44;Zkk0=0.78;SB=120;UB=115;If1=SB/(UB*sqrt(3))*Uf0/(Zkk1+Zkk2+Zkk0)If=3*If139

程序运行结果为：If1=0.2792If=0.8376即发生单相直接接地短路时，其短路电流If=0.8376。3.6.2两相直接接地短路⑴边界条件当电力系统中的K点发生单相（B相和C相）直接短路接地故障时，其短路点的边界条件为：图3.6.3两相直接接地短路⑵复合相序图将边界条件用对称分量法表示为：由上式可以作出两相接地短路的复合序网络图如图3.6.4所示。图3.6.4单相接地短路的复合序网络(Zf=Zg=0)(设各序阻抗为纯阻抗)：由此图直接可以求其序电流为进而推出：40

MATLAB程序如下：%两相接地短路时的短路电流计算Uf0=0.95;Zkk1=0.83;Zkk2=0.44;Zkk0=0.78;SB=120;UB=115;IB=SB/(UB*sqrt(3));Zzeta=Zkk2*Zkk0/(Zkk2+Zkk0);m(1,1)=sqrt(3)*sqrt(1-(Zkk2*Zkk0/((Zkk2+Zkk0)*(Zkk2+Zkk0))));If1=SB/(UB*sqrt(3))*Uf0/(Zkk1+Zzeta);disp('两相接地短路电流的正序分量If1='),disp(If1)If=m(1,1)*If1;disp('两相接地短路电流If='),disp(If)程序运行结果为：两相接地短路电流的正序分量If1=0.5150两相接地短路电流If=41

0.7824即发生两相直接接地短路时，其短路电流正序分量If1=0.5150kA，短路电流If=0.7824kA。3.6.3两相短路Matlab程序如下：%两相短路时的短路电流计算Uf0=0.95;Zkk1=0.83;Zkk2=0.44;Zkk0=0.78;SB=120;UB=115;If1=SB/(UB*sqrt(3))*Uf0/(Zkk1+Zkk2);disp('两相短路电流If=sqrt(3)*If1;disp('两相短路电流的正序分量If1='),disp(If1)If='),disp(If)程序运行结果为：两相短路电流的正序分量If1=0.4507两相短路电流If=42

0.7806注释：以上程序中的计算公式都是根据正序等效定则得到的。3.7正序等效定则的内容三种简单不对称短路时短路电流正序分量的通式为：式中，称附加阻抗。Z(n)su正序等效定则：在简单不对称短路的情况下，短路点电流的正序分量与在短路点后每一相中加入附加阻抗而发生三相短路的电流相Z(n)等。su表3.7.1各种类型短路时附加阻抗值由于故障相短路点短路电流的绝对值与它的正序分量的绝对值成正比，即：式中，是比例系数。其值视短路的种类而异。各种简单短路的m4.7.2m(n)(n)值见表。表4.7.2各种类型短路时比例系数值根据以上的讨论，可以得到一个结论：简单不对称短路电流的计算，归根结底，不外乎先求出系统对短路点正序、负序和零序等值阻抗（或电抗）；再根据短路种类的不同而组成附加阻抗，将它接入短路点Z(n)的正序等值电抗；然后就像计算三相短路一样，计算出短路点的正序su电流，从而可以算出其他各序电流、各序电压，及短路点的三相电流和三相电压。这样三相短路电流的各种计算方法，也适用于不对称短路时正序电流的计算。3.8短路计算的matlab/simulink模型如下：3.9所用到的模块及参数设置：(1)三相电源120MVA，10.5KV；(2)PI型输电线路；(3)双绕组变压器；模块(4)三相故障Three-PhaseFault;(5)三相电压电流测量模块(充当母线A，B)；(6)电力系统图形用户界面模块powergui;(7)短路电流multimeter,scope等。获取模块仿真参仿真算数法设置为ode15s：在0.2s。时刻发生短路故障，仿真结束时间为1.0s，3.9.1变压器和线路参数设置：变压器的参数计算方法：43

V%V*100S10.510.52*100.19293753一次侧电抗：2*103XT1*10060000s1NNVs%V*100S10.51152*10323.14375*10060000二次侧电抗：2*103XT22NN即电感：LX/(2*pi*f)23.14375/(2*3.14*50)0.07371T2图3.9.1变压器模块参数设置图3.9.2输电线路模块参数设置3.9.2短路模块和负载模块的参数设置输电线路L：正序电抗；零序电抗X0.4/kmL1；线路长度L1=105km；频率为50Hz，X3*X3*0.41.2/kmL0L1其他采用默认值；线路2的参数设置方法相同。负载1参数：所接电压等级V=10.5kv,容量60MW；负载2的参数设置方法相同图3.9.3短路模块参数设置图3.9.4负载模块参数设置3.9.3故障相短路相电流和相电压波形分析：由波形可知：短路故障相电流If的幅值约为1200A，转换为果If=0.8376非常接近，从而验证了计算的正确性。其他的两相短路，两相接地短路，三相短路的波形在此不再赘述，方法一样，只是需要修改一下，三相短路模块的参数设置。有效值If==848A=0.848kA。此结果与上面编程计算的结1200/2设计总结本课程设计解决的核心问题有两个：一是对给定系统进行潮流计算，二是对简单系统进行短路故障计算，其中包括对称短路和不对称短路。在计算三相对称短路时，根据节点阻抗矩阵计算短路电流；在计算不对称短路故障时，根据对称分量法进行解析，即把该网络分解为44

正，负，零序三个对称序网，这三组对称序分量可分别按对称的三相电路分解，然后将其结果叠加起来。最后根据正序等效定则得到各种类型不对称短路故障的短路电流。进行潮流计算时，要会用AutoCAD或PSASP画出电力系统图，在看懂读懂程序的基础上，修改参数，明白个输入矩阵元素的物理意义，弄清楚个计算公式的来源，然后一步步调试程序，最终得到自己想要的结果。求解对称短路时，先用multisim10.0软件画出等效电路图，然后用matlab编程计算。求解不对称短路，后制定各序网络。根据不同的故障类型，确定出以相分量界条件，进而列出以序分量表示的边界条件，按边界条件将三个序网联合成复合网，由复合网求出故障处各序电流和电压，进而合成三相电流电压。最后并且用matlab/simulink建模仿真方法得到单相短路电流及电压，验证了程序计算结果的正确性。首先应该计算各元件的序参数和画出等值电路。然法表示的边通过上面两种方法放入计算进行对比可以得出，计算机程序法比较解析法具有计算过程简单及结果更精确，通用性更强的优点。而建模仿真的方法虽然只管，实现容易，方法简便，当具体参数设置困难，导致计算精确度不高。总之，这次设计过程中也曾在Matlab编程中，matpower编程和画电力系统图的过程中遇到了一些问题，在老师的指导和同学的帮助下以及网上查询资料的完成了，在查阅资料的过程中也丰富了识，但不可避免还是存在一些未解决的问题。最明显的问题是AutoCAD2012设计电气主接线的方法不够熟练，用PSASP软件进自己的知行电力系统潮流计算和短路计算时，总出现运行错误，有待解决。在此希望通过自己的进一步学习，丰富自己的专业技能和动手能力，去解决更多的实际问题，同时也很感谢这次课程设计的机会，让我受益匪浅。参考文献[1]何仰赞，温增银.《电力系统分析》.武汉：华中科技大学出版社(第三版)2011年；[2]汪建.《电路原理》.武汉：清华大学出版社,2007年；[3]于群，曹娜.《Matlab/Simulink电力系统建模与仿真》北京：机械工业出版社,2012年1月；[4]王葵.《电力系统分析学习指导书》北京：中国电力出版社,2009年1月；45

[5]胡晓东，董辰辉.《matlab从入门到精通》北京：人民邮电出版社,2010年6月；[6]米麟书，刘芳宁《电力系统故障的计算机辅助分析》重庆：重庆大学出版社,1990年8月；[7]吴际舜，侯志俭《电力系统潮流计算的计算及方法》上海：上海交通大学出版社,1999年4月；[8]马银莎《基于matlab的电力系统潮流计算》东北电力大学课程设计论文,2011年。附录1.潮流计算的matlab编辑的m语言源程序：%潮流计算源程序如下(交互式)clear;clc;n=input('请输入节点数：n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡母线节点号isb=');pr=input('请输入误差精度pr=');46

B1=input('请输入由之路参数形成的矩阵B1=');B2=input('请输入各节点参数形成的矩阵B2=');X=input('节点号和对地参数:X=');%X=input('节点号和对地参数:X=');Y=zeros(n);Times=1;%置迭代次数为初始值%创建节点导纳矩阵fori=1:n1ifB1(i,6)==0%不含变压器的支路，B1中每行的最后一个元素p=B1(i,1);%节点首端号q=B1(i,2);%节点末端号Y(p,q)=Y(p,q)-1/B1(i,3);Y(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);47

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)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));endenddisp('导纳矩阵Y=');disp(Y)OrgS=zeros(2*n-2,1);%有两个PQ节点，即是DetaS=zeros(2*n-2,1);%将OrgS、DetaS初始化为全零矩阵%创建OrgS，用于存储初始功率参数4个未知量h=0;48

j=0;fori=1:n%对PQ节点的处理ifi~=isb&B2(i,6)==2%B2中每行的最后一个元素表示节点类型，2表示PQ节点h=h+1;forj=1:n%公式来源于电力系统分析下册P57，式11-45OrgS(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))*imag(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(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfori=1:n%对PV节点的处理，注意这时不可再将h初始化为0ifi~=isb&B2(i,6)==3%3表示PV节点49

h=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))*imag(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(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendenddisp('初始功率参数OrgS=')disp(OrgS)%创建PVU用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;fori=1:n50

ifB2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endend%disp('PV节点的初始电压PVU=')%disp(PVU)%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量h=0;fori=1:n%对PQ节点的处理ifi~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);%来源于式11-46DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);end51

endt=0;fori=1:n%对PV节点的处理，注意这时不可再将h初始化为0ifi~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);%来源于式11-47DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;endenddisp('功率和电压的不平衡量DetaS=')disp(DetaS)%创建I，用于存储节点电流参数i=zeros(n-1,1);52

h=0;fori=1:nifi~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endend%disp('节点电流参数I=')%disp(I)%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;fori=1:n%对PQ节点的处理ifB2(i,6)==253

h=h+1;forj=1:nifj~=isbk=k+1;ifi==j%对角元素的处理Jacbi(2*h-1,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)=-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%非对角元素的处理i！=j时，来源于式11-49Jacbi(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);54

Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endifk==(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendk=0;fori=1:n%对PV节点的处理ifB2(i,6)==3h=h+1;forj=1:nifj~=isb55

k=k+1;ifi==j%对角元素的处理Jacbi(2*h-1,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(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;endifk==(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换56

行k=0;endendendendenddisp('第一次迭代的雅克比矩阵Jacbi=')disp(Jacbi)%实际上求出的是-Jacbi%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;disp('第一次迭代的修正方程DetaU=')disp(DetaU)%修正节点电压57

j=0;fori=1:n%对PQ节点处理ifB2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfori=1:n%对PV节点的处理ifB2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB2;forno=1:n58

Vn(no)=B2(no,3);str1=sprintf('节点%d的电压是%d',no);disp(str1);disp(Vn(no));endwhileabs(max(DetaS))>pr%控制计算精度OrgS=zeros(2*n-2,1);%初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0;j=0;fori=1:nifi~=isb&B2(i,6)==2%对PQ节点进行处理h=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))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));59

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(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfori=1:nifi~=isb&B2(i,6)==3%对PV节点进行处理h=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))*imag(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(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end60

endendOrgS%创建功率不平衡量DetaSh=0;fori=1:nifi~=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;fori=1:nifi~=isb&B2(i,6)==361

h=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))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;endendDetaS%创建Ii=zeros(n-1,1);h=0;fori=1:nifi~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end62

endI%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;fori=1:nifB2(i,6)==2h=h+1;forj=1:nifj~=isbk=k+1;ifi==jJacbi(2*h-1,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(63

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);endifk==(n-1)k=0;endendend64

endendk=0;fori=1:nifB2(i,6)==3h=h+1;forj=1:nifj~=isbk=k+1;ifi==jJacbi(2*h-1,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