电力系统潮流计算,DOC

1、电力系统课程设计题目：电力系统潮流计算 院系名称：电气工程学院 专业班级：电气F1206班 学生姓名：学号：指导教师：张孝远2.成绩：目录料指导老师签名：.2日期：4卜算节点介绍4-变量的分类51概述5原始资2潮流2.1 TOC o 1-5 h z 3计算方法简介6 HYPERLINK l bookmark41 o Current Document 3.1牛顿一拉夫逊法原理63.1.1牛顿一拉夫逊法概要6 HYPERLINK l bookmark44 o Current Document 3.1.2牛顿法的框图及求解过程8 HYPERLINK l bookmark55 o Current Do

2、cument 3.2MATLAB 简介94潮流分布计算10104.1系统的一次接线图 HYPERLINK l bookmark61 o Current Document 4.2参数计算10 HYPERLINK l bookmark109 o Current Document 4.3丰大及枯大下地潮流分布情况144.3.1该地区变压器的有功潮流分布数据154.3.2重、过载负荷元件统计表175设计心得17参考文献18附录:程序19原始资料一、系统接线图见附件1。二、系统中包含发电厂、变电站、及其间的联络线路。500kV变电站以外的系统以一个等值发电机代替。各元件的参数见附件2。设计任务手动画出该

3、系统的电气一次接线图，建立实际网络和模拟网络之间的联系。1、2、根据已有资料，先手算出各元件的参数，后再用Matlab表格核算出各元件的参数。3、潮流计算1）对两种不同运行方式进行潮流计算，注意110kV电网开环运行。2）注意将电压调整到合理的范围110kV母线电压控制在106kV117kV之间;220kV母线电压控制在220kV242kV之间。附件一:500kV变电站附件二：1、变压器：两个220kV变电站均采用参数一致的三绕组变压器，具体参数 如下。220kV变电站参数表高压侧绕组中压侧绕组低压侧绕组容量12012060(120)电压220110(121)10.5110kV及以下的变电站的

4、变压器省略，即可将负荷直接挂在110kV母线上。而110kV升压变只计及以下参数。110kV变电站参数表序号变电 站名容量X1X01A162B283C1205D636E207G31.58F12.5+31.59水电站130MW、2*20MVA0.822500.5250010水电站272MW、3*31.5MVA0.344440.3333311水电站324MW、31.5MVA1.033330.3333312水电 站54*7.5MW、2*20MVA1.033330.3333313水电 站418MW、30MVA1.283331.0500014火电 厂50MW、31.5+40MVA0.482850.262

5、502、线路：具体参数如下。220kV线路参数表序号线路名称导线牌号线路长度km1ML2x240102ML2x240103MH2x3005110kV线路参数表序号线路名称导线牌号线路长度km1水电站1水电站2150302AB95803BC9514水电站4C150655BD95636CD240607水电站2C240758ED24014.79水电站3D15010.510水电站5T节点705.511T节点D70512T节点G240513HF1503.914HD185715水电站2L240/213016LG2401017DF150418火电厂D24013、发电机各发电机的参数如下:Xd Xq 装机 容

6、量功率 因数水电 站10.1740.174300.85水电站20.1740.174720.85水电站30.1740.174240.85水电 站50.1740.174300.85水电 站40.1740.174180.85华鑫电厂0.1740.17450.85出力情况：水力发电机丰大出力70%,枯大出力20%。火力发电机丰大出力80%,枯大出力80%。4、负荷各110kV变电站丰大负荷按该站变电容量的50%估算，枯大负荷按该站变电 容量的60%估算。两个220kV变电站的低压侧上各挂10MW的负荷，中压侧各挂20MW负荷。功率因素均为0.95。5、并联电容器两个220kV变电站的低压侧上均装设并联

7、补偿。补偿总量按该站变电容量 的20%装设，分组原则以每组电容器的容量不超过10MVar且经济性较好为准。1概述潮流计算是电力系统最基本最常用的计算。根据系统给定的运行条件，网 络接线及元件参数，通过潮流计算可以确定各母线的电压，包括电压的幅值和 相角，各元件流过的功率，整个系统的功率损耗等一系列系统中的潮流数据。近几年，对潮流算法的研究仍然是如何改善传统的潮流算法，即高斯-塞德 尔法、牛顿法和快速解耦法。牛顿法，由于其在求解非线性潮流方程时采用的 是逐次线性化的方法，为了进一步提高算法的收敛性和计算速度，人们考虑采 用将泰勒级数的高阶项或非线性项也考虑进来，于是产生了二阶潮流算法。后 来又提

