电力系统设计教材.doc

电力系统设计教材

———————————————————————————————— 作者：
———————————————————————————————— 日期：

2

山东交通学院

电力系统分析课程设计
报告书

院(部)别 信息科学与电气工程学院
班 级 电气
学 号 110
姓 名
指导教师
时 间 2014.06.9-2014.06.13

课 程 设 计 任 务 书

题 目 电力系统分析课程设计
学 院 信息科学与电气工程学院
专 业 电气工程及其自动化
班 级
学生姓名
学 号 1108

6 月 9 日至 6 月 13 日 共 1 周

指导教师(签字)
院 长(签字)

2014 年 6 月 13日

一、设计内容及要求
复杂网络牛顿—拉夫逊法潮流分析与计算的设计
电力系统潮流计算是电力系统中一项最基本的计算，设计内容为复杂网络潮流计算的计算机算法——牛顿-拉夫逊法。
首先，根据给定的电力系统简图，通过手算完成计算机算法的两次迭代过程，从而加深对牛顿-拉夫逊法的理解，有助于计算机编程的应用。
其次，利用计算机编程对电力系统稳态运行的各参数进行解析和计算；编程完成复杂网络的节点导纳矩阵的形成；电力系统支路改变、节点增减的程序变化；编程完成各元件的功率损耗、各段网络的电压损耗、各点电压、功率大小和方向的计算。
二、设计原始资料
给出一个4节点、环网和多引出的电力系统；
参数给定，可以选用直角坐标表示的牛拉公式计算，也可以选用极坐标表示的牛拉公式计算。
试用牛顿-拉夫逊潮流计算方法求解教材P88图4-10网络的潮流分布。求解精度为10e-6

三、设计完成后提交的文件和图表
1．计算说明书部分

设计报告和手算潮流的步骤及结果

2．图纸部分：

电气接线图及等值电路；
潮流计算的计算机算法，即程序；运算结果等以图片的形式附在设计报告中。
四、进程安排

第一天上午：选题，查资料，制定设计方案；
第一天下午——第三天下午：手算完成潮流计算的要求；
第四天上午——第五天上午：编程完成潮流计算，并对照手算结果，分析误差
第五天下午答辩，交设计报告。
五、主要参考资料

《电力系统分析（第三版）》 于永源主编，中国电力出版社，2007年
《电力系统分析》，何仰赞 温增银编著，华中科技大学出版社，2002年版；
《电力系统分析》，韩桢祥主编，浙江大学出版社，2001年版；
《电力系统稳态分析》，陈珩 编，水利电力出版社；

成绩评定表

指导教师成绩
答辩小组成绩
总评成绩

摘 要
潮流计算，电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。传统的潮流计算程序缺乏图形用户界面，结果显示不直观，难于与其他分析功能集成，网络原始数据输入工作量大且易于出错。结合电力系统的特点，对于复杂电力系统，根据定条件，应用牛顿-拉夫逊法进行计算，在手工计算中，由于涉及大量变量、微分方程、矩阵计算，求解很繁琐，计算不同系统时需要重新计算。运用 MATLAB软件进行仿真潮流计算，图形界面直观，运行稳定，计算准确，提高了计算速度，各个类的有效封装又使程序具有很好的模块性.可维护性和可重用性。
关键字： 潮流计算 牛拉法 Matlab

目录
1潮流计算 1
1.1潮流计算概述 1
1.2潮流计算的要求 2
1.3潮流计算的优势 2
1.4潮流计算的用途 3
2 MATLAB简介 4
2.1MATLAB概述 4
3牛顿-拉夫逊法概述 5
3.1牛顿-拉夫逊基本原理 5
3.2直角牛顿-拉夫逊法潮流计算求解过程 6
4设计过程 12
4.1题目C 12
4.2牛顿拉夫逊法程序流程图 17
4.3设计程序 18
5心得体会 25
6参考文献 25

