电力系统分析课设.doc

电力系统分析课程设计基于Matlab的电力系统潮流计算专 业 电气工程及其自动化 班 级 电气1404 姓 名 梁智星 学 号 14039404 同 组 成 员 王少军 刘晗 指 导 教 师 李季 完 成 时 间 2016 。12 .30 目录绪论11 电力系统潮流计算要求21.1潮流计算具体分析22 潮流计算的简介32.1潮流计算的发展32。2潮流计算节点的类型33 潮流计算牛顿拉夫逊法53。1牛顿拉夫逊法简介53.2牛顿拉夫逊法原理53.3牛顿拉夫逊法潮流程序原理总框图74 牛顿拉夫逊法潮流计算实例84.1 实例电力系统84。2原始数据的输入125 运行结果126 心得体会15附 录：16参考文献：21 绪论 牛顿拉夫逊法（简称牛顿法）潮流计算是目前使用最为广泛、效果最好的一种潮流计算方法，这种方法把非线性方程式的求解方法过程变成反复对相应的线性方程式的求解过程，即主次线性化的过程提出一种基于牛顿-拉夫逊法的潮流计算方法，该方法改进了牛顿拉夫逊法中修正方程式的建立.将牛顿拉夫逊法中的雅可比矩阵进行了简化，并将修正量的求解分步进行,从而在计算中节约了存储单元的使用量，减少了每次迭代的计算量.。根据牛顿法本身的特点及电力系统本身的运行特性，我们可对直角坐标的牛顿法潮流计算作一些改进以求减少雅可比矩阵元素的计算，并改变计算过程从而提高计算速度，降低对计算机贮存容量的要求。本文首先对直角坐标牛顿拉夫逊计算正方程式作一简单介绍,然后根据电力系统运行特点提出假设,从而简化雅可比矩阵元素的计算，再根据改进法的特点提出新的计算框图，然后对给出的计算结果进行比较分析.潮流计算在数学上属于多元非线性方程组的求解问题,在用迭代法求解时涉及到大量的矢量和矩阵运算，因而寻求一个能高效地处理矩阵运算的语言，将给潮流计算程序的编制带来巨大的方便.而于二十世纪八十年代出现的科学与工程计算语言MATLAB,因其以复数矩阵为基本运算单元等特点，使它在潮流程序计算中具有独特的优势。总之在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算.因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。11 电力系统潮流计算要求 1.1潮流计算具体分析 (1)高斯赛德尔法潮流计算 潮流计算高斯赛德尔迭代法(Gauss一Seidel method)是求解电力系统潮流的方法。潮流计算高 斯赛德尔迭代法又分导纳矩阵迭代法和阻抗矩阵迭代法两种。前者是以节点导纳矩阵为基础建立的赛德尔迭代格式;后者是以节点阻扰矩阵为基础建立的赛德尔迭代格式。高斯赛德尔迭代法这是数学上求解线性或非线性方程组的一种常用的迭代方法。本实验通过对电力网数学模型形成的计算机程序的编制与调试,获得形成电力网数学模型：高斯-赛德尔法的计算机程序，使数学模型能够由计算机自行形成，即根据已知的电力网的接线图及各支路参数由计算程序运行形成该电力网的节点导纳矩阵和各节点电压、功率。通过实验教学加深学生对高斯-赛德尔法概念的理解，学会运用数学知识建立电力系统的数学模型，掌握数学模型的形成过程及其特点,熟悉各种常用应用软件，熟悉硬件设备的使用方法,加强编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。 (2）牛顿拉夫逊法潮流计算牛顿-拉夫逊法(简称牛顿法）潮流计算是目前使用最为广泛、效果最好的一种潮流计算方法，这种方法把非线性方程式的求解方法过程变成反复对相应的线性方程式的求解过程，即主次线性化的过程提出一种基于牛顿拉夫逊法的潮流计算方法,该方法改进了牛顿-拉夫逊法中修正方程式的建立.将牛顿拉夫逊法中的雅可比矩阵进行了简化,并将修正量的求解分步进行，从而在计算中节约了存储单元的使用量,减少了每次迭代的计算量.根据牛顿法本身的特点及电力系统本身的运行特性，我们可对直角坐标的牛顿法潮流计算作一些改进以求减少雅可比矩阵元素的计算,并改变计算过程从而提高计算速度，降低对计算机贮存容量的要求。 （3）PQ分解法潮流计算 P-Q分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思想是：把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据,把有功功率和无功功率迭代分开来进行。2 潮流计算的简介2.1潮流计算的发展 利用电子计算机进行潮流计算从20世纪50年代中期就已经开始。此后，潮流计算曾采用了各种不同的方法，这些方法的发展主要是围绕着对潮流计算的一些基本要求进行的。对潮流计算的要求可以归纳为下面几点： (1)算法的可靠性或收敛性 （2）计算速度和内存占用量 （3）计算的方便性和灵活性 2.2潮流计算节点的类型主要目的是由这些已知量去求电力系统内的各种电气量。所以，根据电力系统中各节点性质的不同，很自然地把节点分成三类：(1) PQ节点 对这类节点，等值负荷功率GiP、LiQ和等值电源功GiP、GiQ是给定的，从而注入功率iP、iQ是给定的，待求的则是节点电压的大小iU和相位角id。属于这一类节点的有按给定有功无功功率发电的发电厂母线和没有其他电源的变电所母线。 (2) PV节点 对这类节点，等值负荷和等值电源的有功功率LiP、GiP是给定的，从而注入有功功率iP是给定的。等值负荷的无功功率LiQ和节点电压的大小iU也是给定的。待求的则是等值电源的无功功率GiQ,从而注入无功功率iQ和节点电压的相位角id。有一定无功功率储备的发电厂和一定无功功率电源的变电所母线都可选作为PV节点.(3） 平衡节点 潮流计算时，一般只设一个平衡节点。对这节点，等值负荷功率是给定的，节点电压的大小SU和相位角Sd也是给定的，如给SU=1。0、Sd=0。待求的则是等值电源功率GsP、GsQ，从而注入功率sP、sQ。担负调整系统频率任务的发电厂母线往往被选作为平衡节点。例如，为提高计算的收敛性。可以选择出线数多或者靠近电网中心的发电厂母线作平衡节点. 进行计算时，平衡节点是不可少的;PQ节点是大量的；PV节点较少,甚至可能有。2。3节点功率方程 节点电压向量可以表示为极坐标的形式，也可以表示为直角坐标的形式,与此相对应，在潮流计算中节点功率方程也有两种形式.节点功率可表示为： 式2.3.1 如果上式中电压向量表示为极坐标的形式：导纳矩阵中元素表示为： 因此：（i=1，2，3n） 式2。3.2 又由 则可以得到: 式2.3.3式中：为两个节点的相位差。按上式的实部和虚部展开得： 式 2。3。4 这就是功率的极坐标方程式把上式中各节点的电压向量表示为直角坐标： 式2.3.5 式2.3.5带入式： 式2。3.6即可得到： 式2。3.7 式2。4。8式中 式2。4。9这就是功率的直角坐标方程式。3 潮流计算牛顿-拉夫逊法3。1牛顿拉夫逊法简介在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性和经济性。此外,在进行电力系统静态及暂态稳定计算时,要利用潮流计算的结果作为其计算的基础;一些故障分析以及优化计算也需要有相应的潮流计算作配合；潮流计算往往成为上述计算程序的一个重要组成部分.以上这些,主要是在系统规划设计及运行方式安排中的应用,属于离线计算范畴.3。2牛顿拉夫逊法原理 牛顿拉夫逊迭代法是常用的解非线性方程组的方法，也是当前广泛采用的计算潮流的方法，其标准模式如下. 设有非线性方程组 式3.2.1 设近似解与精确解分别相差则如下的关式应该成立： 式3.2。2上式中任何一式都可按泰勒级数展开，由此可得： 式3。2.3带入这些偏导数表示式时的计算所得 式3.2。4可简写为： 式3。2。5 运用这种方法计算时，ix的初值要选择比较接近他们的精确解，否则迭代过程可能不收敛。与运用高斯赛德尔法时不同，运用牛顿法拉夫逊法时，可以直接用以求解功率方程： 式3。2。6而为此需将有关值代入 式3。2。7并将实数部分和虚数部分分列 式3.2。8此外，由于系统中还有电压大小给定的PV节点，还应补充一组方程式 式3.2。9 牛顿法的核心便是反复形成并求解修正方程。牛顿法当初始估计值和方程的精确解足够接近时，收敛速度非常快。3.3牛顿拉夫逊法潮流程序原理总框图图3。2.14 牛顿拉夫逊法潮流计算实例 4。1 实例电力系统发电厂一发电厂二电变所1变电所2变电所3变电所4图4。1。1 电力系统电所1、2、3、4低压母线的电压等级分别为：10kv、 10kv、 10kv 、10kv；变电所负荷分别为：60Mw 、50Mw 、60Mw 、60Mw。功率因数,变电所分别配有两台容量为75WVA的变压器,短路损耗414KW,短路电压UK =16.7。发电厂和变电所之间的输电线路，单位长度的电阻是0.17，单位长度的电抗为0。402,单位长度的电纳为2.7810-6S。发电厂一总装机容量为400MW，发电厂二总装机容量为145MW。我们首先对给定的程序输入部分作了简要的分析，程序开始需要我们确定输入节点数、支路数、平衡母线号、支路参数矩阵、节点参数矩阵.为了保证整个系统潮流计算的完整性,我们把凡具有母线处均选作节点，这样,共有10条母线即选定10个节点，我们确定发电厂一母线为平衡节点，节点号为，发电厂二母线为PV节点，节点号为,其余节点均为PQ节点，节点号见等值电路图。 确定完节点及编号后，各条支路也相应确定了，网络中总计有10条支路，我们对各支路参数进行了计算。根据所给实际电路图和题中的已知条件，有以下公式计算各输电线路的阻抗和对地支路电容的标幺值和变压器的阻抗标幺值。选择电压基准值为 和功率基准值，所以 。 计算各线路参数：支路阻抗 支路对地电容 支路60km输电线路： 支路50km输电线路：支路70km输电线路：支路60km同杆架设双回线：支路70km同杆架设双回线：支路70km铜杆架设双回线： 计算各变器参数：75MVA的变压器：两台变压器并联，得阻抗 假设我们所采用的变压器有5个抽头,电压调节范围为22。5%， 对应的分接头开始时设变压器高压侧主分接头，设四个变电所变压器的非标准变比，降压变压器5个分接头的非标准变比如下，以备调压时选用变压器非标准变比按照上面数据输入后到得到各个接点电压标幺值，具体检验过程如下：可得出变压器二次侧的实际电压:形成节点参数矩阵时,还需要我们计算出各节点所接发电机功率,节点负荷功率（可根据负荷有功功率，发电机有功功率及已知的功率因数计算）变电所1 变电所2 变电所3 变电所4 平衡节点为节点，所设节点电压的初始值为，给定值为1.05PV节点为节点，给定电压值为1.05，在进行潮流计算前先规定三个矩阵 V矩阵：各节点电压初始值、节点分类号（1为平衡节点，2为PQ节点，3为PV节点）、节点电压给定值，节点所接的无功补偿设备容量故其V矩阵为 b矩阵：支路始端号、支路末端号、支路阻抗、线路电容、支路变比故b矩阵为 S矩阵：节点的注入功率为正，输出功率为负故S矩阵为4。2原始数据的输入 为了让程序更有通用性，以及方便初学用户使用,近几年，MATLAB编程大多都用到GUI功能。本次设计也不例外，用到了简单的人机对话界面作为原始数据输入界面。 在这次设计中，用到了一个Excel表格作为原始数据输入界面，通过该界面，用户不用依照矩阵的形式，将一连串的数据输入,而是，按照图表的提示，在图表中填入要求的电力系统节点、支路的参数，节点、支路个数以及要求精度即可。数据输入界面如图53所示. 在此,对数据的输入有以下几点说明： 在节点信息里,节点电压为迭代计算时所设的初值。 在节点信息里，节点类型一栏中，1表示平衡节点，2表示PQ节点，3表示PV节点. KT一栏要求输入的是变压器的变比，非标准变比变压器，KT=k（k1）,标准变压器 KT=1，若该线路无变压器，KT=0。 输入变压器电阻、电抗时,如无特殊说明，均采用归算到低压侧的数值，再进行计算. 本条程序默认节点数,支路数均在在100以下，可以解决绝大部分电力系统的潮流问题，若遇到超大系统,可对程序做稍加调整，仍然适用。5 运行结果运行结果nl = 10isb = 1pr = 1.0000e-05b = 1。0000 + 0。0000i 2。0000 + 0.0000i 0.0190 + 0。0460i 0。0000 + 0.0880i 1。0000 + 0.0000i 1。0000 + 0.0000i 3。0000 + 0。0000i 0.0220 + 0。0530i 0.0000 + 0.1020i 1.0000 + 0.0000i 2.0000 + 0。0000i 3。0000 + 0。0000i 0.0160 + 0.0380i 0.0000 + 0。0740i 1。0000 + 0.0000i 3.0000 + 0.0000i 4.0000 + 0.0000i 0.0096 + 0。0230i 0。0000 + 0.1760i 1.0000 + 0。0000i 4。0000 + 0。0000i 5.0000 + 0。0000i 0。0110 + 0.0270i 0。0000 + 0.2060i 1。0000 + 0。0000i 5。0000 + 0。0000i 6.0000 + 0。0000i 0.0110 + 0。0270i 0.0000 + 0。2060i 1。0000 + 0。0000i 7。0000 + 0.0000i 2.0000 + 0.0000i 0。0037 + 0.1115i 0.0000 + 0。0000i 1.0500 + 0.0000i 8.0000 + 0。0000i 3。0000 + 0。0000i 0.0037 + 0。1115i 0.0000 + 0。0000i 1。0500 + 0。0000i 9。0000 + 0.0000i 4。0000 + 0.0000i 0。0037 + 0。1115i 0。0000 + 0.0000i 1。0500 + 0.0000i 10。0000 + 0。0000i 5。0000 + 0。0000i 0。0037 + 0.1115i 0.0000 + 0.0000i 1。0500 + 0。0000iV = 1.0500 1。0000 1。0500 0 1.0000 2。0000 0 0 1.0000 2.0000 0 0 1。0000 2.0000 0 0 1。0000 2。0000 0 0 1.0500 3。0000 1。0500 0 1。0000 2.0000 0 0 1。0000 2.0000 0 0 1.0000 2.0000 0 0 1.0000 2。0000 0 0导纳矩阵Y： 1 至 6 列 14。3514 -34.5706i 7。6706 +18.5709i -6。6808 +16.0947i 0.0000 + 0.0000i 0。0000 + 0.0000i 0。0000 + 0.0000i 7。6706 +18。5709i 17。3796 -49。8015i -9.4118 +22。3529i 0。0000 + 0.0000i 0。0000 + 0。0000i 0。0000 + 0。0000i -6。6808 +16。0947i -9。4118 +22。3529i 31.8448 84。2579i -15.4550 +37。0275i 0.0000 + 0.0000i 0.0000 + 0。0000i 0.0000 + 0。0000i 0.0000 + 0.0000i -15。4550 +37。0275i 28.6934 77.5599i 12。9412 +31。7647i 0。0000 + 0.0000i 0.0000 + 0.0000i 0。0000 + 0.0000i 0.0000 + 0。0000i -12.9412 +31.7647i 26.1796 72.2822i -12.9412 +31.7647i 0。0000 + 0。0000i 0。0000 + 0.0000i 0。0000 + 0。0000i 0。0000 + 0.0000i 12。9412 +31.7647i 12。9412 31。6617i 0。0000 + 0.0000i -0.2831 + 8.5321i 0。0000 + 0.0000i 0。0000 + 0。0000i 0。0000 + 0。0000i 0。0000 + 0。0000i 0.0000 + 0。0000i 0.0000 + 0.0000i -0。2831 + 8.5321i 0。0000 + 0.0000i 0。0000 + 0.0000i 0。0000 + 0.0000i 0。0000 + 0。0000i 0.0000 + 0.0000i 0.0000 + 0。0000i -0.2831 + 8.5321i 0.0000 + 0.0000i 0.0000 + 0.0000i 0。0000 + 0.0000i 0.0000 + 0。0000i 0.0000 + 0。0000i 0.0000 + 0。0000i -0.2831 + 8。5321i 0.0000 + 0。0000i 7 至 10 列 0。0000 + 0。0000i 0。0000 + 0。0000i 0.0000 + 0.0000i 0。0000 + 0。0000i 0.2831 + 8。5321i 0。0000 + 0。0000i 0。0000 + 0。0000i 0.0000 + 0。0000i 0。0000 + 0.0000i -0.2831 + 8。5321i 0。0000 + 0.0000i 0。0000 + 0。0000i 0。0000 + 0。0000i 0.0000 + 0。0000i -0。2831 + 8。5321i 0。0000 + 0.0000i 0.0000 + 0.0000i 0。0000 + 0。0000i 0。0000 + 0.0000i 0.2831 + 8.5321i 0。0000 + 0。0000i 0.0000 + 0.0000i 0。0000 + 0.0000i 0。0000 + 0.0000i 0。2696 8.1258i 0.0000 + 0。0000i 0。0000 + 0.0000i 0。0000 + 0。0000i 0.0000 + 0.0000i 0.2696 8。1258i 0。0000 + 0。0000i 0.0000 + 0.0000i 0。0000 + 0。0000i 0.0000 + 0.0000i 0。2696 8.1258i 0。0000 + 0。0000i 0.0000 + 0.0000i 0.0000 + 0。0000i 0。0000 + 0。0000i 0。2696 8。1258i迭代次数= 3每次不满足个数= 17 18 13 0各节点电压大小= 1。0500 1。0230 1.0228 1。0195 1.0283 1。0500 1.0248 1。0336 1.0209 1.0307各节点电压相角= 0 0.8605 0.6340 -0.2327 0。9810 2。8984 -4。6234 3。7428 4.0231 2。7408节点电压= 1 至 6 列 1。0500 + 0.0000i 1.0229 0.0154i 1.0228 - 0。0113i 1.0194 0.0041i 1。0282 + 0.0176i 1.0487 + 0。0531i 7 至 10 列 1。0215 0.0826i 1.0314 - 0。0675i 1.0184 - 0.0716i 1。0295 - 0.0493i各节点功率= 1 至 6 列 0。9003 + 0.6813i 0.0000 + 0。0000i 0。0000 + 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i 1。4500 + 0.1613i 7 至 10 列 -0.6000 - 0。3720i -0。5000 0。3100i -0.6000 - 0。3720i 0。6000 - 0。3720i各支路首端功率= 0。5182 + 0.3570i 0。3822 + 0。3243i 0.0912 + 0。0032i 0。2163 + 0.1503i 0.8193 0.0993i -1。4283 0.3304i -0.6000 0。3720i 0.5000 0.3100i 0。6000 - 0。3720i 0.6000 0。3720i各支路末端功率= -0.5107 - 0.4335i -0.3764 - 0。4199i 0.0914 0.0803i 0.2173 - 0。3315i 0。8264 0.0993i 1。4500 + 0。1613i0.6019 + 0.4303i 0.5013 + 0。3498i 0。6020 + 0。4308i 0。6019 + 0。4297i各支路功率损耗= 0.0075 0.0765i 0。0058 0。0956i 0.0002 - 0。0771i 0.0010 0。1812i 0.0071 - 0.1985i 0.0217 0.1691i 0.0019 + 0.0583i 0。0013 + 0。0398i 0.0020 + 0.0588i网络总损耗= 0。0503 0。5834i6 心得体会 由于电力系统分析我学的还可以，所以在潮流计算理论方面没有太大题，这里并没有投入多少精力。但是程序设计对于我确实一大难题。拿到题目后,我就开始学习MATAB。说实话MATAB语言简易通俗，很容易理解各种指令的意思。但是对于没有编过程序的我来说，程序的行文逻辑很是棘手。通过多看程序，一句一句的注释,终于在元旦之前理解了整个潮流计算程序的脉络,有了一些进展。前段时间一直在调试程序,每次调试,都要从键盘输入电力网络的信息,诸如节点数、支路数、各节点信息、各支路信息等，一次次的输入很是麻烦。总之，通过这次课设设计，是我学会了如何更好的学习，如何更好的查资料，如何筛选好的信息.附 录：clear;clc；n=10；nl=10isb=1pr=0.00001b=1 2 0.0190+0.0460i 0。0880i 1。0000；1 3 0.0220+0.0530i 0.1020i 1;2 3 0。0160+0。0380i 0。0740i 1；3 4 0。0096+0。0230i 0.1760i 1；4 5 0.0110+0。0270i 0。2060i 1；5 6 0.0110+0.0270i 0.2060i 1；7 2 0。0037+0.1115i 0 1。05;8 3 0.0037+0.1115i 0 1.05；9 4 0.0037+0。1115i 0 1。05;10 5 0。0037+0.1115i 0 1.05V=1.05 1 1.05 0；1 2 0 0；1 2 0 0；1 2 0 0；1 2 0 0;1。05 3 1.05 0；1 2 0 0；1 2 0 0；1 2 0 0；1 2 0 0S=0;0;0；0；0;1。4500；-0.60000。3720i；0.5000-0。3100i；0。6000-0.3720i;-0.6000-0.3720i;%各节点的注入功率S=0；0；0；0；0;1.4500;-0.6000-0.3720i;0.5000-0.3100i;0。6000-0.3720i；-0.6000-0.3720i；各节点的注入功率w1=zeros（2n2,1）；P=real(S);Q=imag(S）；e=zeros（1，n）；f=zeros(1，n)；E=zeros（1，n）；Y=zeros（n);for i=1：nl%导纳矩阵生成 p=b（i，1);q=b（i，2)； Y(p，q)=Y(p,q)1./（b(i,3）b(i，5)）； Y（q,p）=Y（p，q）; Y(q，q）=Y(q,q)+1。/b（i,3)+b（i,4)。/2； Y（p,p)=Y（p，p)+1./(b(i,3）*b(i,5）2）+b（i,4)./2；enddisp（导纳矩阵Y：)disp(Y）U=zeros(1，n）；G=real（Y）；B=imag(Y);for i=1：ne（i)=real(V（i,1)）；f（i)=imag(V(i，1）)；U(i）=V（i,3）；B（i,i)=B（i,i）+V(i,3）;endT=0;co=0;d=0;while T=0 A=0；co=co+1;for i=2：n 生成雅可比矩阵和功率修正量 for j=2:n x=0；x1=0； if V(i，2）=2 for r=1:n x=x+（e（i)*(G（i，r）e(r）-B(i，r）f（r）+f（i）*（G（i，r)*f(r）+B(i，r)e(r）； x1=x1+（f(i)（G（i,r)e（r)-B(i，r)*f（r）)-e（i)（G(i，r)*f(r）+B（i,r）*e（r）; end w(2*i-1)=P(i)-x； w（2*i）=Q（i)x1； else if V(i，2)=3 for r=1:n x=x+（e（i)*(G(i，r）e(r)-B(i，r)f(r)）+f(i）(G（i，r)*f（r)+B（i,r)e(r）); end w(2i-1)=P（i）-x; w(2i)=U(i）2(e(i)2+f(i)2）; end end h=0;h1=0； if V（i,2）=2 if i=j for r=1:n if r=i continue end h=h+(G（i,r)f（r)+B(i，r)*e(r)； h1=h1+（G（i，r)*e（r）B（i,r）*f(r)）; end J(2*i-1，2*j1）=2*G（i，i）f（i）+h； J(2*i1,2j)=2G（i,i)e（i)+h1; J（2i，2*j-1）=-2*B(i，i)*f（i)+h1； J（2*i，2*j）=-2B(i,i)e（i）h； else J(2i-1，2*j-1）=-B(i，j)*e(i)+G（i，j）*f(i)； J（2*i-1，2j)=G(i,j）e(i)+B(i,j）*f(i）; J(2*i,2*j1）=G(i,j）e（i)-B(i，j）*f(i）; J（2i，2j）=B（i,j）*e（i)+G(i,j)f（i)； end else if V（i,2）=3 if i=j for r=1:n if r=i continue end h=h+(G（i,r）f（r）+B(i,r）*e（r)）； h1=h1+(G（i,r)e（r）-B（i,r）*f（r)）； end J（2*i1,2j1)=2*G(i,i)f（i)+h; J（2i1，2j)=2G(i，i）e（i)+h1; J(2*i，2j-1)=2*f（i)； J(2i,2j)=2e（i)； else J(2i1，2j1）=B（i,j）*e（i）+G（i,j）f(i）； J（2*i-1,2j）=G(i，j)*e（i）+B（i,j）*f（i)； J（2*i，2*j-1）=0； J（2*i，2*j）=0； end end end endend%disp（J)disp(w)for i=3:2n高斯消去法求电压修正量 for j=3：2n J1(i2,j2）=J(i,j); endendfor i=3：2*n w1（i-2）=w（i)；endu=zeros(2*n2，1)；N=2*n2；for k=1:N m=0;