电力系统潮流计算

1、精选优质文档-倾情为你奉上电力系统课程设计题 目: 电力系统潮流计算 院系名称： 电气工程学院 专业班级： 电气F1206班 学生姓名： 学 号： 指导教师： 张孝远 成绩：指导老师签名：日期： 目录原始资料21 概述42 潮流计算节点介绍42.1 变量的分类52.2 节点的分类53 计算方法简介63.1 牛顿拉夫逊法原理63.1.1 牛顿拉夫逊法概要63.1.2 牛顿法的框图及求解过程83.2 MATLAB简介94 潮流分布计算104.1 系统的一次接线图104.2 参数计算104.3 丰大及枯大下地潮流分布情况144.3.1 该地区变压器的有功潮流分布数据154.3.2 重、过载负荷元件统

2、计表175 设计心得17参考文献18附录:程序19原始资料一、系统接线图见附件1。二、系统中包含发电厂、变电站、及其间的联络线路。500kV变电站以外的系统以一个等值发电机代替。各元件的参数见附件2。设计任务1、手动画出该系统的电气一次接线图，建立实际网络和模拟网络之间的联系。2、根据已有资料，先手算出各元件的参数，后再用Matlab表格核算出各元件的参数。3、潮流计算1）对两种不同运行方式进行潮流计算，注意110kV电网开环运行。2）注意将电压调整到合理的范围110kV母线电压控制在106kV117kV之间;220kV母线电压控制在220 kV242kV之间。附件一：专心-专注-专业附件二：

3、1、变压器：两个220kV变电站均采用参数一致的三绕组变压器，具体参数如下。220kV变电站参数表高压侧绕组中压侧绕组低压侧绕组容量12012060（120）电压220110(121)10.5110kV及以下的变电站的变压器省略，即可将负荷直接挂在110kV母线上。而110kV升压变只计及以下参数。110kV变电站参数表序号变电站名容量X1X01A162B283C1205D636E207G31.58F12.5+31.59水电站130MW、2*20MVA0.82250 0.52500 10水电站272MW、3*31.5MVA0.34444 0.33333 11水电站324MW、31.5MVA1.

4、03333 0.33333 12水电站54*7.5MW、2*20MVA1.03333 0.33333 13水电站418MW、30MVA1.28333 1.05000 14火电厂50MW、31.5+40MVA0.48285 0.26250 2、线路：具体参数如下。220kV线路参数表序号线路名称导线牌号线路长度km1ML2x240102ML2x240103MH2x3005110kV线路参数表序号线路名称导线牌号线路长度km1水电站1水电站2150302AB95803BC9514水电站4C150655BD95636CD240607水电站2C240758ED24014.79水电站3D15010.51

5、0水电站5T节点705.511T节点D70512T节点G240513H F1503.914H D185715水电站2L240/213016LG2401017DF150418火电厂D24013、发电机各发电机的参数如下：XdXq装机容量功率因数水电站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变电站

6、丰大负荷按该站变电容量的50%估算，枯大负荷按该站变电容量的60%估算。两个220kV变电站的低压侧上各挂10MW的负荷，中压侧各挂20MW负荷。功率因素均为0.95。5、并联电容器两个220kV变电站的低压侧上均装设并联补偿。补偿总量按该站变电容量的20%装设，分组原则以每组电容器的容量不超过10MVar且经济性较好为准。1 概述潮流计算是电力系统最基本最常用的计算。根据系统给定的运行条件，网络接线及元件参数，通过潮流计算可以确定各母线的电压，包括电压的幅值和相角，各元件流过的功率，整个系统的功率损耗等一系列系统中的潮流数据。近几年，对潮流算法的研究仍然是如何改善传统的潮流算法，即高斯-塞德

7、尔法、牛顿法和快速解耦法。牛顿法，由于其在求解非线性潮流方程时采用的是逐次线性化的方法，为了进一步提高算法的收敛性和计算速度，人们考虑采用将泰勒级数的高阶项或非线性项也考虑进来，于是产生了二阶潮流算法。后来又提出了根据直角坐标形式的潮流方程是一个二次代数方程的特点，提出了采用直角坐标的保留非线性快速潮流算法。潮流计算在数学上是多元非线性方程组的求解问题，求解的方法有很多种，牛顿拉夫逊Newton-Raphson法是数学上解非线性方程组的有效方法，有较好的收敛性。将N-R法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性，稀疏性及节点编号顺序优划等技巧，使N-R法在收敛性，占用内存，