1潮流计算
1.1潮流计算概述
潮流计算是研究电力系统稳态运行情况的一种基本电气计算，常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。通过潮流计算可以判断电网母线电压、支路电流和功率是否越限，如果有越限，就应采取措施，调整运行方式。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置定整计算、电力系统故障计算和稳定计算等提供原始数据。
具体表现在以下方面：
(1)在电网规划阶段,通过潮流计算,合理规划电源容量及接入点,合理规划网架,选择无功补偿方案,满足规划水平的大、小方式下潮流交换控制、调峰、调相、调压的要求。
(2)在编制年运行方式时,在预计负荷增长及新设备投运基础上,选择典型方式进行潮流计算,发现电网中薄弱环节,供调度员日常调度控制参考,并对规划、基建部门提出改进网架结构,加快基建进度的建议。
(3)正常检修及特殊运行方式下的潮流计算,用于日运行方式的编制,指导发电厂开机方式,有功、无功调整方案及负荷调整方案,满足线路、变压器热稳定要求及电压质量要求。
(4)预想事故、设备退出运行对静态安全的影响分析及作出预想的运行方式调整方案。
总结为在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。所以潮流计算是研究电力系统的一种很重要和基础的计算。
1.2潮流计算的要求
电力系统运行必须满足一定技术和经济上的要求。这些要求够成了潮流问题中某些变量的约束条件，常用的约束条件如下：
1.节点电压应满足

从保证电能质量和供电安全的要求来看，电力系统的所有电气设备都必须运行在额定电压附近。PU节点电压幅值必须按上述条件给定。因此，这一约束条件对PQ节点而言。
2.节点的有功功率和无功功率应满足

PQ节点的有功功率和无功功率，以及PU节点的有功功率，在给定是就必须满足上述条件，因此，对平衡节点的P和Q以及PU节点的Q应按上述条件进行检验。
3.节点之间电压的相位差应满足

为了保证系统运行的稳定性，要求某些输电线路两端的电压相位不超过一定的数值。这一约束的主要意义就在于此。

因此，潮流计算可以归结为求解一组非线性方程组，并使其解答满足一定的约束条件。常用的方法是迭代法和牛顿法，在计算过程中，或得出结果之后用约束条件进行检验。如果不能满足要求，则应修改某些变量的给定值，甚至修改系统的运行方式，重新进行计算。
1.3潮流计算的优势
电力系统潮流计算是电力系统分析中的一种最基本的计算，是对复杂电力系统正常和故障条件下稳态运行状态的计算。潮流计算的目标是求取电力系统在给定运行状态的计算。即节点电压和功率分布，用以检查系统各元件是否过负荷。各点电压是否满足要求，功率的分布和分配是否合理以及功率损耗等。对现有电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和暂态稳定分析都是以潮流计算为基础。潮流计算结果可用于电力系统稳态研究，安全估计或最优潮流等，实际电力系统的潮流计算主要采用牛顿-拉夫逊法。借由MATLAB可以轻松实现计算复杂的电力系统潮流分布。
1.4潮流计算的用途
流计算是电力系统非常重要的分析计算，用以研究系统规划和运行中提出的各种问题。对规划中的电力系统，通过潮流计算可以检验所提出的电力系统规划方案能否满足各种运行方式的要求；对运行中的电力系统，通过潮流计算可以预知各种负荷变化和网络结构的改变会不会危及系统的安全，系统中所有母线的电压是否在允许的范围以内，系统中各种元件(线路、变压器等)是否会出现过负荷，以及可能出现过负荷时应事先采取哪些预防措施等。
潮流计算是电力系统分析最基本的计算。除它自身的重要作用之外，在《电力系统分析综合程序》(PSASP)中，潮流计算还是网损计算、静态安全分析、暂态稳定计算、小干扰静态稳定计算、短路计算、静态和动态等值计算的基础