8、出了根据直角坐标形式的潮流方程是一个二次代数方程的特点，提出了 采用直角坐标的保留非线性快速潮流算法。潮流计算在数学上是多元非线性方程组的求解问题，求解的方法有很多种， 牛顿一拉夫逊Newton-Raphson法是数学上解非线性方程组的有效方法，有较好 的收敛性。将N-R法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩 阵的对称性，稀疏性及节点编号顺序优划等技巧，使N-R法在收敛性，占用内 存，计算速度等方面的优点都超过了阻抗法总结为在电力系统运行方式和规划方案的研究中都需要进行潮流计算以比 较运行方式或规划供电方案的可行性、可靠性和经济性。同时为了实时监控电 力系统的运行状态也需要进行大

9、量而快速的潮流计算。因此潮流计算是电力系 统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系 统的运行方式时采用离线潮流计算在电力系统运行状态的实时监控中则采用在 线潮流计算。2潮流计算节点介绍常规的电力系统潮流计算中一般具有三种类型的节点：PQ、PV及平衡节点。 一个节点有四个变量，即注入有功功率、注入无功功率，电压大小及相角。常 规的潮流计算一般给定其中的二个变量：PQ节点（注入有功功率及无功功率）， PV节点（注入有功功率及电压的大小），平衡节点（电压的大小及相角）。2.1变量的分类负荷消耗的有功、无功功率一一匕、qli、P22、ql2电源发出的有功、无功功率匕、qgi

10、、pg2、Qg2母线或节点的电压大小和相位一一U 1、U2、8 1、5 2在这十二个变量中，负荷消耗的有功和无功功率无法控制，因它们取决于 用户，它们就称为不可控变量或是扰动变量。电源发出的有功无功功率是可以 控制的自变量，因此它们就称为控制变量。母线或节点电压的大小和相位角一 是受控制变量控制的因变量。其中，U、U主要受Q、Q的控制，5、5主 1 2G1G 21 2要受PG1、PG2的控制。这四个变量就是简单系统的状态变量。为了保证系统的正常运行必须满足以下的约束条件:对控制变量对没有电源的节点则为对状态变量匕的约束条件则是对某些状态变量5还有如下的约束条件i2.2节点的分类第一类称PQ节点

11、。等值负荷功率P、Q和等值电源功率P、Q是给 Li LiGi Gi定的，从而注入功率P、Q是给定的，待求的则是节点电压的大小气和相位角 5,。属于这类节点的有按给定有功、无功率发电的发电厂母线和没有其他电源 的变电所母线。第二类称pv节点。等值负荷和等值电源的有功功率pl匕.是给定的，从而注入有功功率P是给定的。等值负荷的无功功率q和节点电压的大小U也 ini是给定的。待求的则是等值电源的无功功率Q，从而注入无功功率Qt和节点电 压的相位角5。有一定无功功率储备的发电厂和有一定无功功率电源的变电所 母线都可以作为PV节点；第三类平衡节点。潮流计算时一般只设一个平衡节点。等值负荷功率Pl、 Ql

12、?是给定的，节点电压的大小和相位也是给定的。担负调整系统频率任务的发 电厂母线往往被选作为平衡节点。3计算方法简介3.1牛顿一拉夫逊法原理3.1.1牛顿一拉夫逊法概要首先对一般的牛顿一拉夫逊法作一简单的说明。已知一个变量X函数为：到此方程时，由适当的近似值x出发，根据：反复进行计算，当X满足适当的收敛条件就是上面方程的根。这样的方 法就是所谓的牛顿一拉夫逊法。这一方法还可以做下面的解释，设第n次迭代得到的解语真值之差，即X 的误差为时，则：把f （ X （n） +e ）在X（ n）附近对用泰勒级数展开上式省略去 2以后部分X（n）的误差可以近似由上式计算出来。比较两式，可以看出牛顿一拉夫逊法的

13、休整量和X（n）的误差的一次项相等。用同样的方法考虑，给出n个变量的n个方程：对其近似解X得修正量AX可以通过解下边的方程来确定： TOC o 1-5 h z 11式中等号右边的矩阵 f 都是对于X1，X；,., X：的值。这一矩阵称为雅可比n（JACOBI）矩阵。按上述得到的修正向量AX, AX，,AX后，得到如下关系12n这比X。X，,X更接近真实值。这一步在收敛到希望的值以前重复进行， 12n一般要反复计算满足为预先规定的小正数，Xnn+1是第n次迭代Xn的近似值。3.1.2牛顿法的框图及求解过程1、用牛顿法计算潮流时，有以下的步骤：（1）给这各节点电压初始值e（o）, f（o）；（2）

