电气工程及其自动化专业综合实训一{新}.doc

电气工程及其自动化专业综合实训一
一． 原始资料
网络接线如图，各支路阻抗和各节点功率均以标幺值标于图中，其中节点2连接的实际是发额定功率的发电厂。设节点1的电压保持为1.06，试计算图示系统中的潮流分布。（计算精度10-5）

二、等值电路图

电力系统潮流计算
二、报告要求：
1、写出潮流计算的原理
2、写出功率方程（极坐标或直角坐标）
3、写出雅克比矩阵元素的表达式或PQ分解法中的系数矩阵
4、利用MATLABLE语言或其他语言编程调试
5、心得或体会
6、附程序代码，打印
潮流计算步骤：
1、写出潮流计算的原理

潮流是指在发电机母线上功率被注入网络，而在变（配）电站的母线上接入负荷，其间，功率在网络中流动。对于这种流动的功率，电力生产部门称之为潮流。以电力网络潮流、电压计算为主要内容的电力网络稳态行为特性计算的目的在于估计对用户电力供应的质量以及为电力网运行的安全性与经济性评估提供基础数据。配电网潮流计算是配电网络分析的基础，配电网的网络重构、无功功率优化、状态估计和故障处理都需要用到配电网潮流数据。
　　电力系统稳态运行应满足以下要求:
　　1）满足系统经济性运行的要求，每一台发电机的输出必须接近于预先设定值；
　　2）必须确保联络线潮流低于线路热极限和电力系统稳定极限；
　　3）必须保持某些中枢点母线上的电压水平在容许范围内，必要时用无功功率补偿计划来达到；
　　4）区域电网是互联系统的一部分，必须执行合同规定的输送至邻网的联络线功率计划；
　　5）用故障前的潮流控制策略使事故扰动效应最小化。
　　通常情况下，输电线路电压在轻载时会较高，重载时会较低，电压调整是指在负载由轻载到满载变化过程中实时调整线路电压满足运行要求；对于超高压输电线路，线路电压维持在额定电压的±5%之内, 实际运行时，通常电压调整约为10% 。对于低压输电线路，电压调整数值为10%，包含了变压器本身的电压降落。
　　3.1.1 潮流计算的基本物理量
　　潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，就是在三相平衡稳态状态下计算电力系统中每条母线的电压幅值和相角，其中每一设备如传输线和变压器中的有功和无功潮流，以及各设备的损耗都需要计算出来。
　　潮流计算采用电力系统的单线图，对于任意一条母线i，需要以下四个变量描述：电压幅值Ui、相角，电网供给母线的有功Pi、无功Qi。若某一电力系统有N个节点，则共有4N个变量，对于每条母线，这些变量中的两个指定为输入数据，其它的两个是潮流程序所要计算的未知量。为方便起见，在图3.1中传送给母线i的功率可分为发电机发出和负载吸收两部分。也就是
　　Pi = PGi – PLi
　　Qi = QGi – QLi