2 MATLAB简介
2.1MATLAB概述
MATLAB的含义是矩阵实验室（Matrix Laboratory），是美国MathWork公司于1982推出的一套高性能的数值计算可视化软件，，包括MATLAB主程序、SIMULINK 动态系统仿真包和各种专业工具箱它集数值分析，矩阵计算，信号处理和图形显示于一体，构成一个方便的，界面友好的用户环境，具有极强大的计算功能和极高的编程效率，特别适合于科学计算、数值分析、系统仿真和信号处理等任务。MATLAB程学设计语言结构完整，且具有优良的移植性，它的基本数据元素是不需要定义的数组。它可以高效率的解决工业计算问题，特别是关于矩阵和矢量的计算。通过MATLAB语言，可以用类似数学公式的方式来编写算法，大大降低了程序需要的难度别难并节省了时间，从而可把主要的经历集中在算法的构思而不是编程上。学习运用MATLAB计算电力系统潮流分布是本次课程设计的重点，可以说，作为工科学生，会运用MATLAB来解决工程问题已成为必须。到目前为止，MATLAB已发展成为国际上最优秀的科技应用软件之一。它的功能十分强大，不仅仅可以实现计算潮流分布，还可以模拟仿真各式各样的数值系统，工程。这里将借助MATLAB来完成用直角牛顿-拉夫逊法计算电力系统潮流分布。

3牛顿-拉夫逊法概述
3.1牛顿-拉夫逊基本原理
潮流计算的目标是求取电力系统在给定运行状态的计算。即节点电压和功率分布，用以检查系统各元件是否过负荷。各点电压是否满足要求，功率的分布和分配是否合理以及功率损耗等。对现有电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和暂态稳定分析都是以潮流计算为基础。潮流计算结果可用如电力系统稳态研究，安全估计或最优潮流等对潮流计算的模型和方法有直接影响。实际电力系统的潮流技术那主要采用牛顿-拉夫逊法。
牛顿--拉夫逊法(简称牛顿法)在数学上是求解非线性代数方程式的有效方法。其要点是把非线性方程式的求解过程变成反复地对相应的线性方程式进行求解的过程。即通常所称的逐次线性化过程。
对于非线性代数方程组：
即 (3-1-1)
在待求量x的某一个初始估计值附近，将上式展开成泰勒级数并略去二阶及以上的高阶项，得到如下的经线性化的方程组：
(3-1-2)
上式称之为牛顿法的修正方程式。由此可以求得第一次迭代的修正量
(3-1-3)
将和相加，得到变量的第一次改进值。接着就从出发，重复上述计算过程。因此从一定的初值出发，应用牛顿法求解的迭代格式为：

(3-1-4)
(3-1-5)