14、将以上电压初始值代入公式，求修正方程的常数项向量AP（0）, AQ （o） ,（AV 2 ）（o）；（3）将电压初始值在带入上述公式，求出修正方程中系数矩阵的各元素。（4）解修正方程式Ae（0）,Af （0）；(5)修正各节点电压e=e(0) + Ae(0)，f =f (0) + Af (0)；（6）将e，f在带入方程式，求出AP（i）,AQ（i）,（AV2）（1）；（7）检验是否收敛，即maxAp（k）,AQ（k） Ksii（8）如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率， 并打印输出结果。如果不收敛，转回（2）进行下次迭代计算，直到收敛为止。2、程序框图如下3.2MATLA

15、B 简介MATLAB是用于算法开发、数据可视化、数据分析以及数值计算的高级技术 计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。是由美国 mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技 计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的 建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、 工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方 案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编 辑模式，代表了当今国际科学计算软件的先进水平。MATLAB是一种交互式、

16、面向对象的程序设计语言广泛应用于工业界与学 术界主要用于矩阵运算同时在数值分析、自动控制模拟、数字信号处理、 动态分析、绘图等方面也具有强大的功能。MATLAB程序设计语言结构完整 且具有优良的移植性 它的基本数据元素 是不需要定义的数组。它可以高效率地解决工业计算问题 特别是关于矩阵和 矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言 可以用类似数学公式的方式来编写算法大大降低了程序所需的难度并节省了 时间，从而可把主要的精力集中在算法的构思而不是编程上。目前电子计算机已广泛应用于电力系统的分析计算潮流计算是其基本应用软件之一。现有很多潮流计算方法。对潮流计算

17、方法有五方面的要求（1）计算 速度快（2）内存需要少（3）计算结果有良好的可靠性和可信（4）适应性好 亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强。4潮流分布计算4.1系统的一次接线图图4.1系统的一次连接图4.2参数计算设定基准值SB 100MVA , Ub二Uav.n,则各参数如下。（1）发电机的次暂态电抗：X=X*Sb/Sn, ZB=Ub2/Sn发电机参数单位（MW）电厂装 机 容 量枯水 出力 比例丰水 出力 比例丰大有 功丰大无 功枯大有 功枯大无 功短路X*水电 站1300.20.721.00013.0156.0003.7180.493水电 站2720.20

18、.750.40031.23514.4008.9240.205417水电 站3240.20.716.80010.4124.8002.9750.61625水电 站5300.20.721.00013.0156.0003.7180.493华鑫 电厂500.80.840.00024.79040.00024.7900.2958（2）110KV升压变压器的参数:电阻：R=P *U 2/1000S2；KNN电抗：X=Uk（%）*Un2/100Sn；电导：G=P0/1000Un2；电纳：B=I0（%）*Sn/100Un2；式中U以KV为单位，S以MVA为单位，P、P以KW为单位。NN0 K110KV变压器参数变

19、电站 名容量X1X0SbR*X*G*B*A16B28C120D63E20G31.5F12.5+31.5水电站130MW、2*20MVA0.82250.5251000.0373640.5812120.0001990.001445水电站272MW、3*31.5MVA0.344440.333331000.0211990.3690230.0002800.001992水电站324MW、31.5MVA1.033330.333331000.0211990.3690230.0002800.001992水电站54*7.5MW、2*20MVA1.033330.333331000.0373640.5812120.0

20、001990.001445火电厂50MW、31.5+40MVA0.482850.26251000.0211990.3690230.0002800.0019921000.0138380.2906060.0003790.002529（3）220KV三绕组变压器的参数:电阻：R=P *U 2/1000S2；KNN电抗：X=Uk(%)*Un2/100Sn；电导：G=P0/1000UN2；电纳：B=I0(%)*Sn/100Un2220kV变电站参数表SBR*X*G*B*高压侧绕 组中压侧绕 组低压侧绕 组高压 侧绕 组 （1）1000.0029200.1166500.0016500.012000容量12

21、0120120中压 侧绕 组 （2）100-0.0021380.0716670.0004990.003630电压22012110.5低压侧绕组（3）1000.003580-0.0045750.0000040.000027（4） 110KV线路参数:r=p/S；二=0Sig； 一二5 Xl*=x1*l*SB/U2av.n。110KV线路参数标么值序号线路名称导线牌号线路长度 kmUB(KV)SB(MW)R*X*B*1水电 站1水电 站2150301151000.0476370.0907370.0115062AB95801151000.1996220.2522500.0297303BC951115