图3-1 节点的变量
　　每条母线被归分为以下三种母线类型中的某一类：
　　1）平衡节点，一般一个系统只有一个平衡节点。在潮流分布算出以前，网络中的功率损耗是未知的，因此，至少有一个节点的有功功率P和无功功率Q不能给定。另外必须选定一个节点，制定其电压相角为零，作为其它节点电压相位的参考，这个节点叫基准节点。为了计算方便，常将平衡节点和基准节点设在同一个节点上。为方便起见在本书中把它标号为母线1。平衡节点是电压参考节点，该母线的是给定值，作为输入数据，典型取标幺值。潮流程序计算P1和Q1。因为平衡节点的P、Q事先无法确定，为使潮流计算结果符合实际，常把平衡节点选在有较大调节裕量的发电机节点，潮流计算结束时若平衡节点的有功功率、无功功率和实际情况不符，就要调整其他节点的边界条件以使平衡节点的功率满足实际允许范围。
　　2）PQ节点，Pi和Qi是输入数据。这类节点的有功功率Pi和无功功率Qi是给定的，潮流计算程序计算节点电压幅值Ui和相角。负荷节点和无功功率注入的联络节点都属于这类节点。有些情况下，系统中某些发电厂送出的功率在一定时间内为固定时，该发电厂母线也可以作为PQ节点。在一个典型的潮流程序中绝大多数母线作为PQ节点。
　　3）PU节点（电压控制母线），Pi和Ui是输入数据。这类节点的有功功率Pi和节点电压幅值Ui是给定的，潮流程序计算节点的无功功率Qi和电压相角。这类节点必须具有足够的无功可调容量，用以保持给定的节点电压幅值。在电力系统中这类节点的数目较少。例如与发电机、并联补偿电容器或者静止无功系统相连的母线。设备无功功率最大值QGimax和最小值QGimin都是输入数据。另一个例子是与抽头可调节变压器相连的母线，用潮流程序计算抽头的位置。
　　注意，当母线i是无发电机相连接的负载母线时，Pi=-PLi为负值；也就是说，在图3.1中给母线i提供的有功为负值。如果负荷是感性的，Qi=-QLi为负值。
　　综上所述，若系统中有n个节点，n为网络中除去参考节点外的节点数，本书中以大地为参考节点，选第1个节点为平衡节点，剩下的n-1个节点中有r个PU节点，则有n-r-1个PQ节点。因此，除了平衡节点外，有n-1个节点的注入有功功率、n-r-1个PQ节点的注入无功功率和r个PU节点电压幅值为已知量。
　　3.1.2 潮流计算的数学模型
　　在稳态潮流计算中，电力系统各元件（参数）等效成一个有源网络。将发电机和负荷用无阻抗线从网络中抽出，剩下的是由接地和不接地支路组成的无源线性网络，可以用导纳矩阵（Y）或阻抗矩阵（Z）来描述。
　　采用导纳矩阵时，节点电流和节点电压构成以下方程：
　　（3-1）
　　其中：Y为n×n阶导纳矩阵，其阶数n为网络中除去参考节点外的节点数，如果不考虑网络元件的非线性及变压器的相位偏移，Y为对称矩阵，为n×1维节点注入电流列向量，在电力系统计算中，节点注入电流可理解为该节点电源电流与负荷电流之和，并规定流入节点电流为正。因此仅有负荷的节电电流就为负值，某些仅起联络作用的节点，图3-2中节点n=3，其注入电流为零。为n×1维节点电压列向量。网络中有接地支路时，节点电压通常指该节点的对地电压，以大地作为参考节点，并规定其编号为零。

图3-2 运用节点电压法时的电力网路等值电路
　　对于第i个节点，展开为如下形式：
　　（3-2）
　　若采用阻抗矩阵可表示为：
　　展开为：
　　（3-3）
　　在潮流计算时一般以节点电压方程进行。节点导纳矩阵与阻抗矩阵互为逆阵，在短路计算时可直接利用导纳矩阵求逆得到阻抗矩阵以求得短路点的短路电流。
　　由于实际系统中一般不给出节点电流而是节点功率，因此将式（3-2）中的节点注入电流用节点注入功率来表示为：
　　（3-4）
　　如果节点电压用极坐标表示，令

　　n个节点电力系统的潮流方程的一般形式是
　　（3-5）

　　或
　　（3-6）
　　若采用直角坐标系，节点电压可表示为

　　导纳矩阵元素可以表示为

　　将上述表达式带入式（2-8）的右端，展开并分出实部和虚部，便得
　　（3-7）
　　可见，原来电流电压的线性方程组变换为功率和电压的非线性方程组，式（3-6）（3-7）就是潮流计算的基本方程。它是一组共有n个非线性方程组成的复数方程组，如果把实部和虚部分开便得到2n个实数方程，由该方程组可解出2n个运行参数。但是每一个节点都有P、Q、U、 四个运行变量，共有4n个运行参数，所以要事先给定其余2n个参数。这就要根据节点的分类，将每个节点的4个运行参数中的两个作为原始数据，另外两个作为待求量。
　　3.1.3 潮流计算的约束条件
　　为了保证电力系统的正常运行，潮流问题中某些变量应满足一定的约束条件，常用的约束条件有：
　　(1) 所有节点电压必须满足

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

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

为了保证系统运行的稳定性，要求某些输电线两端的电压相位差不超过一定的数值。因此，潮流计算可以归结为求解一组非线性方程组，并使其解答满足一定的约束条件。如果不能满足，则应修改某些变量。甚至修改系统的运行方式，重新进行计算。
2、写出功率方程（极坐标或直角坐标）
3、写出雅克比矩阵元素的表达式或PQ分解法中的系数矩阵