8、计算速度等方面的优点都超过了阻抗法总结为在电力系统运行方式和规划方案的研究中都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时为了实时监控电力系统的运行状态也需要进行大量而快速的潮流计算。因此潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时采用离线潮流计算在电力系统运行状态的实时监控中则采用在线潮流计算。2 潮流计算节点介绍常规的电力系统潮流计算中一般具有三种类型的节点：PQ、PV及平衡节点。一个节点有四个变量，即注入有功功率、注入无功功率，电压大小及相角。常规的潮流计算一般给定其中的二个变量：PQ节点（注入有功功率及

9、无功功率），PV节点（注入有功功率及电压的大小），平衡节点（电压的大小及相角）。2.1 变量的分类负荷消耗的有功、无功功率、电源发出的有功、无功功率、母线或节点的电压大小和相位、 、在这十二个变量中，负荷消耗的有功和无功功率无法控制，因它们取决于用户，它们就称为不可控变量或是扰动变量。电源发出的有功无功功率是可以控制的自变量，因此它们就称为控制变量。母线或节点电压的大小和相位角是受控制变量控制的因变量。其中, 、主要受、的控制, 、主要受、的控制。这四个变量就是简单系统的状态变量。为了保证系统的正常运行必须满足以下的约束条件：对控制变量 对没有电源的节点则为对状态变量的约束条件则是对某些状态变

10、量还有如下的约束条件2.2 节点的分类 第一类称PQ节点。等值负荷功率、和等值电源功率、是给定的，从而注入功率、是给定的，待求的则是节点电压的大小和相位角。属于这类节点的有按给定有功、无功率发电的发电厂母线和没有其他电源的变电所母线。 第二类称PV节点。等值负荷和等值电源的有功功率、是给定的，从而注入有功功率是给定的。等值负荷的无功功率和节点电压的大小 也是给定的。待求的则是等值电源的无功功率，从而注入无功功率和节点电压的相位角。有一定无功功率储备的发电厂和有一定无功功率电源的变电所母线都可以作为PV节点； 第三类平衡节点。潮流计算时一般只设一个平衡节点。等值负荷功率、是给定的，节点电压的大小

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

12、下边的方程来确定：式中等号右边的矩阵都是对于的值。这一矩阵称为雅可比（JACOBI）矩阵。按上述得到的修正向量后，得到如下关系这比更接近真实值。这一步在收敛到希望的值以前重复进行，一般要反复计算满足为预先规定的小正数，是第n次迭代的近似值。3.1.2 牛顿法的框图及求解过程1、用牛顿法计算潮流时，有以下的步骤：（1） 给这各节点电压初始值；（2） 将以上电压初始值代入公式，求修正方程的常数项向量 ；（3） 将电压初始值在带入上述公式，求出修正方程中系数矩阵的各元素。（4） 解修正方程式；（5） 修正各节点电压，；（6） 将，在带入方程式，求出；（7） 检验是否收敛，即（8） 如果收敛，迭代到此

13、结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果。如果不收敛，转回（2）进行下次迭代计算，直到收敛为止。2、程序框图如下3.2 MATLAB简介 MATLAB是用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度

14、上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB是一种交互式、面向对象的程序设计语言广泛应用于工业界与学术界主要用于矩阵运算同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。MATLAB程序设计语言结构完整且具有优良的移植性它的基本数据元素是不需要定义的数组。它可以高效率地解决工业计算问题特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言可以用类似数学公式的方式来编写算法大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是

15、编程上。目前电子计算机已广泛应用于电力系统的分析计算潮流计算是其基本应用软件之一。现有很多潮流计算方法。对潮流计算方法有五方面的要求（1）计算速度快（2）内存需要少（3）计算结果有良好的可靠性和可信（4）适应性好亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强。4 潮流分布计算4.1 系统的一次接线图图4.1 系统的一次连接图4.2 参数计算设定基准值，Ub=Uav.n，则各参数如下。（1）发电机的次暂态电抗：X=X*Sb/Sn，ZB=UB2/SN发电机参数 单位（MW）电厂装机容量枯水出力比例丰水出力比例丰大有功丰大无功枯大有功枯大无功短路X*''水电站1

