程序有个错误但是调解不成功，请强人帮忙，先谢谢各位
#include "stdafx.h"#include <iostream>
#include <fstream>
#include "NEquation.h"
using namespace System;
using namespace std;
#include "math.h"
////////////////////////////定义结构体数组
struct NodeData
{
int ID;
int Type;
double P, Q, U, A;
}ND[9];
struct BranchData
{
int ID;
int Type;
int HeadNode;
int EndNode;
double R,X,B,k;
}BD[9];
int main(array<System::String ^> ^args)
{
///////////////////////////////////////////////////输入节点、支路数据
ifstream ob2("F:\\node.txt");
for(int i=0;i<9;i++)
{
double f1, f2;
ob2 >> ND[i].ID >> ND[i].Type >> f1 >> f2;
switch(ND[i].Type)
{
case 1: //1为PQ节点
ND[i].P = f1;
ND[i].Q = f2;
ND[i].U = 1.0;
ND[i].A = 0.0;
break;
case 2: //2为PV节点
ND[i].P = f1;
ND[i].Q = 0;
ND[i].U = f2;
ND[i].A = 0.0;
break;
case 3: //3为平衡节点
ND[i].P = 0;
ND[i].Q = 0;
ND[i].U = f1;
ND[i].A = f2;
break;
}
}
ob2.close();
ifstream ob3("F:\\line.txt");
for(int i=0;i<9;i++)
{
float f1;
ob3 >> BD[i].ID >> BD[i].Type >> BD[i].HeadNode >> BD[i].EndNode >> BD[i].R
>> BD[i].X >> f1;
switch(BD[i].Type)
{
case 1: //1为输电线路
BD[i].B = f1;
BD[i].k = 1.0;
break;
case 2: //2为变压器
BD[i].B = 0;
BD[i].k = f1;
}
}
ob3.close();
////////////////////////////////////////形成节点导纳阵
double YG[9][9], YB[9][9];
for(int i=0;i<9;i++)
for(int j=0;j<9;j++)
YG[i][j]=YB[i][j]=0;
for (int i=0;i<9;i++)
{ double G,B;
int ii=BD[i].HeadNode-1;
int jj=BD[i].EndNode-1;
G=BD[i].R/(BD[i].R*BD[i].R+BD[i].X*BD[i].X);
B=-BD[i].X/(BD[i].R*BD[i].R+BD[i].X*BD[i].X);
YG[ii][ii]+=G;
YG[jj][jj]+=G;
YG[ii][jj]-=G;
YG[jj][ii]-=G;
YB[ii][ii]+=B+BD[i].B/2;
YB[jj][jj]+=B+BD[i].B/2;
YB[ii][jj]-=B;
YB[jj][ii]-=B;
}
////////////////////////////////////////输出节点导纳阵
ofstream ob4("F:\\YMatrix.txt");
for(int i=0;i<9;i++)
{ for(int j=0;j<9;j++)
{if(YB[i][j]>=0)
ob4<<YG[i][j]<<"+"<<YB[i][j]<<"j"<<"\t";
else
ob4<<YG[i][j]<<YB[i][j]<<"j"<<"\t";
}
ob4<<"\n";
}
ob4.close();
/////////////////////////////////////////形成Jacobi矩阵
double E[9],F[9],PS[9],QS[9],US[9]; //给各节点付初值
for(int i=0;i<9;i++)
{
E[i]=F[i]=0;
switch(ND[i].Type)
{case 1:E[i]=1.0;
F[i]=0.0;
break;
case 2:E[i]=ND[i].U;
F[i]=0.0;
break;
case 3:E[i]=1.04;
F[i]=0.0;
break;
}
PS[i]=QS[i]=US[i]=0;
}
int k=0; //迭代次数
int flag=0; //迭代结束标志符
double sum;
do{
for(int i=1;i<9;i++) //计算deltaP
{
sum=0;
for(int j=0;j<9;j++)
sum+=E[i]*(YG[i][j]*E[j]-YB[i][j]*F[j])+F[i]*(YG[i][j]*F[j]+YB[i][j]*E[j]);
PS[i]=ND[i].P-sum;
}
for(int i=3;i<9;i++) //计算PQ节点的deltaQ
{
sum=0;
for(int j=0;j<9;j++)
sum+=F[i]*(YG[i][j]*E[j]-YB[i][j]*F[j])-E[i]*(YG[i][j]*F[j]+YB[i][j]*E[j]);
QS[i]=ND[i].Q-sum;
}
for(int i=1;i<3;i++) //计算PV节点的deltaU
US[i]=ND[i].U*ND[i].U-(E[i]*E[i]+F[i]*F[i]);
double a[9],b[9],sum1,sum2; //求解a[i]、b[i]
for(int i=1;i<9;i++)
{
sum1=0;
sum2=0;
for(int j=0;j<9;j++)
{
if(i!=j)
{
sum1+=YG[i][j]*E[j]-YB[i][j]*F[j];
sum2+=YG[i][j]*F[j]+YB[i][j]*E[j];
}
}
a[i]=YG[i][i]*E[i]-YB[i][i]*F[i]+sum1;
b[i]=YG[i][i]*F[i]+YB[i][i]*E[i]+sum2;
}
///////////////////////////////////////////////////////////求解H、N、J、L、R、S
double H[8][8],N[8][8],J[8][8],L[8][8],R[8][8],S[8][8];
for(int i=1;i<9;i++)
for(int j=1;j<9;j++)
{
H[i][j]=0;
N[i][j]=0;
J[i][j]=0;
L[i][j]=0;
R[i][j]=0;
S[i][j]=0;
}
for(int i=1;i<9;i++) //求解H、N
for(int j=1;j<9;j++)
{
if(i=j)
{
H[i][j]=YG[i][j]*F[i]-YB[i][j]*E[i]+b[i];
N[i][j]=YG[i][j]*E[i]+YB[i][j]*F[i]+a[i];
}
else
{
H[i][j]=YG[i][j]*F[i]-YB[i][j]*E[i];
N[i][j]=YG[i][j]*E[i]+YB[i][j]*F[i];
}
}
for(int i=3;i<9;i++) //求解J、L
for(int j=1;j<9;j++)
{
if(i=j)
{
J[i][j]=-YG[i][j]*E[i]-YB[i][j]*F[i]+a[i];
L[i][j]=YG[i][j]*F[i]-YB[i][j]*E[i]-b[i];
}
else
{
J[i][j]=-YG[i][j]*E[i]-YB[i][j]*F[i];
L[i][j]=YG[i][j]*F[i]-YB[i][j]*E[i];
}
}
for(int i=1;i<3;i++) //求解R、S
for(int j=1;j<9;j++)
{
if(i=j)
{
R[i][i]=2*F[i];
S[i][i]=2*E[i];
}
else
{
R[i][j]=0;
S[i][j]=0;
}
}
///////////////////////////////////////组合生成JAC矩阵
double JAC[16][16];
for(int i=0;i<16;i++)
for(int j=0;j<16;j++)
JAC[i][j]=0;
for(int i=1;i<9;i++)
for(int j=1;j<9;j++)
{
JAC[2*(i-1)][2*(j-1)]=H[i][j];
JAC[2*(i-1)][2*j-1]=N[i][j];
}
for(int i=1;i<3;i++)
for(int j=1;j<9;j++)
{
JAC[2*i-1][2*(j-1)]=R[i][j];
JAC[2*i-1][2*j-1]=S[i][j];
}
for(int i=3;i<9;i++)
for(int j=1;j<9;j++)
{
JAC[2*i-1][2*(j-1)]=J[i][j];
JAC[2*i-1][2*j-1]=L[i][j];
}
//////////////////////////////////////////////////////输出deltaP、deltaQ、deltaU及JAC矩阵
ofstream ob5("F:\\deltaPQU.txt");
for(int i=1;i<9;i++)
ob5<<PS[i]<<"\t"<<QS[i]<<"\t"<<US[i]<<"\t"<<endl;
ob5.close();
ofstream ob6("F:\\JAC.txt");
for(int i=0;i<16;i++)
{
for(int j=0;j<16;j++)
ob6<<JAC[i][j]<<"\t";
ob6<<endl;
}
ob6.close();
///////////////////////////////////////////////////调用方程，求解修正量deltaf、deltae
NEquation rs;
rs.SetSize(16);
for(int i=0;i<16;i++)
for(int j=0;j<16;j++)
rs.Data(i,j)=JAC[i][j];
for(int i=1;i<9;i++)
rs.Value(2*(i-1))=PS[i];
rs.Value(1)=US[1];
rs.Value(3)=US[2];
for(int i=3;i<9;i++)
rs.Value(2*i-1)=QS[i];
rs.Run();
///////////////////////////////////////////////////求出deltae、deltaf最大绝对值，并判断是否达到收敛精度
double temp1=fabs(rs.Value(0));
double temp2=fabs(rs.Value(1));
double U[9],A[9];
for(int i=0;i<8;i++)
{
if(temp1<fabs(rs.Value(2*i)))
temp1=fabs(rs.Value(2*i));
if(temp2<fabs(rs.Value(2*i+1)))
temp2=fabs(rs.Value(2*i+1));
}
if(temp1<0.000001&&temp2<0.000001)
flag=1;
if(flag==1)
{
for(int i=1;i<9;i++)
{
F[i]=F[i]+rs.Value(2*(i-1));
E[i]=E[i]+rs.Value(2*i-1);
U[i]=sqrt(F[i]*F[i]+E[i]*E[i]);
A[i]=atan(F[i]/E[i])*180/3.1415926;
}
ofstream ob7("F:\\result1.txt"); //输出节点电压、相角
for(int i=1;i<9;i++)
ob7<<U[i]<<"\t"<<A[i]<<endl;
ob7<<"迭代次数为："<<k;
ob7.close();
}
else
{
for(int i=1;i<9;i++)
{
F[i]=F[i]+rs.Value(2*(i-1));
E[i]=E[i]+rs.Value(2*i-1);
}
k++;
}
}while(flag==0);
///////////////////////////////////////////求解平衡节点功率
/*
double sump=0,sumq=0,P1=0,Q1=0;
for(int i=0;i<9;i++)
{
sump+=YG[0][i]*E[i]-YB[0][i]*F[i];
sumq+=YG[0][i]*F[i]+YB[0][i]*E[i];
}
P1=E[0]*sump;
Q1=-E[0]*sumq;
ofstream ob8("F:\\result2.txt"); //输出1、2、3节点功率
ob8<<P1<<"\t"<<Q1<<endl;
ob8.close();
*/
double sump[3],sumq[3],P[3],Q[3];
for(int i=0;i<3;i++)
{
sump[i]=0;
sumq[i]=0;
P[i]=0;
Q[i]=0;
}
for(int i=0;i<3;i++)
for(int j=0;j<9;j++)
{
sump[i]+=YG[i][j]*E[j]-YB[i][j]*F[j];
sumq[i]+=YG[i][j]*F[j]+YB[i][j]*E[j];
}
P[0]=E[0]*sump[0];
Q[0]=-E[0]*sumq[0];
P[1]=E[1]*sump[1]+F[1]*sumq[1];
Q[1]=F[1]*sump[1]-E[1]*sumq[1];
P[2]=E[2]*sump[2]+F[2]*sumq[2];
Q[2]=F[2]*sump[2]-E[2]*sumq[2];
ofstream ob8("F:\\result2.txt"); //输出1、2、3节点功率
for(int i=0;i<3;i++)
ob8<<P[i]<<"\t"<<Q[i]<<endl;
ob8.close();
/////////////////////////////////////////输出各支路功率
double S[9][9]={0},A=0,B=0,Pij=0,Qij=0;
ofstream ob9("F:\\result3.txt");
for(int i=0;i<9;i++)
{
int ii=BD[i].HeadNode-1;
int jj=BD[i].EndNode-1;
A=(F[ii]-F[jj])*YB[ii][jj]-(E[ii]-E[jj])*YG[ii][jj];
B=(F[ii]-F[jj])*YG[ii][jj]+(E[ii]-E[jj])*YB[ii][jj];
Pij=E[ii]*A-F[ii]*B;
Qij=E[ii]*B-F[ii]*A;
if(BD[i].Type==1)
{
Qij-=(E[ii]*E[ii]+F[ii]*F[ii])*BD[i].B/2;
}
if(Qij<0)
ob9<<ii+1<<"\t"<<jj+1<<"\t"<<"S"<<"["<<ii+1<<"]"<<"["<<jj+1<<"]"<<"="<<Pij<<Qij<<"j"<<endl;
else
ob9<<ii+1<<"\t"<<jj+1<<"\t"<<"S"<<"["<<ii+1<<"]"<<"["<<jj+1<<"]"<<"="<<Pij<<"+"<<Qij<<"j"<<endl;
}
return 0;
}
