电力系统分析课程设计1111

是否

信息工程学院课程设计任务书学号030740625030740626学生姓名专业（班级）0307406班设计题目电力系统潮流计算设计技术参数变压器的阻抗R； X，励磁损耗△S1。线路的阻抗R、X:'无功损耗△Q°o变压器及线路的功率损耗△S 、AS'1变压器及线路的电压降落Ak；、AV；各节点的功率和电压 T1 '1设计要求作给出电路的潮流计算，即算出各节点的电压和功率。两人合作，一人手算（夏丹丹），一人机算（陈志彪）。工作量需求分析（2页）总体设计（1页）详细设计（14页）总结（1页）工作计划整体分析电路，向老师和同学请教，对电路图有宏观把握，并且有初步的设计思路和方法。作出等值电路，计算出各参数，并且有清晰的思路和准确的方法。进行潮流计算。与MATLAB仿真比较，找出问题所在并分析。总结与反思参考资料何仰赞温增银《电力系统分析》第三版/view/627420.王守相刘玉田 电力系统潮流计算研究现状--《山东电力技术》1996年05期飞思科技产品研发中心.MATLAB基础与提高.北京：电子工业出版社，2005.5.StephenJ.Chapman.MATLABProgrammingforEngineersSecondEdition.国外高校电子信息类优秀教材6.刘同娟.MATLAB在电路分析中的应用.电气电子教学学报.2002指导教师签字教研室主任签字潮流计算是电力网络设计及运行中最基本的计算，对电力网络的各种设计方案及各种运行方式进行潮流计算，可以得到各种电网各节点的电压，并求得网络的潮流及网络中各元件的电力损耗，进而求得电能损耗。本文以一个具体例子分析潮流计算的具体方法，并给出MATLAB仿真。关键词：电力系统潮流计算Mmatlab仿真TOC\o"1-5"\h\z\o"CurrentDocument"1需求分析 61.1所给电路图及设计要求 6\o"CurrentDocument"1.2潮流计算的意义 6\o"CurrentDocument"2总体设计 8\o"CurrentDocument"2.1等值电路图 8\o"CurrentDocument"2.2设计思路及方法 8\o"CurrentDocument"3详细设计 9\o"CurrentDocument"3.1参数计算 9\o"CurrentDocument"3.2潮流计算 123.2MATLAB仿真 14\o"CurrentDocument"4总结 23\o"CurrentDocument"参考文献 2411需求分析1.1二级标题10MW

