课 程 设 计 任 务 书

1、 山东交通学院电力系统分析课程设计报告书 院(部)别 信息科学与电气工程学院 班 级 电气114 学 号 110812436 姓 名 赵文兵 指导教师 栗玉霞 时 间 2014.06.9-2013.06.13 课 程 设 计 任 务 书题 目 电力系统分析课程设计 题目C 学 院 信息科学与电气工程学院 专 业 电气工程及其自动化 班 级 电气114（电力方向） 学生姓名 赵文兵 学 号 110812436 6 月 9 日至 6 月 13 日 共 1 周 6 月 16 日至 6 月 20 日 共 1 周指导教师(签字) 负责人(签字) 年 月 日一、 设计内容及要求1.设计目的 (1). 掌握

2、电力系统潮流计算的基本原理； (2). 掌握并能熟练运用一门计算机语言（MATLAB语言或FORTRAN或C语言或C+语言）； (3). 采用计算机语言对潮流计算进行计算机编程计算。2 .设计要求 (1). 程序源代码； (2). 给定题目的输入，输出文件； ( 3). 程序说明； (4). 给定系统的程序计算过程； 3 .设计内容 (1).根据电力系统网络推导电力网络数学模型，写出节点导纳矩阵； (2).赋予各节点电压变量（直角坐标系形式）初值后，求解不平衡量； (3).形成雅可比矩阵； (4).求解修正量后，重新修改初值，从2开始重新循环计算； ( 5).求解的电压变量达到所要求的精度时，

3、再计算各支路功率分布、功率损耗和平衡节 点功率; (6).上机编程调试；连调； ( 7).计算分析给定系统潮流分析并与手工计算结果作比较分析。 (8).准备计算机演示答辩，书写该课程设计说明书.二、 设计原始资料题目C：试用牛顿-拉夫逊潮流计算方法求解教材P88图4-10网络的潮流分布。求解精度为10e-6。求解思路 此电力系统是一个5节点，4支路的电力网络。其中包含3个PQ节点，一个PV节点，和一个平衡节点。综合比较牛顿拉夫逊法（直角坐标、极坐标）、PQ分解法等多种求解方法的特点，最后确定采用牛顿拉夫逊法（极坐标）。因为此方法所需解的方程组最少。三、设计完成后提交的文件和图表1计算说明书部分

4、附录有设计报告及手算结果2图纸部分：（1）设计计算程序以及程序注解（2）计算详细结果可见附录四、进程安排1.第一天上午：选题，查资料，制定设计方案；2.第一天下午第三天下午：手算完成潮流计算的要求；3.第四天上午第五天上午：编程完成潮流计算，并对照手算结果，分析误差4.第五天下午：答辩，交设计报告五、主要参考资料1.电力系统分析（第三版） 于永源主编，中国电力出版社，2007年2.电力系统分析何仰赞 温增银编著，华中科技大学出版社，2002年版；3.电力系统分析韩桢祥主编，浙江大学出版社，2001年版；4.电力系统稳态分析陈珩 编，水利电力出版社； 成绩评定表指导教师成绩答辩小组成绩总评成绩

5、摘要潮流计算是电力系统最基本最常用的计算。根据系统给定的运行条件，网络接线及元件参数，通过潮流计算可以确定各母线的电压，包括电压的幅值和相角，各元件流过的功率，整个系统的功率损耗等一系列数据。牛顿拉夫逊Newton-Raphson法是数学上解非线性方程组的有效方法，有较好的收敛性。将N-R法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性，稀疏性及节点编号顺序优划等技巧，使N-R法在收敛性，占用内存，计算速度等方面的优点都超过了阻抗法。本文运用MATLAB软件进行潮流计算，对给定题目进行分析计算，并应用了电力系统潮流计算仿真软件DDTR与利用程序计算的结果进行比较，使计算的结果更加

6、准确。利用成形的程序对系统中出现的各种情况，例如负荷的变化以及线路上所发生的变化进行计算，并对母线上不满足范围的数据进行调控，最终进行合理的系统潮流。目录1潮流计算11.1潮流计算概述11.2潮流计算的要求21.3潮流计算的优势21.4潮流计算的用途32 MATLAB简介42.1MATLAB概述43牛顿-拉夫逊法概述53.1牛顿-拉夫逊基本原理53.2直角牛顿-拉夫逊法潮流计算求解过程64设计过程124.1题目C124.2牛顿拉夫逊法程序流程图244.3设计程序25心得体会28参考文献29MATLAB概述目前电子计算机已广泛应用于电力系统的分析计算，潮流计算是其基本应用之一。现有很多潮流计算方

