电力系统分析课程设计-潮流计算.doc

课程设计 电力系统分析课 程 设 计成绩评定表设计课题 ： 潮流计算 学院名称 ： 电气工程学院 专业班级 ： 组 别 ： 学生姓名 ： 学 号 ： 指导教师 ： 设计时间 ： 2016.12.262016.12.31 指导教师意见：成绩: 签名： 年 月 日电力系统分析课 程 设 计课程设计名称： 潮流计算 专 业 班 级 ： 组 别 ： 学 生 姓 名 ： 学 号 ： 指 导 教 师 ： 课程设计时间： 2016.12.262016.12.31 摘 要电力系统潮流计算是电力系统最基本的计算，也是最重要的计算。所谓潮流计算，就是已知电网的接线方式与参数及运行条件，计算电力系统稳态运行各母线电压、个支路电流与功率及网损。对于正在运行的电力系统，通过潮流计算可以判断电网母线电压、支路电流和功率是否越限，如果有越限，就应采取措施，调整运行方式。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置定整计算、电力系统故障计算和稳定计算等提供原始数据。在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。关键词：电力系统潮流计算牛顿-拉夫逊法目 录1 设计要求41.1 原始资料41.2 设计任务41.3 该地区电力系统描述42 潮流计算节点介绍53 计算方法简介63.1 牛顿拉夫逊法的一般原理63.2 用牛顿法计算潮流时，有以下的步骤83.3 系统程序框图94 潮流分布计算114.1 系统的一次接线图114.2 参数计算114.3 丰大及枯大下地潮流分布情况155 心得体会18参考文献19附件一20附件二21附件三 程序211 设计要求1.1 原始资料（1）系统接线图见附件1。（2）系统中包含发电厂、变电站、及其间的联络线路。500kV变电站以外的系统以一个等值发电机代替。各元件的参数见附件2。1.2 设计任务（1）手动画出该系统的电气一次接线图，建立实际网络和模拟网络之间的联系。（2）根据已有资料，先手算出各元件的参数，后再用Matlab表格核算出各元件的参数。（3）潮流计算1）对两种不同运行方式进行潮流计算，注意110kV电网开环运行。2）注意将电压调整到合理的范围110kV母线电压控制在106kV117kV之间;220kV母线电压控制在220 kV242kV之间。1.3 该地区电力系统描述该地区包括该系统有2个220kV变电站，220kV线路总长25km；7个110kV变电站，110kV线路总长635.6km。该地区主要是110KV的变电站先与220KV的电站连接后，再通过500KV的M变电站再并入电网。拥有电源装机206MW，其中水电装机156MW，火电装机50M。2 潮流计算节点介绍常规的电力系统潮流计算中一般具有三种类型的节点：PQ、PV及平衡节点。一个节点有四个变量，即注入有功功率、注入无功功率，电压大小及相角。常规的潮流计算一般给定其中的二个变量：PQ节点（注入有功功率及无功功率），PV节点（注入有功功率及电压的大小），平衡节点（电压的大小及相角）。1、变量的分类：负荷消耗的有功、无功功率、电源发出的有功、无功功率、母线或节点的电压大小和相位、 、在这十二个变量中，负荷消耗的有功和无功功率无法控制，因它们取决于用户，它们就称为不可控变量或是扰动变量。电源发出的有功无功功率是可以控制的自变量，因此它们就称为控制变量。母线或节点电压的大小和相位角是受控制变量控制的因变量。其中, U1、U2主要受QG1、QG2的控制, 1、2主要受PG1、PG2的控制。这四个变量就是简单系统的状态变量。为了保证系统的正常运行必须满足以下的约束条件：对控制变量 对没有电源的节点则为对状态变量的约束条件则是 对某些状态变量还有如下的约束条件 2、节点的分类： 第一类称PQ节点。等值负荷功率、和等值电源功率、是给定的，从而注入功率、是给定的，待求的则是节点电压的大小和相位角。属于这类节点的有按给定有功、无功率发电的发电厂母线和没有其他电源的变电所母线。 第二类称PV节点。等值负荷和等值电源的有功功率、是给定的，从而注入有功功率是给定的。等值负荷的无功功率和节点电压的大小 也是给定的。待求的则是等值电源的无功功率，从而注入无功功率和节点电压的相位角。有一定无功功率储备的发电厂和有一定无功功率电源的变电所母线都可以作为PV节点； 第三类平衡节点。潮流计算时一般只设一个平衡节点。等值负荷功率、是给定的，节点电压的大小和相位也是给定的。担负调整系统频率任务的发电厂母线往往被选作为平衡节点。3 计算方法简介3.1 牛顿拉夫逊法的一般原理已知一个变量X函数为： ,到此方程时，由适当的近似值出发，根据：反复进行计算，当满足适当的收敛条件就是上面方程的根。这样的方法就是所谓的牛顿拉夫逊法。这一方法还可以做下面的解释，设第次迭代得到的解语真值之差，即的误差为时，则： 把在附近对用泰勒级数展开 上式省略去以后部分得： 的误差可以近似由上式计算出来。 比较两式，可以看出牛顿拉夫逊法的休整量和的误差的一次项相等。 用同样的方法考虑，给出个变量的个方程： 对其近似解得修正量可以通过解下边的方程来确定：式中等号右边的矩阵都是对于的值。这一矩阵称为雅可比（JACOBI）矩阵。按上述得到的修正向量后，得到 这比更接近真实值。这一步在收敛到希望的值以前重复进行，一般要反复计算满足为预先规定的小正数，是第n次迭代的近似值。3.2 用牛顿法计算潮流时，有以下的步骤输入线路，电气元件参数，形成节点导纳矩阵 。 给这各节点电压初始值 。 将以上电压初始值代入式（438a）式（438c）或式（445c）、（445a），求出修正方程式中的不平衡量 。 将各节点电压的初值代入式（4-41a）、式（4-41b）或式（4-49a）式（4-49d）,求修正方程式的系数矩阵雅克比矩阵的各个元素。 解修正方程式，求各节点电压的变化量，即修正量。 计算各节点电压的新值，即修正后值运用各节点电压的新值自第三步开始进入下一次迭代。计算平衡节点功率和线路功率。其中，平衡节点功率为 线路功率为 从而，线路上损耗的功率为 输入原始数据形成导纳矩Y阵给定电压初值，置R对于PQ结点，按式计算P(r),Q(r)对于PV结点，按式计算，启动是否|P(r),Q(r)| ？计算雅克比矩阵各元素解修正方程式，求e(r),f(r)用e(r+1)=e(r) e(r), f (r+1)=f(r) f (r)修正结点电压以e(r+1)e(r), f (r+1)f (r)以(r+1)r按系统的潮流分布计算平衡节点功率及线路功率输出3.3 系统程序框图系统程序框图如图3-1所示。置r=0图3-1 系统程序框图3.4 MATLAB简介MATLAB是用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB是一种交互式、面向对象的程序设计语言广泛应用于工业界与学术界主要用于矩阵运算同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。MATLAB程序设计语言结构完整且具有优良的移植性它的基本数据元素是不需要定义的数组。它可以高效率地解决工业计算问题特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言可以用类似数学公式的方式来编写算法大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。目前电子计算机已广泛应用于电力系统的分析计算潮流计算是其基本应用软件之一。现有很多潮流计算方法。对潮流计算方法有五方面的要求（1）计算速度快（2）内存需要少（3）计算结果有良好的可靠性和可信（4）适应性好亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强。4 潮流分布计算4.1 系统的一次接线图系统的一次接线图如图4-1所示。图4-1 系统的一次接线图4.2 参数计算设定基准值，Ub=Uav.n，则各参数如下。（1） 发电机的次暂态电抗如表4-2所示，计算公式如下式。 X=X*Sb/Sn，ZB=UB2/SN表4-2 各电厂发电机次暂态电抗表 电厂装机容量枯水出力比例丰水出力比例丰大有功丰大无功枯大有功枯大无功水电站1300.20.721.00013.0156.0003.718水电站2720.20.750.40031.23514.4008.924水电站3240.20.716.80010.4124.8002.975水电站5300.20.721.00013.0156.0003.718华鑫电厂500.80.840.00024.79040.00024.790 （2）110KV升压变压器的参数如表4-3所示。电阻：R=PK*UN2/ 1000SN2;电抗：X=UK(%)*UN2/ 100 SN;电导：G=P0/ 1000 UN2;电纳：B=I0(%)*SN/ 100 UN2;式中UN以KV为单位，SN以MVA为单位，P0、PK以KW为单位。表4-3 110KV变压器参数 变电站名容量X1X0SbR*X*G*B*A16B28C120D63E20G31.5F12.5+31.5水电站130MW、2*20MVA0.82250.5251000.037364 0.581212 0.000199 0.001445 水电站272MW、3*31.5MVA0.344440.333331000.021199 0.369023 0.000280 0.001992 水电站324MW、31.5MVA1.033330.333331000.021199 0.369023 0.000280 0.001992 水电站54*7.5MW、2*20MVA1.033330.333331000.037364 0.581212 0.000199 0.001445 火电厂50MW、31.5+40MVA0.482850.26251000.021199 0.369023 0.000280 0.001992 1000.013838 0.290606 0.000379 0.002529 （3）220KV三绕组变压器的参数如表4-4所示。电阻：R=PK*UN2/ 1000SN2;电抗：X=UK(%)*UN2/ 100 SN;电导：G=P0/ 1000 UN2;电纳：B=I0(%)*SN/ 100 UN2表4-4 220kV变电站参数表SBR* X* G*B*高压侧绕组中压侧绕组低压侧绕组高压侧绕组（1）1000.002920 0.116650 0.001650 0.012000 容量120120120中压侧绕组（2）100-0.002138 0.071667 0.000499 0.003630 电压22012110.5低压侧绕组（3）1000.003580 -0.004575 0.000004 0.000027 （4）110KV线路参数如表4-5所示。 r=/S；lg ；Xl*=x1*l*SB/U2av.n。表4-5 110KV线路参数表序号线路名称导线牌号线路长度kmUB(KV)SB（MW)R*X*B*1水电站1水电站2150301151000.047637 0.090737 0.011506 2AB95801151000.199622 0.252250 0.029730 3BC9511151000.002495 0.003153 0.000372 4水电站4C150651151000.103214 0.196597 0.024929 5BD95631151000.157202 0.198647 0.023412 6CD240601151000.059887 0.175123 0.023964 7水电站2C240751151000.074858 0.218904 0.029955 8ED24014.71151000.014672 0.042905 0.005871 9水电站3D15010.51151000.016673 0.031758 0.004027 10水电站5T节点705.51151000.018715 0.017758 0.001986 11T节点D7051151000.017013 0.016144 0.001805 12T节点G24051151000.004991 0.014594 0.001997 13H F1503.91151000.006193 0.011796 0.001496 14H 008998 0.021066 0.002740 15水电站2L2401301151000.129754 0.379433 0.051921 16LG240101151000.009981 0.029187 0.003994 17D006352 0.012098 0.001534 18火电厂D24011151000.000998 0.002919 0.000399 110KV变电站负荷参数及其补偿电容参数如表4-6、表4-7所示。表4-6 110KV变电站负荷参数 变电站总容量丰变站比枯站变比丰大有功丰大无功枯大有功枯大无功A0002.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.6773G15.75005.176818.90006.2121表4-7 各变电站补偿电容参数站别运行方式ABCDEFGH低压侧L低压侧丰大34.609.805.901.92.0枯大34.40005.901.92.04.3 丰大及枯大下地潮流分布情况电压是衡量电力系统电能质量的标准之一。电压过高或过低，都将对人身及其用电设备产生重大的影响。保证用户的电压接近额定值是电力系统调度的基本任务之一。当系统的电压偏离允许值时，电力系统必须应用电压调节技术调节系统电压的大小，使其维持在允许值范围内。本文经过手算形成了等值电路图，并编写好了程序得出节点电压标幺值，使其满足所要求的调整范围。我们首先对给定的程序输入部分作了简要的分析，程序开始需要我们确定输入节点数、支路数、平衡母线号、支路参数矩阵、节点参数矩阵。（1）为了保证整个系统潮流计算的完整性，我们把凡具有母线及发电机处均选作节点，这样，我们确定电厂一母线上的发电机作为平衡节点，节点号为，其它机组作为PV节点，节点号为，其余节点均为PQ节点，节点号见等值电路图。（2）确定完节点及编号后，各条支路也相应确定了，我们对各支路参数进行了计算。根据所给实际电路图和题中的已知条件，有以下公式计算各输电线路的阻抗和对地支路电容的标幺值和变压器的阻抗标幺值。4.3.1 该地区变压器的有功潮流分布数据(1)该图为丰大潮流模型图：图4.3.1(1) 丰大潮流模型图该运行方式的电压合理，负荷分配也均匀，但是有些线路的负载率偏低。比如水电厂2-L站的负载率仅仅为21.1% ，M-L站的线路的负载率也只是8.8%。(2)该图为枯大潮流模型图：图4.3.1(2) 枯大潮流模型图该运行方式的电压合理，负荷分配也均匀，有些线路的负载率偏低但是有些线路的负载率则偏高。比如水电厂2-L站的负载率仅仅为9.4%，M-L站的线路的负载率也只是 7.6% D-H站的线路却重载到84%。（1） 该地区变压器的有功潮流分布数据如表4-8所示。 表4-8 有功潮流分布数据表主变丰大(mw)枯大(mw)水1主变10.63水1主变210.43水二主变116.54.6水二主变216.54.6水二主变316.54.6水3主变16.84.8水五主变110.53水五主变210.53火电厂主变1(31.5mva)17.617.6火电厂主变2(40.0mva)22.422.4L站高压9.260.6H站1变高压46.875.1H站2变高压46.875.1（2）重、过载负荷元件统计表丰大运行方式下无重、过载元件。枯大重、过载负荷元件统计表如下表4-9。 表4-9 枯大重、过载负荷元件统计表类型负载率实际视在功率额定容量建议线路0.83870.486.2限D站负荷到34MW以下 5 心得体会 电力系统课程设计是本专业的理论知识与实际引用的桥梁，在这次课程设计中，我强烈感觉到自己在很多方面的不足，对别人的依赖性比较强。我想我会在以后的学习中不断去发现自己的不足，并一一改正，希望在以后的学习中不要犯同样的错误。在参数计算中复习了各个参数的计算方法，并且在电网元件的绘制过程中，须将实际的电网地理分布图转画成与之大致相同的电气一次接线图。在此过程中学会了正确的使用电气元件符号对实际设备的替代，并在此在此同时输入个元件的参数，为下一步做好准备。在正确绘制的条件下，学会正确布局，使一次接线图清清晰、美观、无线路交叉，达到了对地理分布图的大致描述。在十几天的课程设计中，我们尽量按照老师的要求做，但在具体的操作过程中，还是出现了很多的问题。课程设计让我感觉电力系统分析是一门很有用的课程。因为我对它的学到的知识比较少。在很多时候我很多东西都不了解，并且走了很多的弯路。而且我感觉自己的知识不够连贯。这次课程设计后，我一定要重新对电力系统分析这门课程做进一步的了解。对在此过程中遗留下的问题做好好的研究。争取早点对电力系统分析这门课程有个全方位的了解。为在以后的毕业课程设计中多些方案。参考文献1何仰赞,温增银.电力系统分析M.华中科技大学出版社,20062何仰赞,温增银.电力系统分析题解M.华中科技大学出版社,20063王锡凡.现代电力系统分析M.科学出版社,20044于永源,杨绮雯.电力系统分析（第二版）M.中国电力出版社,20075薛定宇.控制系统计算机辅助设计M.清华大学出版社,1996附件一附件二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.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.510水电站5T节点705.511T节点D70512T节点G240513H F1503.914H D185715水电站2L240/213016LG2401017DF150418火电厂D2401 3、发电机各发电机的参数如下：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变电站丰大负荷按该站变电容量的50%估算，枯大负荷按该变电站容量的60%估算。两个220kV变电站的低压侧上各挂10MW的负荷，中压侧各挂20MW负荷。功率因素均为0.95。5、并联电容器两个220kV变电站的低压侧上均装设并联补偿。补偿总量按该变电站容量的20%装设，分组原则以每组电容器的容量不超过10MVar且经济性较好为准。附件三 程序%潮流计算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 the 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)=data(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(number,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)=2sub=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; 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; elseY(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,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,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,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:14foriii=1:numberfor jj=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;for 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 to test used to be 1.0end; e(number,1)=new_data1(number,12); f(number,1)=0;% e(64,1)=0.88;%test 118ieee% f(64,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;% 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)=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+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; 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+1PE(ii,jj)=-qq-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);%TEST+1 end; end;end;来源：(/s/blog_4c4e78800100cmlq.html) - Matlab潮流计算程序（改进）_小周_新浪博客 copy=3.14159;disp(=copy=)for 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+1QF(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+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+1end; end;end;%disp(QF);%disp(QF);%disp(QE);%disp(