解：（1）、将阻抗转换为导纳矩阵：

（2）、形成导纳矩阵：

(3)、B′，B〞形成以及其逆矩阵。
由题可知除1为平衡节点外，其他节点均为PQ节点，系数B′，B〞阶数相同，又因对该等值网络不存在去除与有功功率和电压或无功功率和电压大小关系较小因数的可能性，这两个矩阵B′，B〞完全相同，它们就由导纳矩阵的虚数部分中除第一行和第二行的各个元素组成。

B′=B〞=

由此可见，网络的节点导纳矩阵为奇异矩阵，但他的虚数部分的子矩阵B′和B〞，则是非记异矩阵，可以求逆矩阵，其逆矩阵：
=
(4)、计算各节点有功功率不平衡量。
取 ,按下式计算各节点有功率不平衡：

相似的可得：
；
(5)、计算各节点电压的相位角 ：
由下列矩阵方程式：

=

从而可得：

（6）、计算各节点无功功率不平衡 。
按下列计算各节点无功功率不平衡：

相似可得：
（7）、计算各节点电压的大小 。
由下列矩阵方程式：

4、利用MATLABLE语言或其他语言编程调试
PQ分解法潮流程序
//文件输入格式：节点总数n(包括联络节点),支路数zls //
//节点数(发电机和负荷)nb,接地电抗数mdk,迭代精度eps //
//考虑负荷静特性标志kk2(0考虑),平衡节点号，优化标志(0不优化) //
//最大迭代次数it1,支路左右节点号izl[],jzl[],支路电阻zr[],电抗zx[] //
//支路容纳zyk[],节点号nob[]及标志nobt[](0-PQ,-1-PV) //
//发电机和负荷有功、无功pg[],qg[],pl[],ql[] //
//电压v0[](pv节点输入实际值，PQ节点任输入一值) //
//电抗节点号idk[],电抗值dkk[] //

#include "math.h"
#include "stdio.h"

#define NS 2000 //最大节点数
#define NS2 NS * 2
#define NS4 1000 //NS4、NS必须大于2*zls。
#define ZS 3000 //最大支路数
#define ZS2 ZS * 2
#define DKS 200 //最大电抗器数
#define N2 ZS * 4
#define N3 ZS * 8 + NS * 4

FILE *fp1, *fp2;
char inname[12], outname[12];
// fp1输入数据文件指针 fp2输出文件指针
// inname[]输入数据文件名 outname[]输出数据文件名
int n, zls, nb, mdk, mpj, bnsopton, it1, dsd, kk2, nzls;
// 节点总数n(包括联络节点) 支路数(回路数)zls 节点数nb(发电机和负荷)
// 接地电抗数mdk 精度eps 平衡节点号mpj
// 节点优化(标志)bnsopton(=0节点不优化，!=0节点优化)
// 最大迭代次数it1 最低电压或最大功率误差节点号dsd
// 负荷静特性标志(=0考虑负荷静特性)
// 支路数(双回线算一条支路)
int izl[ZS], jzl[ZS], idk[DKS], yds[NS], ydz[NS], iy[ZS2];
// izl[],jzl[],idk[]:分别存放左、右节点号和电抗器节点号。
// yds[]存放各行非零非对角元素的个数。
// ydz[i]是第 i 行第一个非零非对角元素的首地址，
// 即在所有非零非对角元素中的次序号
// iy[]存放列足码。
int nnew[NS4], old[NS], nob[NS], nobt[NS];
// nnew[],old[]存放的是新、旧节点号。
// nnew[i]中为i对应的新号
// nob[]存放的是节点号。nobt[]存放的是节点类型, 0: pq节点, -1: pv节点。
double eps, dsm, vmin, dph, dqh, af[3];
// eps迭代收敛精度，dsm最大功率误差
// vmin:系统最低电压值。dph,dqh:系统有、无功损耗。
// af[0]和af[1]分别是负荷有功功率、无功功率静态特性系数。
double v00;
// v00: 系统平均电压 ci,cj分别作为节点i,j的电压相角的临时存储单元。
double zr[ZS], zx[ZS], zyk[ZS], dkk[DKS], gii[NS], bii[NS], yg[ZS2], yb[ZS2];
double pg[NS], qg[NS], pl[NS], ql[NS], v0[NS], v[NS], va[NS];
// 支路电阻zr[] 支路电抗zx[] 输电线路充电容纳zyk[](y0/2)
// 接地电抗dkk[] 对角元实部gii[] 对角元虚部
// 非对角元实部yg[] 非对角元虚部yb[]
// pg[],qg[],pl[],ql[]:发电机，负荷功率实、虚部
// v[]是电压幅值,va[]是电压相角。
double w[NS2], kg[3], b[NS2];
int newsort[NS4];
// newsort[i]存放i对应的老号