7、法。对潮流计算方法有五方面的要求：（1）计算速度快；（2）内存需要少；（3）计算结果有良好的可靠性和可信性；（4）适应性好，亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强；（5）简单。MATLAB是一种交互式、面向对象的程序设计语言，广泛应用于工业界与学术界，主要用于矩阵运算，同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。MATLAB程序设计语言结构完整，且具有优良的移植性，它的基本数据元素是不需要定义的数组。它可以高效率地解决工业计算问题，特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言

8、，可以用类似数学公式的方式来编写算法，大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。另外，MATLAB提供了一种特殊的工具：工具箱（TOOLBOXES）.这些工具箱主要包括：信号处理（SIGNALPROCESSING）、控制系统（CONTROLSYSTEMS）、神经网络（NEURALNETWORKS）、模糊逻辑(FUZZYLOGIC)、小波(WAVELETS)和模拟（SIMULATION）等等。不同领域、不同层次的用户通过相应工具的学习和应用，可以方便地进行计算、分析及设计工作。MATLAB设计中，原始数据的填写格式是很关键的一个环节，它与程序使用的方便

9、性和灵活性有着直接的关系。原始数据输入格式的设计，主要应从使用的角度出发，原则是简单明了，便于修改。牛顿-拉夫逊基本原理潮流计算的目标是求取电力系统在给定运行状态的计算。即节点电压和功率分布，用以检查系统各元件是否过负荷。各点电压是否满足要求，功率的分布和分配是否合理以及功率损耗等。对现有电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和暂态稳定分析都是以潮流计算为基础。潮流计算结果可用如电力系统稳态研究，安全估计或最优潮流等对潮流计算的模型和方法有直接影响。实际电力系统的潮流技术那主要采用牛顿-拉夫逊法。牛顿-拉夫逊法(简称牛顿法)在数学上是求解非线性代数方程式的有效方

10、法。其要点是把非线性方程式的求解过程变成反复地对相应的线性方程式进行求解的过程。即通常所称的逐次线性化过程。对于非线性代数方程组： 即 (3-1-1)在待求量x的某一个初始估计值附近，将上式展开成泰勒级数并略去二阶及以上的高阶项，得到如下的经线性化的方程组： (3-1-2)上式称之为牛顿法的修正方程式。由此可以求得第一次迭代的修正量 (3-1-3)将和相加，得到变量的第一次改进值。接着就从出发，重复上述计算过程。因此从一定的初值出发，应用牛顿法求解的迭代格式为： (3-1-4) (3-1-5)上两式中：是函数对于变量x的一阶偏导数矩阵，即雅可比矩阵J;k为迭代次数。有上式可见，牛顿法的核心便是

11、反复形式并求解修正方程式。牛顿法当初始估计值和方程的精确解足够接近时，收敛速度非常快，具有平方收敛特性。牛顿潮流算法突出的优点是收敛速度快，若选择到一个较好的初值，算法将具有平方收敛特性，一般迭代45次便可以收敛到一个非常精确的解。而且其迭代次数与所计算网络的规模基本无关。牛顿法也具有良好的收敛可靠性，对于对以节点导纳矩阵为基础的高斯法呈病态的系统，牛顿法也能可靠收敛。牛顿法所需的内存量及每次迭代所需时间均较高斯法多。牛顿法的可靠收敛取决于有一个良好的启动初值。如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行的节点上。对于正常运行的系统，各节点电压一般均在额定值附近，偏移不会太大，并

12、且各节点间的相位角差也不大，所以对各节点可以采用统一的电压初值(也称为平直电压)，如假定： 或 (3-1-6) 这样一般能得到满意的结果。但若系统因无功紧张或其它原因导致电压质量很差或有重载线路而节点间角差很大时，仍用上述初始电压就有可能出现问题。解决这个问题的办法可以用高斯法迭代12次，以此迭代结果作为牛顿法的初值。也可以先用直流法潮流求解一次以求得一个较好的角度初值，然后转入牛顿法迭代。直角牛顿-拉夫逊法潮流计算求解过程以下讨论的是用直角坐标形式的牛顿拉夫逊法潮流的求解过程。当采用直角坐标时，潮流问题的待求量为各节点电压的实部和虚部两个分量由于平衡节点的电压向量是给定的，因此待求两共需要2