16、300.20.721.000 13.015 6.000 3.718 0.493水电站2720.20.750.400 31.235 14.400 8.924 0.水电站3240.20.716.800 10.412 4.800 2.975 0.61625水电站5300.20.721.000 13.015 6.000 3.718 0.493华鑫电厂500.80.840.000 24.790 40.000 24.790 0.2958（2）110KV升压变压器的参数：电阻：R=PK*UN2/ 1000SN2;电抗：X=UK(%)*UN2/ 100 SN;电导：G=P0/ 1000 UN2;电纳：B=I0

17、(%)*SN/ 100 UN2;式中UN以KV为单位，SN以MVA为单位，P0、PK以KW为单位。110KV变压器参数 变电站名容量X1X0SbR*X*G*B*A16B28C120D63E20G31.5F12.5+31.5水电站130MW、2*20MVA0.82250.5251000. 0. 0. 0. 水电站272MW、3*31.5MVA0.344440.333331000. 0. 0. 0. 水电站324MW、31.5MVA1.033330.333331000. 0. 0. 0. 水电站54*7.5MW、2*20MVA1.033330.333331000. 0. 0. 0. 火电厂50MW

18、、31.5+40MVA0.482850.26251000. 0. 0. 0. 1000. 0. 0. 0. （3）220KV三绕组变压器的参数：电阻：R=PK*UN2/ 1000SN2;电抗：X=UK(%)*UN2/ 100 SN;电导：G=P0/ 1000 UN2;电纳：B=I0(%)*SN/ 100 UN2220kV变电站参数表SBR* X* G*B*高压侧绕组中压侧绕组低压侧绕组高压侧绕组（1）1000. 0. 0. 0. 容量120120120中压侧绕组（2）100-0. 0. 0. 0. 电压22012110.5低压侧绕组（3）1000. -0. 0. 0. （4）110KV线路参数

19、： r=/S； lg ； Xl*=x1*l*SB/U2av.n。110KV线路参数标么值序号线路名称导线牌号线路长度kmUB(KV)SB（MW)R*X*B*1水电站1水电站2150301151000. 0. 0. 2AB95801151000. 0. 0. 3BC9511151000. 0. 0. 4水电站4C150651151000. 0. 0. 5BD95631151000. 0. 0. 6CD240601151000. 0. 0. 7水电站2C240751151000. 0. 0. 8ED24014.71151000. 0. 0. 9水电站3D15010.51151000. 0. 0.

20、10水电站5T节点705.51151000. 0. 0. 11T节点D7051151000. 0. 0. 12T节点G24051151000. 0. 0. 13H F1503.91151000. 0. 0. 14H 0. 0. 15水电站2L2401301151000. 0. 0. 16LG240101151000. 0. 0. 17D 0. 0. 18火电厂D24011151000. 0. 0. 110KV变电站负荷参数 变电站总容量丰变站比枯站变比丰大有功丰大无功枯大有功枯大无功A160.50.68.0000 2.6295 9.6000

21、 3.1554 B280.50.614.0000 4.6016 16.8000 5.5219 C1200.50.660.0000 19.7210 72.0000 23.6653 D630.50.631.5000 10.3535 37.8000 12.4243 E200.50.610.0000 3.2868 12.0000 3.9442 F440.50.622.0000 7.2311 26.4000 8.6773 G31.50.50.615.7500 5.1768 18.9000 6.2121 各变电站补偿电容参数 站别运行方式ABCDEFGH低压侧L低压侧丰大34.609.805.901.92

22、.0枯大34.40005.901.92.04.3 丰大及枯大下地潮流分布情况电压是衡量电力系统电能质量的标准之一。电压过高或过低，都将对人身及其用电设备产生重大的影响。保证用户的电压接近额定值是电力系统调度的基本任务之一。当系统的电压偏离允许值时，电力系统必须应用电压调节技术调节系统电压的大小，使其维持在允许值范围内。本文经过手算形成了等值电路图，并编写好了程序得出节点电压标幺值，使其满足所要求的调整范围。我们首先对给定的程序输入部分作了简要的分析，程序开始需要我们确定输入节点数、支路数、平衡母线号、支路参数矩阵、节点参数矩阵。（1）为了保证整个系统潮流计算的完整性，我们把凡具有母线及发电机处