void initial();
void pqflow();
void out();
void dataio();
void bnsopt();
void zlsort(int* nnew);
void printo();
void printy();
void y2();
void ya0();
void yzb();
void jdgl(int kq0);
void bbhl(int kq0);
void calc();
int iabs(int a);
void branch_output();
void newval(double* aa);
void printc();
void iswap();
void swap();
void printf2(double* aa, double* bb, int n);
void calc(int* iu, double* u, double* di, int* nfd, double* b);
void printi(int* aa, int n);
void printf1(double* aa, int n);
int find(int k, int a[], int* z);
void yzb(int t, int* iu, double* u, double* di, int* nfd);
int isgn(int a, int b);
void yy1();
void y3();
void newtoold();

int main(void)

{
initial(); //初始化
pqflow(); //pq潮流计算
out(); //输出节点和支路数据
return 1;
}

int isgn(int a, int b)

{
//**** 本函数功能返回值为a的绝对值b的符号 ****//
//参数1提供值，参数2提供符号//

if (b < 0)
if (a > 0)
a = -a;
return a;
}

int find(int k, int a[], int* z)

{
//**** 本函数查找a[]中是否有fabs(k)有则返回0，无则返回1 ****//
//参数1为待查找量，参数2待搜索数组，参数3返回k在a[]中的次序号//

int i;
for (i = 1; i <= n; i ++)
if(iabs(k) == a[i])
{
*z = i;
return 1;
}
return 0;
}

void oldtonew()

{
//**** 本函数将输入数据中的节点号变成从1开始的连续节点号 ****//

int i, j, k, ii1, ii2, zls2, k1, k2, k3, k4, ip;
zls2 = zls + zls;

for (i = 1; i <= zls2; i ++)
newsort[i] = 0;
ii1 = 0;

for (i = 1; i <= zls; i ++)
{
k = izl[i];
if (!find(k, newsort, &ii2))
{
ii1 ++;
newsort[ii1] = iabs(k);
}
k = jzl[i];
if (!find(k, newsort, &ii2))
{
ii1 ++;
newsort[ii1] = iabs(k);
}
}
for (i = 1; i <= ii1-1; i ++)
{
for (j = i+1; j <= ii1; j ++)
{
if (newsort[i] > newsort[j])
{
k = newsort[i];
newsort[i] = newsort[j];
newsort[j] = k;
}
}
}
for (i = 1; i <= zls; i ++)
{
k = izl[i];
if (find(k, newsort, &ii2))
{
izl[i] = isgn(ii2, k);
}
else
printf("error!");
k = jzl[i];
if (find(k, newsort, &ii2))
{
jzl[i] = isgn(ii2, k);
}
else
printf("error!");
printf("izl[%d] = %d, jzl[%d] = %d\n", i, izl[i], i, jzl[i]);
}
for (i = 1; i <= nb; i ++)
{
for (j = 1; j <= n; j ++)
if (nob[i] == newsort[j])
{
nob[i] = j;
break;
}
printf("nob[%d] = %d\n", i, nob[i]);
}
for (j = 1; j <= n; j ++)
{
if (mpj == newsort[j])
{
mpj = j;
break;
}
}
//电抗器节点号转变
for (j = 1; j <= mdk; j ++)
{
for (i = 1; i <= n; i ++)
{
if (idk[j] == newsort[i])
{
idk[j] = i;
break;
}
}
}
}

void initial()