13、(n-1)个方程式。事实上，除了平衡节点的功率方程式在迭代过程中没有约束作用以外，其余每个节点都可以列出两个方程式。对PQ节点来说，是给定的，因而可以写出 （3-2-1）对PV节点来说，给定量是，因此可以列出 （3-2-2）求解过程大致可以分为以下步骤：（1）形成节点导纳矩阵Y（2）将各节点电压设初值U，（3）将节点初值代入式(2-2-1)或式(2-2-2)，求出修正方程式的常数项向量（4）将节点电压初值代入求式，求出雅可比矩阵元素（5）求解修正方程，求修正向量（6）求取节点电压的新值（7）检查是否收敛，如不收敛，则以各节点电压的新值作为初值自第3步重新开始进行狭义次迭代，否则转入下一步（8）

14、计算支路功率分布，PV节点无功功率和平衡节点柱入功率。以直角坐标系形式表示. 迭代推算式 采用直角坐标时,节点电压相量及复数导纳可表示为: (3-2-3)将以上二关系式代入上式中,展开并分开实部和虚部;假定系统中的第1,2,m号为PQ节点,第m+1,m+2,n-1为PV节点,根据节点性质的不同,得到如下迭代推算式: 对于PQ节点 (3-2-4)对于PV节点 (3-2-5)对于平衡节点 平衡节点只设一个,电压为已知,不参见迭代,其电压为: (3-2-6). 修正方程式迭代式共包括2(n-1)个方程.选定电压初值及变量修正量符号,代入方程并按泰勒级数展开,略去二次方程及以后各项,得到一组线性方程组

15、或线性化了的方程组，常称修正方程组: (3-2-7) (3-2-8).雅可比矩阵各元素的算式式(3-2-8)中, 雅可比矩阵中的各元素可通过对式(3-2-4)和(3-2-5)进行偏导而求得.当时, 雅可比矩阵中非对角元素为 (3-2-9)当时,雅可比矩阵中对角元素为: (3-2-10)由式(2-2-9)和(2-2-10)看出,雅可比矩阵的特点:矩阵中各元素是节点电压的函数,在迭代过程中,这些元素随着节点电压的变化而变化;导纳矩阵中的某些非对角元素为零时,雅可比矩阵中对应的元素也是为零.若,则必有;雅可比矩阵不是对称矩阵;雅可比矩阵各元素的表示如下式(2-2-11): 牛顿拉夫逊法程序流程图潮流

16、计算的计算机算法程序：%本程序 的功能是用牛顿拉夫逊法进行潮流计算% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路首段对地电纳% 5、支路的变比；6、支路首端处于K侧为1，1侧为0不含变压器；% 7、支路末段对地电纳% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；% 3为PV节点；clear;n=4;%input(请输入节点数：n=);nl=4;%input(请输入支路数：nl=);isb=1;%input(请输入平衡母线节

17、点号：isb=);pr=0.00001;%input(请输入误差精度：pr=);B1=1 2 0.08+0.40i 0.33333i 1 0 0.02941i; 2 3 0.10+0.40i 0.02941i 1 0 0.03448i; 1 3 0.12+0.50i 0.33333i 1 0 0.03448i; 3 4 0.30i 0 0.9 1 0B2=0 0 1.05 0 0 1; 0 0.55+0.13i 1 0 0 2; 0 0.3+0.18i 1 0 0 2; 0.5 0 1.1 1.1 0 3 ;%input(请输入各节点参数形成的矩阵： B2=);Y=zeros(n);e=zer

18、os(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 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,7);%对角元K侧 Y(p,p)=Y(p,p)+1./B1(i,3

19、)+B1(i,4);%对角元1侧 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);%PV节点电压给定模值 endfor i=1:n%给定各节点注入功率 S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SL B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量 end%=P=real(S);Q=imag(S); %分解出各节点注入的有

