牛顿-拉夫逊法潮流计算方案说明-毕业论文.doc

武汉理工大学《电力系统分析》课程设计说明书
牛顿-拉夫逊法潮流计算课程设计说明书

目录
摘要 1
1.设计意义与要求 2
1.1设计意义 2
1.2设计要求 2
2.牛顿—拉夫逊算法 3
2.1牛顿算法数学原理： 3
2.2 直角坐标系下牛顿法潮流计算的原理 4
3 详细设计过程 9
3.1节点类型 9
3.2待求量 9
3.3导纳矩阵 9
3.4潮流方程 10
3.5修正方程 11
4.程序设计 14
4.1 节点导纳矩阵的形成 14
4.2 计算各节点不平衡量 15
4.3 雅克比矩阵计算 - 17 -
4.4 LU分解法求修正方程 - 19 -
4.5 计算网络中功率分布 - 22 -
5.结果分析 - 22 -
6.小结 - 25 -
参考文献 - 26 -
附录： - 27 -

0
0
摘要
潮流计算是电力网络设计及运行中最基本的计算，对电力网络的各种设计方案及各种运行方式进行潮流计算，可以得到各种电网各节点的电压，并求得网络的潮流及网络中各元件的电力损耗，进而求得电能损耗。
在数学上是多元非线性方程组的求解问题，求解的方法有很多种。牛顿—拉夫逊法是数学上解非线性方程式的有效方法，有较好的收敛性。将牛顿法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性、稀疏性及节点编号顺序优化等技巧，使牛顿法在收敛性、占用内存、计算速度等方面都达到了一定的要求。
本文以一个具体例子分析潮流计算的具体方法，并运用牛顿—拉夫逊算法求解线性方程

关键词： 电力系统 潮流计算 牛顿—拉夫逊算法

1.设计意义与要求
1.1设计意义
潮流计算是电力系统分析中的一种最基本的计算，他的任务是对给定运行条件确定系统运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布及功率损耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。
具体表现在以下方面：
　　(1)在电网规划阶段,通过潮流计算合理规划电源容量及接入点,合理规划网架,选择无功补偿 方案,满足规划水平的大、小方式下潮流交换控制、调峰、调相、调压的要求。
　　(2)在编制年运行方式时,在预计负荷增长及新设备投运基础上,选择典型方式进行潮流计算,发 现电网中薄弱环节,供调度员日常调度控制参考,并对规划、基建部门提出改进网架结构,加快基建进度的建议。
　　(3)正常检修及特殊运行方式下的潮流计算,用于日运行方式的编制,指导发电厂开机方式,有 功、无功调整方案及负荷调整方案,满足线路、变压器热稳定要求及电压质量要求。
　　(4)预想事故、设备退出运行对静态安全的影响分析及作出预想的运行方式调整方案。
总结为在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方 案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和 最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。
1.2设计要求
1）根据给定的运行条件，确定图中电力系统潮流计算时各节点的类型、待求量；
2）求节点导纳矩阵；
3）给出潮流方程或功率方程的表达式；
4）当用牛顿—拉夫逊法计算潮流时，给出修正方程和迭代收敛条件；

2.牛顿—拉夫逊算法
2.1牛顿算法数学原理：
牛顿法 （Newton Method）：解非线性方程f(x)=0的牛顿(Newton) 法，就是将非线性方程线性化的一种方法。它是解代数方程和超越方程的有效方法之一。
设有单变量非线性方程，给出解的近似值，它与真解的误差为，则将满足，即

将上式左边的函数在附近展成泰勒级数，如果差值很小，二次及以上阶次的各项均可略去得：

这是对于变量的修正量的线性方程式，成为修正方程，解此方程可得修正量

用所求得的去修正近似解，便得

修正后的近似解同真解仍然有误差。为了进一步逼近真解，可以反复进行迭代计算，迭代计算通式是

迭代过程的收敛判据为
式中，和为预先给定的小正数。
牛顿-拉夫逊法实质上就是切线法，是一种逐步线性化的方法，此法不仅用于求单变量方程，也适用于多变量非线性代数方程的有效方法。
牛顿法至少是二阶收敛的，即牛顿法在单根附近至少是二阶收敛的，在重根附近是线性收敛的。
牛顿法收敛很快，而且可求复根，缺点是对重根收敛较慢，要求函数的一阶导数存在。
2.2 直角坐标系下牛顿法潮流计算的原理
采用直角坐标时，节点电压可表示为