22、1000.0024950.0031530.0003724水电 站4C150651151000.1032140.1965970.0249295BD95631151000.1572020.1986470.0234126CD240601151000.0598870.1751230.0239647水电 站2C240751151000.0748580.2189040.0299558ED24014.71151000.0146720.0429050.0058719水电 站3D15010.51151000.0166730.0317580.00402710水电 站5T节 点 八、705.51151000.018

23、7150.0177580.00198611T节 点 八、D7051151000.0170130.0161440.00180512T节 点 八、G24051151000.0049910.0145940.00199713HF1503.91151000.0061930.0117960.00149614H0089980.0210660.00274015水电 站2L2401301151000.1297540.3794330.05192116LG240101151000.0099810.0291870.00399417D0063520.0120980.0

24、0153418火电 厂D24011151000.0009980.0029190.000399110KV变电站负荷参数变电站总容 量丰变站比枯站变比丰大有功丰大无 功枯大有 功枯大无功A160.50.68.00002.62959.60003.1554B280.50.614.00004.601616.80005.5219C1200.50.660.000019.721072.000023.6653D630.50.631.500010.353537.800012.4243E200.50.610.00003.286812.00003.9442F440.50.622.00007.231126.40008.

25、6773G31.50.50.615.75005.176818.90006.2121各变电站补偿电容参数、站别运行 方式、ABCDEFGH低压 侧L低压 侧丰大34.609.805.901.92.0枯大34.40005.901.92.04.3丰大及枯大下地潮流分布情况电压是衡量电力系统电能质量的标准之一。电压过高或过低，都将对人身 及其用电设备产生重大的影响。保证用户的电压接近额定值是电力系统调度的 基本任务之一。当系统的电压偏离允许值时，电力系统必须应用电压调节技术 调节系统电压的大小，使其维持在允许值范围内。本文经过手算形成了等值电 路图，并编写好了程序得出节点电压标幺值，使其满足所要求的调

26、整范围。我们首先对给定的程序输入部分作了简要的分析，程序开始需要我们确定 输入节点数、支路数、平衡母线号、支路参数矩阵、节点参数矩阵。（1）为了保证整个系统潮流计算的完整性，我们把凡具有母线及发电机处 均选作节点，这样，我们确定电厂一母线上的发电机作为平衡节点，节点号为 ，其它机组作为PV节点，节点号为，其余节点均为PQ节点，节点号 见等值电路图。（2）确定完节点及编号后，各条支路也相应确定了，我们对各支路参数进行 了计算。根据所给实际电路图和题中的已知条件，有以下公式计算各输电线路的阻抗和对地支路电容的标幺值和变压器的阻抗标幺值。4.3.1该地区变压器的有功潮流分布数据该图为丰大潮流模型图：

27、丰大潮流模型图该运行方式的电压合理，负荷分配也均匀，但是有些线路的负载率偏低。 比如水电厂2L站的负载率仅仅为21.1%，ML站的线路的负载率也只是 8.8%。该图为枯大潮流模型图：枯大潮流模型图该运行方式的电压合理，负荷分配也均匀，有些线路的负载率偏低但是有 些线路的负载率则偏高。比如水电厂2L站的负载率仅仅为9.4%，ML站的线路的负载率也只是7.6%D H站的线路却重载到84%。4.3.2重、过载负荷元件统计表类型负载率实际是在功率额定容量建议线路0.83870.486.2限D站负荷到4MW以下5设计心得通过这次课程设计，我发现自己有很多不足的地方，如基础知识掌握不牢固, 很多知识点都忘

28、记了，计算速度慢及准确性低，分析问题能力不够全面等等。同 时，在设计的过程中遇到很多问题，如怎样使用WORD的工具，计算公式输入， 画图等。明白了有些东西看起来很简单，但一旦做起来却需要很多心思，要注 意到很多细节问题。要做到能好好理解课本的内容，一定要认认真真做一次计 算。因此，完成课程设计使我对课本的内容加深了理解。总体来说，这次的课 程设计不单在专业基础方面反映了我的学习还要加倍努力，还在对一些软件的应用需要加强。计算在各种情况下的潮流分布，对于丰大和枯大情况下的潮流分布有了明 确的认识。在本次课程设计过程当中，锻炼了自己实际操作分析能力，理论联 系实际，对运行中的电力系统，通过潮流计算