{
//**** 本函数进行初始化工作 ****//

int i, k1;

dataio();//输入原始数据
oldtonew();//转化为新号
if (bnsopton == 0) //节点不优化，新节点号即为老节点号。
for (i = 1; i <= n; i ++)
{
old[i] = i;
nnew[i] = i;
}
else
bnsopt(); //节点优化
mpj = nnew[mpj];//mpj:平衡节点
zlsort(nnew); // sort the r,x and b
for (i = 1; i <= mdk; i ++)
{
k1 = idk[i];
idk[i] = nnew[k1];
}
for (i = 1; i <= n; i ++)
{
v[i] = v00;
va[i] = 0.0;
} // 所有节点的电压幅值初值都为1.000(v00)，电压相角初值都为0 。
// exchange the node before and after sort
for (i = 1; i <= n; i ++)
yds[i] = 0; // the immediate
for (i = 1; i <= nb; i ++)
{
k1 = nnew[nob[i]];
yds[k1] = nobt[i];
}
for (i = 1; i <= n; i ++)
nobt[i] = yds[i];
newval(pg);
newval(qg);
newval(pl);
newval(ql);
newval(v0);
for (i = 1; i <= n; i ++) // nobt[] is type of node
if (nobt[i] == -1)
v[i] = v0[i]; // nob[] is serials numbe
//nobt[] = -1: pv节点，v0[]存放的是最后一个节点数据，
//对于pv节点，即为该点应维持的电压值。
//nobt[] = 0: pq节点，v0[]存放的是最后一个节点数据，
//对于pq节点，即为系统平均电压值。
printo();
//输出af[]、v00和节点排序后的支路、节点和
//接地电抗数据(仅仅查看中间结果)
ya0();//获得yds[]、ydz[]、列足码iy[]。( P407 )
}

void printo()

{
//**** 输出af[]、v00和节点排序后的支路、节点和接地电抗数据 ****//

int i;
fprintf(fp2, "\n ******AF AND V0 ******\n");
fprintf(fp2, "\n %7.3f%7.3f%7.3f\n", af[0], af[1], v00);
printc('-', 78);
fprintf(fp2, "\n\n *******ZLB*******\n");
for (i = 1; i <= zls; i ++)
{
fprintf(fp2, "\n");
fprintf(fp2, "%8d%8d%8d%8d", izl[i], jzl[i], old[abs(izl[i])], old[abs(jzl[i])]);
fprintf(fp2, "%9.4f%9.4f%9.4f", zr[i], zx[i], zyk[i]);
}
printc('-', 78);
fprintf(fp2, "\n\n*******BUS*******\n");
for (i = 1; i <= nb; i ++)
{
fprintf(fp2, "\n");
fprintf(fp2, "%8d%8d%8d", nob[i],old[nob[i]], nobt[i]);
fprintf(fp2, "%9.4f%9.4f%9.4f%9.4f%9.4f", pg[i], qg[i], pl[i], ql[i], v0[i]);
}
printc('-', 78);
fprintf(fp2,"\n\n******DKK******\n");
for (i = 1; i <= mdk; i ++)
{
fprintf(fp2, "\n");
fprintf(fp2, "%8d%8d%7.4f", idk[i], old[idk[i]], dkk[i]);
}
}
void dataio()
{
//**** 系统数据初始化 ****//

int i;
af[0] = 0.6;
af[1] = 2.0;//af[0]和af[1]分别是负荷有功功率、无功功率静态特性系数。
v00 = 1.000;//系统平均电压

printf("\nplease input the name of data file\n");
scanf("%s", inname);
fp1 = fopen(inname, "r");
printf("\nplease output the name of data file\n");
scanf("%s", outname);
fp2 = fopen(outname, "w");
fscanf(fp1, "%d %d %d %d", &n, &zls, &nb, &mdk);
// the number of node ,branches, node
fscanf(fp1, "%lf %d %d %d %d", &eps, &kk2, &mpj,
&bnsopton, &it1);
//precision, swing node,sort the node,iteration numbers
for (i = 1; i <= zls; i ++)
{
fscanf(fp1, "%d %d", &izl[i], &jzl[i]);
fscanf(fp1, "%lf %lf %lf ", &zr[i], &zx[i], &zyk[i]);
}
for (i = 1; i <= nb; i ++)
{
fscanf(fp1, "%d %d", &nob[i], &nobt[i]);
fscanf(fp1, "%lf %lf %lf %lf %lf", &pg[i], &qg[i], &pl[i],
&ql[i], &v0[i]);
}
for (i = 1; i <= mdk; i ++)
{
fscanf(fp1, "%d %lf", &idk[i], &dkk[i]);
}
fclose(fp1);
}

void pqflow()