导纳矩阵元素则表示为

将上述表示式代入的右端，展开并分出实部和虚部，便得

假定系统中的第1，2，3，···，m号节点为PQ节点，第i个节点的给定功率设为和，对该节点可列写方程

（i=1,2，···，m）
假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一个节点可以列写方程

（i=m+1，m+2，···，n-1）
第n号节点为平衡点，其电压是给定的，故不参加迭代。
以上两个方程组总共包含了2（n-1）个方程，待求的变量有也是2（n-1）个。我们还可看到，上面两个方程式已经具备了方程组的形式。
因此，不难写出如下的修正方程式

式中

上述方程中雅克比矩阵的各元素，可以对上式求偏导数获得。
当时

当时

修正方程式还可以写成分块矩阵的形式

式中，和都是二维列向量；是介方阵。

对于PQ节点

对于PV节点

从以上表达式可以看到，雅克比矩阵有以下特点：
（1） 雅克比矩阵各元素都是节点电压的函数，它们的数值将在迭代过程中不断的改变。
（2） 雅克比矩阵的子块中的元素的表达式只用到导纳矩阵中的对应元素。若，则必有。因此，式中分块形式的雅克比矩阵同节点导纳矩阵一样稀疏，修正方程的求解同样可以用稀疏矩阵的求解技巧
（3） 雅克比矩阵的元素或子块都不具有对称性。
用牛顿-拉夫逊法计算潮流的流程：首先要输入网络的原始数据以及各节点的给定值并形成节点导纳矩阵。输入节点电压初值和，置迭代计数k=0。然后开始进入牛顿法的迭代过程。在进行第k+1次迭代时，其计算步骤如下：
（1） 按上一次迭代计算出的节点电压值和，计算各类节点的不平衡量、和。
（2） 按条件校验收敛，即
<
如果收敛，迭代到此结束，转入计算各线路潮流和平衡节点的功率，并打印输出计算结果。不收敛则继续计算。
（3）计算雅克比矩阵的各元素。
（4）解修正方程式，求节点电压的修正量和。
（5）修正各节点的电压

（6）迭代计数加1，返回第一步继续迭代过程。

输入原始数据

形成节点导纳矩阵

给定节点电压初值

k=0

计算

<

是

否
计算雅克比矩阵各元素

解修正方程式，求，

计算平衡节点功率

输出

图1牛顿-拉夫逊法潮流计算程序框图

3 详细设计过程
3.1节点类型
电力系统潮流计算中，节点一般分为如下几种类型：
PQ节点：节点注入的有功功率无功功率是已知的
PV节点：节点注入的有功功率已知，节点电压幅值恒定，一般由无功储备比较充足的电厂和电站充当；
平衡节点：节点的电压为1*exp(0°)，其注入的有功无功功率可以任意调节，一般由具有调频发电厂充当。
更复杂的潮流计算，还有其他节点，或者是这三种节点的组合，在一定条件下可以相互转换。
对于本题目，节点分析如下：
节点1给出有功功率为2，无功功率为1， PQ节点。
节点2给出有功功率为0.5，电压幅值为1.0，PV节点。
节点3电压相位是0，电压幅值为1，平衡节点。
3.2待求量
节点1待求量是V，；
节点2待求量是Q，；
节点3待求量是P，Q。
3.3导纳矩阵
导纳矩阵分为节点导纳矩阵、结点导纳矩阵、支路导纳矩阵、二端口导纳矩阵。
结点导纳矩阵：对于一个给定的电路(网络)，由其关联矩阵A与支路导纳矩阵Y所确定的矩阵。
支路导纳矩阵：表示一个电路中各支路导纳参数的矩阵。其行数和列数均为电路的支路总数。
二端口导纳矩阵：对应y于二端口网络方程，由二端口参数组成
节点导纳矩阵：以导纳的形式描述电力网络节点注入电流和节点电压关系的矩阵。它给出了电力网络连接关系和元件特性的全部信息，是潮流计算的基础方程式。