20、功和无功功率ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; %迭代次数ICT1、a；不满足收敛要求的节点数IT2while IT2=0 % N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列 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); %(Gij*ej-Bij*fj) D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1); %(Gij*fj+Bij*ej) end P1=C(i)*e(i)+f(i

21、)*D(i); %节点功率P计算ei(Gij*ej-Bij*fj)+fi(Gij*fj+Bij*ej) Q1=C(i)*f(i)-e(i)*D(i); %节点功率Q计算fi(Gij*ej-Bij*fj)-ei(Gij*fj+Bij*ej)%求i节点有功和无功功率P,Q的计算值 V2=e(i)2+f(i)2;%电压模平方 %= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 = if B2(i,6)=3%非PV节点 DP=P(i)-P1;%节点有功功率差 DQ=Q(i)-Q1; %节点无功功率差 %= 以上为除平衡节点外其它节点的功率计算 =%= 求取Jacobi矩阵 = for j1=1

22、:n if j1=isb&j1=i%非平衡节点&非对角元 X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de X3=X2; % X2=dp/df X3=dQ/de X4=-X1; % X1=dP/de X4=dQ/df p=2*i-1;q=2*j1-1; J(p,q)=X3;J(p,N)=DQ;m=p+1; % X3=dQ/de J(p,N)=DQ节点无功功率差 J(m,q)=X1;J(m,N)=DP;q=q+1; % X1=dP/de J(m,N)=DP节点有功功率差

23、J(p,q)=X4;J(m,q)=X2; % X4=dQ/df X2=dp/df elseif j1=i&j1=isb%非平衡节点&对角元 X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); % dP/de X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dP/df X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/de X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i); % dQ/df p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ; %扩展列Q m=p+1; J(m,q)=X1;q=

24、q+1;J(p,q)=X4;J(m,N)=DP; %扩展列P J(m,q)=X2; end end else %= 下面是针对PV节点来求取Jacobi矩阵的元素 = DP=P(i)-P1;% PV节点有功误差 DV=V(i)2-V2;% PV节点电压误差 for j1=1:n if j1=isb&j1=i%非平衡节点&非对角元 X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de X2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df X5=0;X6=0; p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; % PV节点电压

25、误差 m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; % PV节点有功误差 J(m,q)=X2; elseif j1=i&j1=isb %非平衡节点&对角元 X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); % dP/de X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dP/df X5=-2*e(i); X6=-2*f(i); p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; % PV节点电压误差 m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

26、% PV节点有功误差 J(m,q)=X2; end end end end end%= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 = for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点） k1=k+1;N1=N; % N=N0+1 即 N=2*n+1扩展列P、Q 或 U for k2=k1:N1% 从k+1列的Jacobi元素到扩展列的P、Q 或 U J(k,k2)=J(k,k2)./J(k,k); % 用K行K列对角元素去除K行K列后的非对角元素进行规格化 end J(k,k)=1; % 对角元规格化K行K列对角元素赋1 %= 回代运算 = if k=3

27、 % 不是第三行 k 3 k4=k-1; for k3=3:k4% 用k3行从第三行开始到当前行的前一行k4行消去 for k2=k1:N1 % k3行后各行上三角元素 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0） end %用当前行K2列元素减去当前行k列元素乘以第k行K2列元素 J(k3,k)=0; %当前行第k列元素已消为0 end if k=N0 %若已到最后一行 break; end %= 前代运算 = for k3=k1:N0 % 从k+1行到2*n最后一行 for k2=k1:N1 % 从k+1列到扩展列消去k+1行后各行

28、下三角元素 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算 end %用当前行K2列元素减去当前行k列元素乘以第k行K2列元素 J(k3,k)=0; %当前行第k列元素已消为0 end else %是第三行k=3 %= 第三行k=3的前代运算 = for k3=k1:N0 %从第四行到2n行（最后一行） for k2=k1:N1 %从第四列到2n+1列（即扩展列） J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0） end %用当前行K2列元素减去当前行3列元素乘以第三行K2列元素 J(k3,k)=0; %当前