A10kV而=157kW度=73WI%=0.5V%=105Is32km25MW25MWYNd11◎I1(MWV^J<12MWYNd114x6MWVA*I%=0.3580冲％=1070km10MW

A10kV而=157kW度=73WI%=0.5V%=105Is32km25MW25MWYNd11◎I1(MWV^J<12MWYNd114x6MWVA*I%=0.3580冲％=1070km110VYN,d112x20MWV2gMW而=186kW度=8束WI%=0.5V%=105Vs10kV而=157kW华°=73W、*-…0110V10kV15MW4x15MWx=0.136xd=0.162x^=0.073cosp=0.8IN4x6MWx=0.134XXd=0.1612 x=0.06cosp=0.85VN慎=50WI%=0.55V%=105vs 一而=132kW寸=63WI%=0.55V%=10580MWV(1「2)%=6.55MWg3)” 、ruJ，V(2)%=1751s(l-3)1(kVYNd112x10MWVA30kmI%=0.5V%=105、s3x12MWx=0.136xd=0.1612x=0.075cosp=0.8110VYNd112x16MWVA6华华°=12kWI%=0.35V%=105Q1x5MWx=0.128x!=0.1542 .x=0.054cosp=0.85IN35MYN^Y,d11•2x10MV410kV2x2MWx-=0.128£=0.1572x=0.0591cosp=0.8N1.2潮流计算的意义潮流计算是研究电力系统稳态运行情况的一种基本电气计算，常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。通过潮流计算可以判断电网母线电压、支路电流和功率是否越限，如果有越限，就应采取措施，调整运行方式。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置定整计算、电力系统故障计算和稳定计算等提供原始数据。具体表现在以下方面：（1）在电网规划阶段，通过潮流计算，合理规划电源容量及接入点，合理规划网架，选择无功补偿方案，满足规划水平的大、小方式下潮流交换控制、调峰、调相、调压的要求。在编制年运行方式时，在预计负荷增长及新设备投运基础上，选择典型方式进行潮流计算，发现电网中薄弱环节，供调度员日常调度控制参考，并对规划、基建部门提出改进网架结构，加快基建进度的建议。正常检修及特殊运行方式下的潮流计算，用于日运行方式的编制，指导发电厂开机方式，有功、无功调整方案及负荷调整方案，满足线路、变压器热稳定要求及电压质量要求。预想事故、设备退出运行对静态安全的影响分析及作出预想的运行方式调整方案。总结为在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。2.1等值电路图• 1 I4bn rZT2T T25MW2总体设计ZL2A2.1等值电路图• 1 I4bn rZT2T T25MW2总体设计ZL2ASo2JAQL2jAQL3,25MWASo430MWe120MWjAQo3jAQo3JAQL21d土ZL4IIIjAQL4*ZT580MWm5MWjAQL5QzT5-335MWm5MWjAQL5QzT5-335MWEjAQL5 gASo6, •Dzt615MW2.2设计思路及方法潮流计算（loadflowcalculation）根据电力系统接线方式、参数和运行条件计算电力系统稳态运行状态下的电气量。通常给定的运行条件有电源和负荷节点的功率、枢纽点电压、平衡节点的电压和相位角。待求的运行状态量包括各节点电压及其相位角和各支路（元件）通过的电流（功率）、网络的功率损耗等。此题先按电网平均额定电压VN计算功率损耗，再利用上面求得的功率分布，从电源点开始，顺着功率的传送方向，依次计算各段线路的电压降落，求出各节点电压。对于首端电压和功率已知的线路，只需顺着电路的方向算各个变压器及线路的电压降落和功率损耗，以此求得各节点的电压和功率。3详细设计3.1参数计算变压器1：-APV273x1152(16x103)2x103=3.771。R SNx103T1SNX=^S^Vx10_10.5x1152-x10=86.79。T1 SN励磁功率：I%AS=AP+j—16x103S=(0.0157+j0.08)MVA01 0 100N线路一设r0=0.21。/km，*=0.42Q/km，b0=2.74x10%S/kmR=rxl=0.21x32=6.72。X=xxl=0.42x32=13.44。B=2.74x10-6x32S/kmL1产生的充电功率：1△Q=—一BV2=-0.5305MVarL1 2L1N同理可得变压器2（2台变压器并联运行）：R=1xAPSVN2x103=1.4715。；T2 2SNX=1xVS%VN2x10=34.715。T2 2SN=2(=2(AP+j4%0 100S)=(0.0372+j0.2)MVA线路2：R=rxR=rxl=0.21x140=29.4。l=0.42x140=58.8。产生的充电功率：AQL2-1BV2=-2.542L2NMVar变压器3：5NX103=0.1008QSN—VS%Vn2X10=5.51。4SNASOASO3线路3（双回路）:=4(AP+j~^S)=(0.716+j0.882)MVA100NTOC\o"1-5"\h\z1 1R = 2 r l = 2X 0.21 X 80 =8.4。X = 1 x l = 1x 0.42 x 80 =16.8 QL3 2 0 2产生的充电功率：1AQ =—-BV2=—2.899MVar变压器4（两台变压器并联）:1APV2R=—X—X103=1.886QT42SN1V%V2X=r―lx10=43.4QT42SNAS=2(AP+j10%S)=(0.0314+j0.16)MVA04 0 100N线路4：R=rl=0.21x70=14.7QX=xl=0.42x70=29.4Q产生的充电功率：AQ=-1BV2=-1.268MVarL4 2L4N变压器5(两台三绕组变压器并联)：RT(5-1)=RT(5-2)=RT(5-3)=4.116QV%=1[V%+V%-V%]=10.75S12S(1-2)S(1-3) S(2-3)1V%=_[V %+V %-V%]=-0.25S22S(1-2)S(2-3) S(1-3)1V%=一[V %+V %-V %]=6.75S32S(1-3)S(2-3) S(1-2)XT5-1ST——LX10=7.11QS7——lx10=—0.1653QS3——lx10=4.465。2(AP2(AP0+j-0100S)=(0.0264+j0.11)MVA线路5（双回路）:rl=1x0.21x30=3.15Q2l=-X0.42X30=6.3Q2产生的充电功率:AQ=-1BV2=-1.087MVarL5 2L5N变压器6(两台变压器并联)：R=1A.匕2X103=0.3422Q；T6 2SNX=1VS%VN2X10=7.185QT6 2SN=2(AP+0j1%S)=2(AP+0100N发电机A：P=4x15=60MWQ=4x15xtanp=45MVar发电机B：P=4x63=252MWQ=4x63xtan中=156MVar发电机C：P=3X12=36MWQ=3x12xtan中=27MVar发电机D：P=63x1=63MWQ=63x1xtan中=39MVar

发电机E：P=2x25=50MWQ=2x25xtan中=38MVar3.2潮流计算变压器1的功率损耗:AS=102(R+jX)=(0.0312+j0.7173)MVAT1 1102 T1 T1a点功率S=10MW+AS+AS+jAQ=10+0.0312+j0.7173+0.0157+j0.08-j0.5435=(10.0469+j0.2538)MVA线路1的功率损耗：AS=P+气(R+jQ)=(0.056+j0.1118)MVAL1 1102L1L2变压器T2的功率损耗：602+602+4521102(1.4715+j34.715)=(0.6840+j16.138)MVA所以b点功率为S=60-20+j45-0.6840-j16.138=(39.316+j28.862)MVAc点功率：S=S-25-S-jAQ-AS=(4.213+j29.04)MVA计算线路I的电压分布：VA=115V;变压器2中的电压降落：=15.27kV40x104715+45=15.27kVAV= T2 115V=115kV-AV=99.7kV线路I中的电压降落：10.1029x6072+(-0.2897x13.44)AV= =0.64kVL1 100

V=V-AV=99.36kV线路II中:1322+1561322+15621102(0.1008+j5.51)=(0.345+j19.02)MVA7.59kV132x0.1008+156x7.59kV115S=(92.684+j108.118)MVAV=115-7.59=107.41kV线路III的功率损耗和电压降落:ASL3ASL392.6842+108.11821102(8.4+j16.8)=(14.08+j28.16)MVA92.684x8.4+108.118x16.8 =22.56kV115V=107.41-22.56=84.85kVS=(113+j98.8)MVA线路IV的损耗:1132+98.821102(14.7+j29.4)=(27.37+j54.74)1102113x14.7+98.8x29.4115AV = =39.70kV115L4点f的功率和电压:S=S-AS =(85.36+j44.06)MVAAV=84.85-39.70=45.15kV变压器6的功率和损耗:152352(0.3422+j7.185)=(0.063+j1.329)MVA15x0.3422 =0.147kV35点g的功率及电压:S=15+AS+AS=(15.085+j1.43)MVAV=35+AV6=35.147kV线路5的功率损耗及电压降落:S+jAQ5=(15.085+j1.43)-j1.087=(15.085+j0.3431)MVAASL5ASL515.0852+0.3432352(3.15+j6.3)=(0.585+j1.117)MVA15.085x3.15+0.343x6.3 =1.42kV35点h的功率:S=S+jAQ+AS +jAQ=(15.67+j0.426)MVAV=V+AV=36.567kV变压器5-2的功率损耗和电压降落:35215.62+0.4262 (4.166-j1.653)=(0.828-j0.321)35215.6x4.166-0.426x1.653 =1.84kV35点m的功率和电压:=(16.498+j=(16.498+j1.06)MVA=38.407kV变压器5-3的功率损耗和电压降落:152+152+382102(4.166+j4.465)=(69.53+j74.37)MVA15x4.166+38x4.465 =23.216kV102点n的功率和电压:S=15+AS=(84.53+j112.37)MVAV=10+AV=33.216kV3.3MATLAB仿真（由陈志彪完成）相关的原始数据输入格式如下：1、 B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写。2、 对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号，将变压器的串联阻抗置于低压侧处理，第三列为支路的串列阻抗参数，第四列为支路的对地导纳参数，第五烈为含变压器支路的变压器的变比，第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，“0”为不含有变压器。3、 B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。4、 X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地参数。程序如下：closeallclearalln=input请输入节点数：n=');n1=input请输入支路数:n1=');isb=inputi青输入平衡节点号：isb=');pr=input请输入误差精度:pr=');B1=input(请输入支路参数:B1=');B2=input(请输入节点参数:B2=');X=input(节点号和对地参数:X=');Y=zeros(n);Times=1; %置迭代次数为初始值%创建节点导纳矩阵fori=1:n1ifB1(i,6)==0 %不含变压器的支路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);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));endendYOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%将OrgS、DetaS初始化%创建OrgS，用于存储初始功率参数h=0;j=0;fori=1:n %对PQ节点的处理ifi~=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))*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)==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))*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)));endendendOrgS%创建PVU用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;fori=1:nifB2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endendPVU%创建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);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=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);DetaS(2*h,1)=real(PVU(t,1))A2+imag(PVU(t,1))A2-real(B2(i,3))A2-imag(B2(i,3))A2;endendDetaS%创建I,用于存储节点电流参数i=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));endendI%创建论由(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;fori=1:n %对PQ节点的处理ifB2(i,6)==2h=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%非对角元素的处理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)=-Jacbi(2*h-1,2*k);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~=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)=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) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendJacbi%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压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%开始循环**********************************************************************whileabs(max(DetaU))>prOrgS=zeros(2*n-2,1); %!!!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0;j=0;fori=1:nifi〜=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))*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:nifi〜=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))*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)));endendendOrgS%创建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)==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))A2-real(B2(i,3))A2-imag(B2(i,3))A2;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));endendI%创建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(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;endendendendendk=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(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));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