23、均选作节点，这样，我们确定电厂一母线上的发电机作为平衡节点，节点号为，其它机组作为PV节点，节点号为，其余节点均为PQ节点，节点号见等值电路图。（2）确定完节点及编号后，各条支路也相应确定了，我们对各支路参数进行了计算。根据所给实际电路图和题中的已知条件，有以下公式计算各输电线路的阻抗和对地支路电容的标幺值和变压器的阻抗标幺值。4.3.1 该地区变压器的有功潮流分布数据(1)该图为丰大潮流模型图： 丰大潮流模型图该运行方式的电压合理，负荷分配也均匀，但是有些线路的负载率偏低。比如水电厂2-L站的负载率仅仅为21.1% ，M-L站的线路的负载率也只是8.8%。(2)该图为枯大潮流模型图： 枯大潮

24、流模型图该运行方式的电压合理，负荷分配也均匀，有些线路的负载率偏低但是有些线路的负载率则偏高。比如水电厂2-L站的负载率仅仅为9.4%，M-L站的线路的负载率也只是 7.6% D-H站的线路却重载到84%。4.3.2 重、过载负荷元件统计表类型负载率实际是在功率额定容量建议线路0.83870.486.2限D站负荷到4MW以下5 设计心得通过这次课程设计,我发现自己有很多不足的地方,如基础知识掌握不牢固,很多知识点都忘记了,计算速度慢及准确性低,分析问题能力不够全面等等。同时，在设计的过程中遇到很多问题，如怎样使用WORD的工具,计算公式输入，画图等。明白了有些东西看起来很简单，但一旦做起来却需

25、要很多心思，要注意到很多细节问题。要做到能好好理解课本的内容，一定要认认真真做一次计算。因此，完成课程设计使我对课本的内容加深了理解。总体来说，这次的课程设计不单在专业基础方面反映了我的学习还要加倍努力，还在对一些软件的应用需要加强。计算在各种情况下的潮流分布，对于丰大和枯大情况下的潮流分布有了明确的认识。在本次课程设计过程当中，锻炼了自己实际操作分析能力，理论联系实际，对运行中的电力系统，通过潮流计算可以预知各种负荷变化和网络结构的改变会不会危及系统的安全，系统中所有母线的电压是否在允许的范围以内，系统中各种元件(线路、变压器等)是否会出现过负荷，以及可能出现过负荷时应事先采取哪些预防措施等

26、。同时采用MATLAB进行计算机的计算，在计算时采用特殊算法使得潮流计算的过程更快，效率更高。而采用计算机的运算应该是未来的一种趋势，所以我会学习一定的编辑语言如C,C+等，以提高运算准确性和快速性。总体而言，这次的课程设计对我们运用所学知识，发现、提出、分析和解决实际问题、锻炼实践能力的考察，使我们更清楚地知道不足之出，从而提高我们。参考文献1于永源 主编.电力系统分析 湖南师范大学出版社M.1992年7月2陈珩编.电力系统稳态分析 水利电力出版社M.1995年1月第二版3邱晓燕.刘天琪电力系统分析的计算机算法北京中国电力出版 2004李光琦.电力系统暂态分析M.北京：水利电力出版社，200

27、2.5?5陆敏政 主编.电力系统习题集 水利电力出版社M.1990年附录:程序%潮流计算MATLAB 粗略程序?%creat a new_datat=0;s=0;r=0;w=0;number=input('How many node are there=');% Convert Pq to a new arrayfor ii=1:number? if data(ii,4)=1? t=t+1;? for jj=1:14? new_data1(t,jj)=data(ii,jj);? end;? a(1,t)=ii;? s=s+1;? %record the number of th

28、e PQ node? end;end;%Convert pv to a new arrayfor ii=1:number? if data(ii,4)=2? t=t+1;? for jj=1:14? new_data1(t,jj)=data(ii,jj);? end;? a(1,t)=ii;? r=r+1;? %record the number of the PV node? end;end;%Convert set_v to a new arrayfor ii=1:number? if data(ii,4)=3? t=t+1;? for jj=1:14? new_data1(t,jj)=d