29、可以预知各种负荷变化和网络结 构的改变会不会危及系统的安全，系统中所有母线的电压是否在允许的范围以 内，系统中各种元件（线路、变压器等）是否会出现过负荷，以及可能出现过负 荷时应事先采取哪些预防措施等。同时采用MATLAB进行计算机的计算，在计算 时采用特殊算法使得潮流计算的过程更快，效率更高。而采用计算机的运算应 该是未来的一种趋势，所以我会学习一定的编辑语言如C,C+等，以提高运算准 确性和快速性。总体而言，这次的课程设计对我们运用所学知识，发现、提出、分析和解 决实际问题、锻炼实践能力的考察，使我们更清楚地知道不足之出，从而提高 我们。参考文献于永源主编.电力系统分析湖南师范大学出版社M

30、.1992年7月陈珩编.电力系统稳态分析水利电力出版社M.1995年1月第二版邱晓燕.刘天琪 电力系统分析的计算机算法 北京 中国电力出版200李光琦.电力系统暂态分析M.北京：水利电力出版社，2002.5?陆敏政主编.电力系统习题集水利电力出版社M.1990年附录：程序%潮流计算MATLAB粗略程序？?？%creatanew_datat=0;s=0;r=0;w=0;number二input(Howmanynodearethere=)；%ConvertPqtoanewarrayforii=1:number?ifdata(ii,4)=1?t=t+1;?forjj=1:14?new_data1(t

31、,jj)=data(ii,jj);?end;?a(1,t)=ii;?s=s+1;?%recordthenumberofthePQnode?end;end;%Convertpvtoanewarrayforii=1:number?ifdata(ii,4)=2?t=t+1;?forjj=1:14?new_data1(t,jj)=data(ii,jj);?end;?a(1,t)=ii;?r=r+1;?%recordthenumberofthePVnode?end;end;%Convertset_vtoanewarrayforii=1:number?ifdata(ii,4)=3?t=t+1;?forjj