本例应用结点导纳矩阵
具体计算时，根据如下公式：

由题给出的导纳可求的节点导纳矩阵如下：

进而节点导纳矩阵为：

3.4潮流方程
网络方程是潮流计算的基础，如果给出电压源或电流源，便可解得电流电压分布。然而，潮流计算中，这些值都是无法准确给定的，这样，就需要列出潮流方程。
对n个节点的网络，电力系统的潮流方程一般形式是

（i=1,2，…，n）

其中,,即PQ分别为节点的有功功率无功功率。

3.5修正方程
计算节点1的不平衡量

计算节点2的不平衡量

节点3是平衡节点，其电压是给定的，故不参加迭代。
根据给定的容许误差，按收敛判据进行校验，以上节点1、2的不平衡量都未满足收敛条件，于是继续以下计算。
修正方程式为： (n=3)

以上雅可比矩阵J中的各元素值是通过求偏导数获得的，对PQ节点来说，
是给定的，因而可以写出

对PV节点来说，给定量是，因此可以列出

当时, 雅可比矩阵中非对角元素为

当时,雅可比矩阵中对角元素为:

代入数值后的修正方程为：

求解修正方程得：

3.6收敛条件

一轮迭代结束，根据收敛条件收敛判据，若等式成立，结果收敛，迭代结束，计算平衡节点的功率和线路潮流计算，否则继续计算雅可比矩阵，解修正方程，直到满足收敛判据。

4.程序设计
4.1 节点导纳矩阵的形成
导纳矩阵元素则表示为

//****************计算导纳矩阵*******************
G[1][1]=1.25;
B[1][1]=-5.5;
G[2][2]=1.3;
B[2][2]=-7;
G[3][3]=1.55;
B[3][3]=-6.5;
G[1][2]=G[2][1]=-0.5;
B[1][2]=B[2][1]=3;
G[1][3]=G[3][1]=-0.75;
B[1][3]=B[3][1]=2.5;
G[2][3]=G[3][2]=-0.8;
B[2][3]=B[3][2]=4;
for(i=1;i<4;i++)
{for(j=1;j<4;j++)
{printf("%f+(%f)j",G[i][j],B[i][j]);
printf(" ");
}
printf("\n");//形成节点导纳矩阵
//*******************************************
}
printf("\n");

4.2 计算各节点不平衡量
假定系统中的第1，2，3···，m号节点为PQ节点，第i个节点的给定功率设为和，对该节点可列写方程

（i=1,2，···，m）
假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一个节点可以列写方程

（i=m+1，m+2，···，n-1）
第n号节点为平衡点，其电压是给定的，故不参加迭代，其计算程序如下：
//计算各节点不平衡量
loop1:
printf("迭代次数k1=%d\n",k1);
for (i=1;i<3;i++)
{float a=0,b=0;
for(j=1;j<4;j++)
{a+=G[i][j]*e[j]-B[i][j]*f[j];
b+=G[i][j]*f[j]+B[i][j]*e[j];
}
P[i]=Ps[i]-(e[i]*a+f[i]*b);//计算有功功率的增量
Q[i]=Qs[i]-(f[i]*a-e[i]*b);//计算无功功率的增量

V22=V2s*V2s-e[2]*e[2];
}
printf("有功功率增量P[1]=%f",P[1]);
printf(" ,");
printf("\n");
printf("有功功率增量P[2]=%f",P[2]);
printf(" ,");
printf("\n");
printf("无功功率增量Q[1]=%f",Q[1]);
printf(" ,")
printf("\n");
printf("电压增量V22=%f",V22);

printf("\n")
16

武汉理工大学《电力系统分析》课程设计说明书
4.3 雅克比矩阵计算
上述方程中雅克比矩阵的各元素，可以对计算各点不平衡量得公式中求偏导数获得。当时

当时