{

int kq0, iu1[N2], nfd1[NS], iu2[N2], nfd2[NS];
int i, t;
double u1[N2], u2[N2], di1[NS], di2[NS];
yy1();
yzb(0, iu1, u1, di1, nfd1); //form the B matrix of P-0 iteration
y2();
yzb(1, iu2, u2, di2, nfd2); //form the B matrix of Q-V iteration
t = 0;
kq0 = 0;
kg[0] = kg[1] = 1;
do
{
jdgl(kq0); // calculating the power
bbhl(kq0); // find out the maxi
if (kq0 == 0)
printf("P: %d\t%d\t%f\n", t, dsd, dsm);
else
printf("Q: %d\t%d\t%f\n", t, dsd, dsm);

if (fabs(dsm) > eps)
{
kg[kq0]=1;
if (kq0 == 0)
calc(iu1, u1, di1, nfd1, b);
if (kq0 == 1)
calc(iu2, u2, di2, nfd2, b);
for (i = 1; i <= n; i ++)
{
if(kq0 == 0 )
va[i] = va[i] - b[i] / v00;
else
v[i] = v[i] - b[i];
}
}
else
kg[kq0] = 0;
if(kq0 == 0)
kq0 = 1;
else
{
kq0 = 0;
t ++;
}
if(t > it1)
break;
}while((fabs(dsm) > eps) || (kg[kq0] != 0));
fprintf(fp2, "\n%s%d", "times = ", t);
}

void out()

{
//**** 本函数输出节点和支路数据 ****//

zlsort(old); // recover the data if sorted
// newtoold();
node_output(); // node data
branch_output(); //branch data
printc('-', 78);
printc('*', 78);
fprintf(fp2, "\n");
}

void newval(double* aa)

{
//**** 本函数将旧号换成新号 ****//

int i, k1;
for (i = 1; i <= n; i ++)
b[i] = 0.0;
for (i = 1; i <= nb; i ++)
{
k1 = nnew[nob[i]];
b[k1] = aa[i];
}
for (i = 1; i <= n; i ++)
aa[i] = b[i];
}

void yzb(int t, int* iu, double* u, double* di, int* nfd)

{
//**** 本函数求因子表 ****//

//参数1为标志(t=0 求B',t=1求B'')//
//参数2因子表上三角矩阵非零非对角元素的列足码
//参数3因子表上三角矩阵非零非对角元素的数值
//参数4因子表上三角矩阵对角元素
//参数5因子表上三角各行非零元素个数

int i, j, k, i1, i2;
int jj, jj1, jj2, im, x, fd[NS];
double ai, b[NS];
nfd[1] = 1;
for (i = 1; i <= n; i ++)
{
//nobt[] 存放的是节点类型, 0: pq节点, -1: pv节点。
if (((t != 1) || (nobt[i] != -1)) && i != mpj) // <---|
{ // |
for (j = i + 1; j <= n; j ++) // |
b[j] = 0.0; // |
b[i] = bii[i]; // |
if ((kk2 == 0) && (t == 1) && (nobt[i] != -1))// 存在(t == 1)的情况，不多余。
b[i] = b[i] + af[1] * ql[i] / v0[i] / v0[i];//af[1]
i1 = ydz[i];
i2 = ydz[i + 1] - 1;
for (j = i1; j <= i2; j ++)
{
k = iy[j];
b[k] = yb[j];
}
b[mpj] = 0.0;
if (t == 1)
for (j = 1; j <= n; j ++)
if (nobt[j] == -1)
b[j] = 0.0;
i1 = i - 1;
for (im = 1; im <= i1; im ++)
{
jj1 = nfd[im];
jj2 = nfd[im + 1] - 1;
for (jj = jj1; jj <= jj2; jj ++)
{
if(iu[jj] == i)
{
ai = u[jj] / di[im];
for(k = jj; k <= jj2; k ++)
{
j = iu[k];
b[j] = b[j] - ai * u[k];
}
break;
}
}
}
x = nfd[i];
di[i] = 1.0 / b[i];
ai = di[i];
k = 0;
i1 = i + 1;
for (j = i1; j <= n; j ++)
{
if (fabs(b[j]) > 1.0e-15)
{
u[x] = b[j] * ai;
iu[x] = j;
k++;
x++;
}
}
fd[i] = k;
}
else
{
fd[i]