32、=1:14?new_data1(t,jj)=data(ii,jj);?end;?a(1,t)=ii;?w=w+1;?end;end;%creatanew_data2x,y=size(data2)?forii=1:x?forjj=1:2?formm=1:number?ifdata2(ii,jj)=a(1,mm)?new_data2(ii,jj)=mm;?end;?end;?end;end;forii=1:x?forjj=3:14?new_data2(ii,jj)=data2(ii,jj);?end;end;%creataYY=zeros(number,number);YY=zeros(numbe

33、r,number);yy=zeros(number,number);forii=1:x ?%forjj=1:14?iii=new_data2(ii,1);?jjj=new_data2(ii,2);?ifnew_data2(ii,5)=2?sub=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)-new_data2(ii,7)./(new_data2(ii,7).*n ew_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;?Y(iii

34、,jjj)=-sub./new_data2(ii,14);?YY(iii,jjj)=sub./new_data2(ii,14);?Y(jjj,iii)=-sub/new_data2(ii,14);?YY(jjj,iii)=sub./new_data2(ii,14);?yy(iii,jjj) = (1.-new_data2(ii,14)./(new_data2(ii,14).*new_data2(ii,14).*sub;?yy(jjj,iii) = (new_data2(ii,14)-1)./(new_data2(ii,14).*sub;?else?Y(iii,jjj)=-new_data2(i

35、i,6)./(new_data2(ii,7).*new_data2(+new_data2(ii,6).*new_data2(ii,6)+new_data2(ii,7)./(new_data2(.*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;?YY(iii,jjj)=new_data2(ii,6)./(new_data2(ii,7).*new_data2+new_data2(ii,6).*new_data2(ii,6)-new_data2(ii,7)./(new_data2.*new_data2(ii,7)+new_data2(ii,6)

36、.*new_data2(ii,6)*i;?Y(jjj,iii)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(+new_data2(ii,6).*new_data2(ii,6)+new_data2(ii,7)./(new_data2(.*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;?YY(jjj,iii)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(+new_data2(ii,6).*new_data2(ii,6)-new_data2(ii,7)

37、./(new_data2(.*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;?yy(iii,jjj)=new_data2(ii,8)./2.*i;?yy(jjj,iii)=new_data2(ii,8)./2.*i;?end;?%end;end;foriii=1:number?Y(iii,iii)=0;end;%forii=1:x?%?forjj=1:14?for?iii=1:number?forjj=1:number?%ifiii=jj?Y(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);?%end;

38、?end;end;%creatB,Gforii=1:number?forjj=1:number?G(ii,jj)=real(Y(ii,jj);?B(ii,jj)=imag(Y(ii,jj);?end;end;%creatInitial_PInitial_QInitial_Vforii=1:(s+r)?set_P(ii,1) = (new_data1(ii,9)-new_data1(ii,7)./100;end;forii=1:s;?set_Q(ii,1) = (new_data1(ii,10)-new_data1(ii,8)./100;end;forii=1:r?set_V(ii,1)=new

39、_data1(ii+s,12).*new_data1(ii+s,12);%trytomodifyforsikeofcorrectingend;Initial_p_q_v=set_P;set_Q;set_V;disp(Initial_p_q_v);%creatInitial_e,Initial_fforii=1:number-1?e(ii,1)=1;?f(ii,1)=0.0;%changeftotestusedtobe1.0end;?e(number,1)=new_data1(number,12);?f(number,1)=0;?%e(64,1)=0.88;%test118ieee?%f(64,

40、1)=0.39395826829394;?%f(14,1)=0;?%e(10,1)=1.045;?%e(11,1)=1.01;?%e(12,1)=1.07;?%e(13,1)=1.09;%StartNEWTOWNCALULATIONfortry_time=1:25%CreateverynodeconsumePQandUn=s;m=r;forii=1:(n+m)?sum1=0;?forjj=1:(n+m+1)?sum1=sum1+e(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1)+f(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*

41、e(jj,1);?end;?p(ii,1)=sum1;end;forii=1:n?sum2=0;?forjj=1:(n+m+1)?sum2=sum2+f(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1)-e(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1);?end;?q(ii,1)=sum2;end;disp(，q=，);disp(q);u=zeros(n+m),1);forii=(n+1):(n+m)?u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);end;forii=n+1:(

42、n+m)?extra_u(ii-n),1)=u(ii,1);end;disp(，extra_u=，);disp(extra_u);sum=p;q;extra_u;disp(sum)disp(s);disp(p);%creatJacobiandisp(n);disp(m);forii=1:(n+m)?forjj=1:(n+m)?if(ii=jj)?PF(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);?PE(ii,jj)=-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);?else?ss=0;?qq=0;?fornum=1:(n+m+

43、1)?ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);?qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);?end;?PF(ii,jj)=-ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1?PE(ii,jj)=-qq-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);%TEST+1?end;?end;end;copy=3.14159;disp(=copy=)forii=1:n?forjj=1:m+n?if(ii=jj)?QE(ii,jj)=B(

44、ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1 ?QF(ii,jj)=G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1?else?ss=0;?qq=0;?fornum=1:(n+m+1)?ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);?qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);?end;?QF(ii,jj)=-qq+G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1?QE(ii,jj)=ss

45、+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1 ?end;?end;end;%disp(QF);%disp(QF);%disp(QE);%disp(QE);UE=zeros(n+m),(n+m);UF=zeros(n+m),(n+m);forii=n+1:n+m?forjj=1:(n+m)?if(ii=jj)?UE(ii,jj)=0;?UF(ii,jj)=0;?else?ss=0;?qq=0;?fornum=1:(n+m+1)?ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);?qq=qq+G(ii,num).

46、*e(num,1)-B(ii,num).*f(num,1);?end;?UF(ii,jj)=-2.*f(ii,1);?UE(ii,jj)=-2.*e(ii,1);?end;?end;end;forii=(n+1):(n+m)?forjj=1:(n+m)?extra_UE(ii-n),jj)=UE(ii,jj);?extra_UF(ii-n),jj)=UF(ii,jj);?end;end;%disp(extra_UE);%disp(extra_UE);%disp(extra_Uf);%disp(extra_UF);Jacobian=PF,PE;QF,QE;extra_UF,extra_UE;%d

47、isp(Jacobian=)；%disp(Jacobian);%creatsubstractresultsubstract_result=Initial_p_q_v-sum;%disp(substract_result);%disp(substract_result);%calculatedelta_f_edelta_f_e=-inv(Jacobian)*substract_result;%disp(delta_f_e);forii=1:number-1;?f(ii,1)=f(ii,1)+delta_f_e(ii,1);?e(ii,1)=e(ii,1)+delta_f_e(ii+number-1,1);end;ifmax(substract_result)1e-4?break;end;end;%disp(substract_result);%disp(substract_result);%disp(e=)；%disp(e);%disp(f=)；%disp(f);forii=1:number?uuu(ii,1)=e(ii