电力系统潮流计算牛拉法C++程序代码

#include “stdafx.h“
#include “stdio.h“
#include “NEquation.h“
#include “math.h“

#define NODE_TOTAL_NUM 9
#define LINE_TOTAL_NUM 9
#define NODE_DATA_FILENAME	 “D:\\desktop\\chaoliu\\chaoliu\\Node.txt“
#define LINE_DATA_FILENAME	 “D:\\desktop\\chaoliu\\chaoliu\\Line.txt“
#define Y_DATA_FILENAME		 “D:\\desktop\\chaoliu\\chaoliu\\Y.txt“
#define U_DATA_FILENAME		 “D:\\desktop\\chaoliu\\chaoliu\\U.txt“
#define UBL_DATA_FILENAME	 “D:\\desktop\\chaoliu\\chaoliu\\UBL.txt“
#define JCB_DATA_FILENAME		“D:\\desktop\\chaoliu\\chaoliu\\JCB.txt“
#define DELEF_DATA_FILENAME		“D:\\desktop\\chaoliu\\chaoliu\\DELEF.txt“
#define NEWVALUE_DATA_FILENAME “D:\\desktop\\chaoliu\\chaoliu\\NEWVALUE.txt“
#define NETCOM_DATA_FINENAME	“D:\\desktop\\chaoliu\\chaoliu\\NETCOM.txt“
#define eps 0.0001
#define pi 3.1415926535

struct  NodeData
{
	unsigned int Index; //node index
	unsigned int Type;  //node type :PQ:1PV:2balance point:0
	float FirstInput;   //PQ or PV:actice powerbalance point:V 
	float SecondInput;  //PQ:reactive powerPV:Vbalance point:angle
};

struct  LineData
{
	unsigned int Index;  //line index
	unsigned int Iftran; //line type: not transformer: 1 transformer: 2
	unsigned int start;
	unsigned int end;
	float R;	 //PQ or PV:actice powerbalance point:V 
	float X;	 //PQ:reactive powerPV:Vbalance point:angle
	float B_2;	 //transformer: k not transformer:B_2
};

struct  NodeData gNodeData[NODE_TOTAL_NUM];
struct  LineData gLineData[LINE_TOTAL_NUM];
double YG[NODE_TOTAL_NUM][NODE_TOTAL_NUM] = { 0 } YB[NODE_TOTAL_NUM][NODE_TOTAL_NUM] = { 0 };   //导纳矩阵
double Gl[LINE_TOTAL_NUM] = { 0 } Bl[LINE_TOTAL_NUM] = { 0 };                                   //支路导纳
double e[NODE_TOTAL_NUM] = { 0 } f[NODE_TOTAL_NUM] = { 0 };						             //电压
double P[NODE_TOTAL_NUM] = { 0 } Q[NODE_TOTAL_NUM] = { 0 };                                     //功率
double DeltaP[NODE_TOTAL_NUM-1] = { 0 } DeltaQU2[NODE_TOTAL_NUM-1] = { 0 };                     //两个小的不平衡列向量
double DeltaPQ[2 * （NODE_TOTAL_NUM - 1）];                                                        //不平衡总列向量
double Jacobi[2 * NODE_TOTAL_NUM][2 * NODE_TOTAL_NUM] = { 0 };                                   //雅可比矩阵                           
double Deltaef[2 * （NODE_TOTAL_NUM -1）] = { 0 };                                                 //总电压偏差量
double Deltaf[NODE_TOTAL_NUM - 1] = { 0 } Deltae[NODE_TOTAL_NUM - 1] = { 0 };                   //横纵电压偏差量
int Count_num=0;                                                                                 //迭代次数


void test（）
{
	NEquation ob1;

	ob1.SetSize（2）;
	ob1.Data（0 0） = 1;
	ob1.Data（0 1） = 2;
	ob1.Data（1 0） = 2;
	ob1.Data（1 1） = 1;

	ob1.Value（0） = 4;
	ob1.Value（1） = 6;
	ob1.Run（）;

	printf（“x1=%f\n“ ob1.Value（0））;
	printf（“x2=%f\n“ ob1.Value（1））;
}
void GetData（）
{
	FILE* fp;
	int i;
	fp = fopen（NODE_DATA_FILENAME “rb“）;//open node.txt
	for （i = 0; i < NODE_TOTAL_NUM; i++）
	{
		fscanf（fp