以下为程序：
//****形成雅克比矩阵**********************
for(j=1;j<3;j++)
{if(1==j)
{float c=0,d=0;
int m;
for(m=1;m<4;m++)
{c+=G[1][m]*e[m]-B[1][m]*f[m];

d+=G[1][m]*f[m]+B[1][m]*e[m];
}
J[1*N-1][j*N-1]=-c-G[1][j]*e[1]-B[1][j]*f[1];
J[1*N-1][j*N]=-d+B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N-1]=d+B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N]=-c+G[1][j]*e[1]+B[1][j]*f[1];
}
else
{J[1*N-1][j*N-1]=-G[1][j]*e[1]-B[1][j]*f[1];
J[1*N][j*N]=G[1][j]*e[1]+B[1][j]*f[1];
J[1*N-1][j*N]=B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N-1]=B[1][j]*e[1]-G[1][j]*f[1];
}
}
for(j=1;j<3;j++)
{if(2==j)
{float c=0,d=0;
int m;
for(m=1;m<4;m++)
{c+=G[2][m]*e[m]-B[2][m]*f[m];
d+=G[2][m]*f[m]+B[2][m]*e[m];
}
J[2*N-1][j*N-1]=-c-G[2][j]*e[2]-B[2][j]*f[2];
J[2*N-1][j*N]=-d+B[2][j]*e[2]-G[2][j]*f[2];
J[2*N][j*N-1]=-2*e[2];
J[2*N][j*N]=-2*f[2];
}
else
{J[2*N-1][j*N-1]=-G[2][j]*e[2]-B[2][j]*f[2];
J[2*N][j*N]=0;
J[2*N-1][j*N]=B[2][j]*e[2]-G[2][j]*f[2];
J[2*N][j*N-1]=0;
}
}
printf("雅克比矩阵是:\n");
for(i=1;i<5;i++)
{for(j=1;j<5;j++)
{printf("%f",J[i][j]);
printf(" ");
}
printf("\n");
}

4.4 LU分解法求修正方程
LU分解,又称Gauss消去法,可把任意方阵分解成下三角矩阵的基本变换形式(行交换)和上三角矩阵的乘积。其数学表达式为:A=LU。其中L为下三角矩阵的基本变换形式,U为上三角矩阵。
若有矩阵Ax=b
把矩阵LU分解，求AX=b的问题就等价于求出A=LU后：因为Ly=b可求y，再因为Ux=y，可求出x。原始的求法x=A^（-1）*b,某些情况下，如果矩阵A中的数非常小，我认为不是因为大数除以小数误差大么，1/A算出的误差会很大。但LU可以把A分解成两个都比A大的矩阵的乘积，1/L的误差比1/A小的多。
求修正方程的程序如下