29、ata(ii,jj);? end;? a(1,t)=ii;? w=w+1;? end;end;%creat a new_data2x,y=size(data2)?for ii=1:x? for jj=1:2? for mm=1:number? if data2(ii,jj)=a(1,mm)? new_data2(ii,jj)=mm;? end;? end;? end;end;for ii=1:x? for jj=3:14? new_data2(ii,jj)=data2(ii,jj);? end;end;%creat a YY=zeros(number,number);YY=zeros(numb

30、er,number);yy=zeros(number,number);for ii=1:x? % for jj=1:14? iii=new_data2(ii,1);? jjj=new_data2(ii,2);? if new_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).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*

31、i;? Y(iii,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)=

32、-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).*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(ii,7)+new_data2(ii,6).*new_data2(ii,6)-new_data2(ii,7)./(new_data2(ii,

33、7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;? Y(jjj,iii)=-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).*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(ii

34、,7)+new_data2(ii,6).*new_data2(ii,6)-new_data2(ii,7)./(new_data2(ii,7).*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;for iii=1:number? Y(iii,iii)=0;end;%for ii=1:x? %? for jj=1:14? for? iii=1:number? for j

35、j=1:number? % if iii=jj? Y(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);? % end;? end;end;%creat B, Gfor ii=1:number? for jj=1:number? G(ii,jj)= real(Y(ii,jj);? B(ii,jj)= imag(Y(ii,jj);? end;end;%creat Initial_P Initial_Q Initial_Vfor ii=1:(s+r)? set_P(ii,1)=(new_data1(ii,9)-new_data1(ii,7)./100;end;fo

36、r ii=1:s;? set_Q(ii,1)=(new_data1(ii,10)-new_data1(ii,8)./100;end;for ii=1:r? set_V(ii,1)=new_data1(ii+s,12).*new_data1(ii+s,12);%try to modify for sike of correctingend;Initial_p_q_v=set_P;set_Q;set_V;disp(Initial_p_q_v);%creat Initial_e,Initial_ffor ii=1:number-1? e(ii,1)=1;? f(ii,1)=0.0;%change f

37、 to test used to be 1.0end;(12,1)=1.07;? %e(13,1)=1.09;% Start NEWTOWN CALULATIONfor try_time=1:25%Creat every node consume P Q and Un=s;m=r;for ii=1:(n+m)? sum1=0;? for jj=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).*e(jj,1);? end;? p(ii,1

38、)=sum1;end;for ii=1:n? sum2=0;? for jj=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);for ii=(n+1):(n+m)? u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);end;for ii=n+1:(n+

39、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);%creat Jacobiandisp(n);disp(m);for ii=1:(n+m)? for jj=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;?

40、qq=0; ? for num=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;? 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=')for

41、 ii=1:n? for jj=1:m+n? if (ii=jj)? QE(ii,jj)=B(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;? for num=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+

42、G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1? QE(ii,jj)=ss+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);for ii=n+1:n+m? for jj=1:(n+m)? if (ii=jj)? UE(ii,jj)=0;? UF(ii,jj)=0;? else? ss=0;

43、? qq=0;? for num=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; ? UF(ii,jj)=-2.*f(ii,1);? UE(ii,jj)=-2.*e(ii,1);? end;? end;end;for ii=(n+1):(n+m)? for jj=1:(n+m)? extra_UE(ii-n),jj)=UE(ii,jj);? extra_UF(ii-n),jj)=UF(ii,jj);? end;end;%

44、disp('extra_UE');%disp(extra_UE);%disp('extra_Uf');%disp(extra_UF);Jacobian=PF,PE;QF,QE;extra_UF,extra_UE;%disp('Jacobian=');%disp(Jacobian);%creat substract resultsubstract_result=Initial_p_q_v-sum;%disp('substract_result');%disp(substract_result);%calculate delta_f_

45、edelta_f_e=-inv(Jacobian)*substract_result;%disp(delta_f_e);for ii=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;if max(substract_result)<1e-4? break;end ;end;%disp('substract_result');%disp(substract_result);%disp('e=');%disp(e);%disp

46、('f=');%disp(f);for ii=1:number? uuu(ii,1)= e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);U_RESULT(ii,1)=sqrt(uuu(ii,1);end;for ii=1:number? for? jj=1:number? if ii=a(1,jj)? Old_Uresult(ii,1)=U_RESULT(jj,1)? end;end;end;for ii=1:number? Old_Uresult(ii,2)=ii;end;%disp('U_result');%disp(U_RESULT);