上两式中：是函数对于变量x的一阶偏导数矩阵，即雅可比矩阵J;k为迭代次数。
有上式可见，牛顿法的核心便是反复形式并求解修正方程式。牛顿法当初始估计值和方程的精确解足够接近时，收敛速度非常快，具有平方收敛特性。
牛顿潮流算法突出的优点是收敛速度快，若选择到一个较好的初值，算法将具有平方收敛特性，一般迭代4~5次便可以收敛到一个非常精确的解。而且其迭代次数与所计算网络的规模基本无关。牛顿法也具有良好的收敛可靠性，对于对以节点导纳矩阵为基础的高斯法呈病态的系统，牛顿法也能可靠收敛。牛顿法所需的内存量及每次迭代所需时间均较高斯法多。
牛顿法的可靠收敛取决于有一个良好的启动初值。如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行的节点上。对于正常运行的系统，各节点电压一般均在额定值附近，偏移不会太大，并且各节点间的相位角差也不大，所以对各节点可以采用统一的电压初值(也称为平直电压)，如假定：
或 (3-1-6)
这样一般能得到满意的结果。但若系统因无功紧张或其它原因导致电压质量很差或有重载线路而节点间角差很大时，仍用上述初始电压就有可能出现问题。解决这个问题的办法可以用高斯法迭代1~2次，以此迭代结果作为牛顿法的初值。也可以先用直流法潮流求解一次以求得一个较好的角度初值，然后转入牛顿法迭代。
3.2直角牛顿-拉夫逊法潮流计算求解过程
以下讨论的是用直角坐标形式的牛顿—拉夫逊法潮流的求解过程。当采用直角坐标时，潮流问题的待求量为各节点电压的实部和虚部两个分量由于平衡节点的电压向量是给定的，因此待求两共需要2(n-1)个方程式。事实上，除了平衡节点的功率方程式在迭代过程中没有约束作用以外，其余每个节点都可以列出两个方程式。对PQ节点来说，是给定的，因而可以写出
（3-2-1）
对PV节点来说，给定量是，因此可以列出
（3-2-2）
求解过程大致可以分为以下步骤：
（1）形成节点导纳矩阵Y
（2）将各节点电压设初值U，
（3）将节点初值代入式(2-2-1)或式(22-2)，求出修正方程式的常数项向量
（4）将节点电压初值代入求式，求出雅可比矩阵元素
（5）求解修正方程，求修正向量
（6）求取节点电压的新值
（7）检查是否收敛，如不收敛，则以各节点电压的新值作为初值自第3步重新开始进行狭义次迭代，否则转入下一步
（8）计算支路功率分布，PV节点无功功率和平衡节点柱入功率。
以直角坐标系形式表示
①. 迭代推算式
采用直角坐标时,节点电压相量及复数导纳可表示为:
(3-2-3)
将以上二关系式代入上式中,展开并分开实部和虚部;假定系统中的第1,2,,m号为P—Q节点,第m+1,m+2,,n-1为P—V节点,根据节点性质的不同,得到如下迭代推算式:
⑴对于PQ节点
(3-2-4)

⑵对于PV节点
(3-2-5)

⑶对于平衡节点
平衡节点只设一个,电压为已知,不参见迭代,其电压为:
(3-2-6)
②. 修正方程式迭代式共包括2(n-1)个方程.选定电压初值及变量修正量符号,代入方程并按泰勒级数展开,略去二次方程及以后各项,得到一组线性方程组或线性化了的方程组，常称修正方程组:

(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):

4设计过程
4.1题目C
试用牛顿-拉夫逊潮流计算方法求解教材P88图4-10网络的潮流分布。求解精度为10e-6

手算结果
节点导纳矩阵
[0.93-j4.26 -0.48+j2.4 -0.45j1.89 0 ]
Yb= [-0.48+j2.4 1.07-j4.71 -0.59+jj2.34 0 ]
[-0.45+j1.89 -0.59+j2.34 1.04-j8.32 j3.7 ]
[ 0 0 j3.7 -j3.33 ]
各节点的实际电压标幺值E
E1=1.05
E2=0.97+0.11j
E3=0.99+0.01j
E4=1.09+0.13j
各节点电压大小
U1=1.05
U2=0.97
U3=0.99
U4=1.1
各节点功率S
S1=0.37+0.45j
S2=-0.55-0.13j
S3=-0.3-0.18j
S4=0.5+0.05j

各条支路的功率损耗Ds
Ds（1,2）=0.0092-0.35j
Ds（2,3）=0.0062-0.036j
Ds（1,3）=0.019-0.39j
Ds（3,4）=0+0.062j
D总=0.0344-0.714j

运算结果

4.2牛顿拉夫逊法程序流程图

4.3程序设计：
%本程序 的功能是用牛顿——拉夫逊法进行潮流计算
% 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('请输入平衡母线节点号：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 0]
B2=[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=zeros(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)+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节点电压给定模值
end
for i=1:n %给定各节点注入功率
S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SL
Bi,i)=B(i,i)+B2(i,5); %i节点无功补偿量
end
%===================================================================
P=real(S);Q=imag(S); %分解出各节点注入的有功和无功功率
ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; %迭代次数ICT1、a；不满足收敛要求的节点数IT2
while IT2~=0 % N0=2*n 雅可比矩阵的阶数；N=N01扩展列
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)*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: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节点有功功率差
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=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(ij1)*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节点电压误差
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; % 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 % 不是第三行 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 %