//********计算修正方程*************
for(i=1;iM;i++)
{
sigma2=0;
for(s=0;s<=n-1;s++)
sigma2+=L[i][s]*U[s][n];
L[i][n]=(J[i][n]-sigma2)/U[n][n];
}
}
b[1]=P[1];
b[2]=Q[1];
b[3]=P[2];
b[4]=V22;
for(i=1;i=1;i--)
{
sigma2=0;
for(n=i+1;n",i);
printf("%f",f[i]);
printf(" ,");
}
printf("\n")
4.5 计算网络中功率分布
最后要计算出平衡节点的功率和网络中的功率分布。

5.结果分析
给定节点电压初值,经过五次迭代过程后，得到程序的显示结果如下（取）：

整理可得节点电压和不平衡功率的变化情况，分别于表1和表2所示：
迭代计数k
节点电压

1
0.745335-j0.361142
1.000000-j0.101537
1
2
0.417760-j0.350987
0.995147-j0.149334
1
3
1.318158-j0.380195
1.019448-j0.019812
1
4
0.832411-j0.364272
0.998876-j0.087212
1
5
0.508103-j0.354140
0.995799-j0.135330
1
表1 迭代过程中节点电压变化情况

迭代计数k
节点不平衡量

0
-2.000000
-1.000000
0.500000
0.000000
1
-0.148189
-0.976948
-0.072642
0.000000
2
-0.086588
-0.595284
-0.049568
0.009683
3
-0.653416
-4.467996
-0.365523
-0.039275
4
-0.191594
-1.288891
-0.101196
0.002247
5
-0.084434
-0.585319
-0.049676
0.008383
表2 迭代过程中节点不平衡量变化情况

进行了五次迭代，结果仍然没有收敛。
经过查找相关的资料得到：
“多年的实践证明，牛顿法具有很好的二次收敛性，是求解多元非线性方程的经典算法，至今仍是电力系统潮流计算的主流。因此，一般认为算法不是导致不收敛的原因，潮流不收敛产生的主要原因是计算的初始条件给得不合理，导致潮流方程无解。”

——中国自动化网. 改善调度员潮流计算收敛性的措施
6.小结
通过本次电力系统分析课程设计，使我了解了自己在哪些方面有缺陷。
首先，在拿到本次课设的题目时，就看不懂题目！这就给了自己一个不小的打击。于是我认真地看电力系统分析下册有关潮流计算的牛顿-拉夫逊法！了解了何谓PQ节点，PV节点，平衡节点等。大体上知道了运用牛顿-拉夫逊法的各个步骤！
其次，读题的过程中也遇到了一些麻烦:1.不知道图中节点二的“0.5”和“1”分别指代的是什么？2.各个节点对应的是什么节点？通过自己上网查资料，在网上看到了一些人有关的说法以及相似的题目的解答，这给我不小的启发！我认识到了查找资料的重要性！经过了自己的努力，我知道了以上的答案：1.节点二中的“0.5”是有功功率，“1”代表的是电压幅值。2.节点1是PQ节点，节点2是PV节点，节点3是平衡节点。
最后，遇到了最难难的地方：用程序来实现牛顿-拉夫逊算法！起初是想用matlab程序来进行程序的编写并进行仿真，但是查找了一些资料以后，由于自己对于matlab的了解并不深而且自己学习的知识也不够扎实，给自己带来了不小的麻烦最终无法完成牛顿-拉夫逊的算法。后来在网上找到了有关C语言的相关程序！经过了自己的仔细研读，了解了各个函数的作用。经过了自己的改造将程序完整的编辑了出来，并实现的预期的功能!
通过本次课程设计，我知道了自己学习的知识还不够扎实！很多方面只是应付考试，到了让你做东西的时候确实还是相当困难的!尤其是在编程方面的缺陷！现如今是一个软硬件相结合的时代，其中软件更具有竞争性。因此在今后的学习过程中要端正学习态度。做好每一个细节，不断完善自我，提高自身的学习的水平。为将来的学习和工作打下良好的基础！

参考文献
[1] 何仰赞等.电力系统分析上册[M]． 武汉：华中理工大学出版社.
[2] 何仰赞等.电力系统分析下册[M]． 武汉：华中理工大学出版社.
[3] 诸俊伟等.电力系统分析[M]. 北京：中国电力出版社，1995.
[4] 周全仁等.电网计算与程序设计[M].长沙：湖南科学技术出版社,1983.
[5]丁化成.单片机应用技术[A]. 北京:北京航空航天大学出版社,2000.

附录：
#include
#include
#include
#define N 2
#define M 5
main()
{double G[4][4],B[4][4],J[5][5];
float e[4]={0,1,1,1},f[4]={0},P[4],Q[4],Ps[4]={0,-2,0.5},xe[3],xf[3];
float Qs[4]={0,-1},V2s=1,
float V22,max,P3,Q3;
float a1=0,b1=0;
int i,j,n,s,k1=0;
float L[M][M]={0},U[M][M]={0},sigma1,sigma2,b[M],y[M],x[M];
//****************计算导纳矩阵*******************
G[1][1]=1.25;
B[1][1]=-5.5;
G[2][2]=1.3;
B[2][2]=-7;
G[3][3]=1.55;
B[3][3]=-6.5;
G[1][2]=G[2][1]=-0.5;
B[1][2]=B[2][1]=3;
G[1][3]=G[3][1]=-0.75;
B[1][3]=B[3][1]=2.5;
G[2][3]=G[3][2]=-0.8;
B[2][3]=B[3][2]=4;
for(i=1;i<4;i++)
{for(j=1;j<4;j++)
{printf("%f+(%f)j",G[i][j],B[i][j]);
printf(" ");
}
printf("\n");//形成节点导纳矩阵
//********
}
printf("\n");
//********************************************
//计算各节点不平衡量
loop1:
printf("迭代次数k1=%d\n",k1);
for (i=1;i<3;i++)
{float a=0,b=0;
for(j=1;j<4;j++)
{a+=G[i][j]*e[j]-B[i][j]*f[j];
b+=G[i][j]*f[j]+B[i][j]*e[j];
}
P[i]=Ps[i]-(e[i]*a+f[i]*b);//计算有功功率的增量
Q[i]=Qs[i]-(f[i]*a-e[i]*b);//计算无功功率的增量

V22=V2s*V2s-e[2]*e[2];
}
printf("有功功率增量P[1]=%f",P[1]);
printf(" ,");
printf("\n");
printf("有功功率增量P[2]=%f",P[2]);
printf(" ,");
printf("\n");
printf("无功功率增量Q[1]=%f",Q[1]);
printf(" ,");
printf("\n");
printf("电压增量V22=%f",V22);

printf("\n");
//************筛选出最大值***********************

max=fabs(P[1])>fabs(P[2])?fabs(P[1]):fabs(P[2]);
max=max>fabs(Q[1])?max:fabs(Q[1]);
max=max>fabs(V22)?max:fabs(V22);
printf("max=%f\n",max);
//********************************************
while (k1<=4)
{

//****形成雅克比矩阵****************
for(j=1;j<3;j++)
{if(1==j)
{float c=0,d=0;
int m;
for(m=1;m<4;m++)
{c+=G[1][m]*e[m]-B[1][m]*f[m];
d+=G[1][m]*f[m]+B[1][m]*e[m];
}
J[1*N-1][j*N-1]=-c-G[1][j]*e[1]-B[1][j]*f[1];
J[1*N-1][j*N]=-d+B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N-1]=d+B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N]=-c+G[1][j]*e[1]+B[1][j]*f[1];
}

else
{J[1*N-1][j*N-1]=-G[1][j]*e[1]-B[1][j]*f[1];
J[1*N][j*N]=G[1][j]*e[1]+B[1][j]*f[1];
J[1*N-1][j*N]=B[1][j]*e[1]-G[1][j]*f[1];
J[1*N][j*N-1]=B[1][j]*e[1]-G[1][j]*f[1];
}
}
for(j=1;j<3;j++)
{if(2==j)
{float c=0,d=0;
int m;
for(m=1;m<4;m++)
{c+=G[2][m]*e[m]-B[2][m]*f[m];
d+=G[2][m]*f[m]+B[2][m]*e[m];
}
J[2*N-1][j*N-1]=-c-G[2][j]*e[2]-B[2][j]*f[2];
J[2*N-1][j*N]=-d+B[2][j]*e[2]-G[2][j]*f[2];
J[2*N][j*N-1]=-2*e[2];
J[2*N][j*N]=-2*f[2];
}
else
{J[2*N-1][j*N-1]=-G[2][j]*e[2]-B[2][j]*f[2];
J[2*N][j*N]=0;
J[2*N-1][j*N]=B[2][j]*e[2]-G[2][j]*f[2];
J[2*N][j*N-1]=0;
}
}
printf("雅克比矩阵是:\n");
for(i=1;i<5;i++)
{for(j=1;j<5;j++)
{printf("%f",J[i][j]);
printf(" ");
}
printf("\n");
}
//********计算修正方程*************
for(i=1;i{

for(j=n;j
{
sigma1=0;
for(n=1;n<=i-1;n++)
sigma1+=L[i][n]*y[n];
y[i]=b[i]-sigma1;
}
for(i=M-1;i>=1;i--)
{
sigma2=0;
for(n=i+1;n]+=xf[i];
}
for(i=1;i<3;i++)
{printf("e[%d]=",i);
printf("%f",e[i]);
printf(" ,");
}
for(i=1;i<3;i++)
{printf("f[%d]=",i);
printf("%f",f[i]);
printf(" ,");
}
printf("\n");
k1=k1+1;
goto loop1;
}

for(j=1;j<4;j++)
{a1+=G[3][j]*e[j]-B[3][j]*f[j];
b1+=G[3][j]*f[j]+B[3][j]*e[j];
}
P3=e[3]*a1+f[3]*b1;
Q3=f[3]*a1-e[3]*b1;
printf("P3+jQ3=%f+j%f",P3,Q3);
printf("\n");

- 33 -