29、行第3列元素已消为0 end end end%=上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵= for k=3:2:N0-1 L=(k+1)./2; e(L)=e(L)-J(k,N); %修改节点电压实部 k1=k+1; f(L)=f(L)-J(k1,N); %修改节点电压虚部 end %-修改节点电压- for k=3:N0 DET=abs(J(k,N); if DET=pr %电压偏差量是否满足要求 IT2=IT2+1; %不满足要求的节点数加1 end end ICT2(a)=IT2; %不满足要求的节点数 ICT1=ICT1+1; %迭代次数end%用高斯消去法解w=-

30、J*Vdisp(迭代次数：);disp(ICT1);disp(没有达到精度要求的个数：);disp(ICT2);for k=1:n V(k)=sqrt(e(k)2+f(k)2); %计算各节点电压的模值 sida(k)=atan(f(k)./e(k)*180./pi; %计算各节点电压的角度 E(k)=e(k)+f(k)*j; %将各节点电压用复数表示end%= 计算各输出量 =disp(各节点的实际电压标幺值E为(节点号从小到大排列)：);disp(E); %显示各节点的实际电压标幺值E用复数表示disp(-);disp(各节点的电压大小V为(节点号从小到大排列)：);disp(V); %显

31、示各节点的电压大小V的模值disp(-);disp(各节点的电压相角sida为(节点号从小到大排列)：);disp(sida); %显示各节点的电压相角for p=1:n C(p)=0; for q=1:n C(p)=C(p)+conj(Y(p,q)*conj(E(q); %计算各节点的注入电流的共轭值 end S(p)=E(p)*C(p); %计算各节点的功率 S = 电压 X 注入电流的共轭值enddisp(各节点的功率S为(节点号从小到大排列)：);disp(S); %显示各节点的注入功率disp(-);disp(各条支路的首端功率Si为(顺序同您输入B1时一致)：);for i=1:n

32、l p=B1(i,1);q=B1(i,2); if B1(i,6)=0 Si(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)+(conj(E(p)*B1(i,5). -conj(E(q)*conj(1./(B1(i,3)*B1(i,5); Siz(i)=Si(p,q); else Si(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)+(conj(E(p)./B1(i,5). -conj(E(q)*conj(1./(B1(i,3)*B1(i,5); Siz(i)=Si(p,q); end disp(Si(p,q); SSi(p,q)=Si(p,q); ZF

33、=S(,num2str(p),num2str(q),)=,num2str(SSi(p,q); disp(ZF); disp(-);enddisp(各条支路的末端功率Sj为(顺序同您输入B1时一致)：);for i=1:nl p=B1(i,1);q=B1(i,2); if B1(i,6)=0 Sj(q,p)=E(q)*(conj(E(q)*conj(B1(i,7)+(conj(E(q)./B1(i,5). -conj(E(p)*conj(1./(B1(i,3)*B1(i,5); Sjy(i)=Sj(q,p); else Sj(q,p)=E(q)*(conj(E(q)*conj(B1(i,7)+(

34、conj(E(q)*B1(i,5). -conj(E(p)*conj(1./(B1(i,3)*B1(i,5); Sjy(i)=Sj(q,p); end disp(Sj(q,p); SSj(q,p)=Sj(q,p); ZF=S(,num2str(q),num2str(p),)=,num2str(SSj(q,p); disp(ZF); disp(-);enddisp(各条支路的功率损耗DS为(顺序同您输入B1时一致)：);for i=1:nl p=B1(i,1);q=B1(i,2); DS(i)=Si(p,q)+Sj(q,p); disp(DS(i); DDS(i)=DS(i); ZF=DS(,n

35、um2str(p),num2str(q),)=,num2str(DDS(i); disp(ZF); disp(-);endfigure(1);subplot(1,2,1);plot(V);xlabel(节点号);ylabel(电压标幺值);grid on;subplot(1,2,2);plot(sida);xlabel(节点号);ylabel(电压角度);grid on;figure(2);subplot(2,2,1);P=real(S);Q=imag(S);bar(P);xlabel(节点号);ylabel(节点注入有功);grid on;subplot(2,2,2);bar(Q);xlab

36、el(节点号);ylabel(节点注入无功);grid on;subplot(2,2,3);P1=real(Siz);Q1=imag(Siz);bar(P1);xlabel(支路号);ylabel(支路首端注入有功);grid on;subplot(2,2,4);bar(Q1);xlabel(支路号);ylabel(支路首端注入无功);grid on;运算结果和数据图片：B1 = 1.0000 2.0000 0.0800 + 0.4000i 0 + 0.3333i 1.0000 0 0 + 0.0294i 2.0000 3.0000 0.1000 + 0.4000i 0 + 0.0294i 1.0000 0 0 + 0.0345i 1.0000 3.0000 0.1200 + 0.5000i 0 + 0.3333i 1.0000 0 0 + 0.0345i 3.0000 4.0000 0 + 0.3000i 0 0.9000 1.0000 0 B2 = 0 0 1