随着敲下明天要交的报告的最后一个句号，本学期的在校学习任务就算正式结束了。
已经快三周没有好好休息了吧，期末考完之后关于ACM的东西也没学几天，还错过了各种想参加练手的比赛。。。。本以为“电力系统潮流上机计算”这种需要上机编程的东西对我来说应该很简单，结果被各种神奇的公式卡住了。但编程方面毕竟还有些优势，最后班里的同学都是参照着我提供的数据把程序调对的，小小地自豪一下。比别人提前3天提交论文，FFT刚看了两行就得知学校要举办美赛数学建模预选赛简直坑爹，还有一个队友不在学校，于是和另一个队友两人爆肝4天总算是搞出了些东西，希望能被选上吧。这是我第一次参加建模大赛，学习了不少姿势，也算是对正式比赛的一个不错的演练，哪天写篇博文总结下。
今天继续爆脑一天研究数电课设，数电学得还行吧，但我设计的电路就是运行不出期望的结果，最后发现是 Multisim 破解不完全导致所有的“非门”电路都不能正常工作……明明用的和别人是同一个版本啊……
后天就要离开帝都回家了，明天休息下吧，回去又要开始新的学习了。。。。
最后贴一下自己的“电力系统潮流上机计算”代码吧，写得很水，仅作纪念。
/*
程序说明：电力系统潮流计算程序,只能处理如下输入格式的数据
先PQ节点,再PV节点,最后平衡节点。
平衡节点的电压相角请使用弧度制
输入样例见最后
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
#define N_Node 9 //节点数
#define N_Line 9
#define NodeLoadPath "node.txt" //节点数据路径
#define LineLoadPath "line.txt" //支路数据路径
#define YMatrixPath "Ymatrix.txt" //节点导纳矩阵路径
#define JMatrixPath "Jmatrix.txt" //迭代过程中雅克比矩阵路径
#define PI acos(-1.0) //圆周率
const double STD=0.0001; //目标迭代精度
#define Upmax(a,b) ((a)=(a)>(b)?(a):(b)) //最大值更新
/*
typr==0,平衡节点,first,second,分别代表电压,相角
type==1,PQ节点,first,second,分别代表有功,无功
type==2,PV节点，first,second,分别代表有功，电压
*/
struct Node //节点信息
{
int id,type; //序号,节点类型
double first,second;
void Get (FILE *fp) //节点读入成员函数,传入参数为文件指针
{
fscanf(fp,"%d%d%lf%lf",&id,&type,&first,&second);
}
void Out () //节点输出函数
{
printf("%d %d %lf %lf\n",id,type,first,second);
}
}NodeData[N_Node];
/*
R为电阻,X为电感
type==1,线路支路,C为对地电容
type==2,变压器支路,忽略励磁支路,C为变比
*/
struct Line //支路信息
{
int id,type; //序号,节点类型
int start,end; //链接的头尾节点
double R,X,C; //输入支路的3个参量
void Get (FILE *fp) //
{
fscanf(fp,"%d%d%d%d",&id,&type,&start,&end);
fscanf(fp,"%lf%lf%lf",&R,&X,&C);
}
void Out ()
{
printf("%d %d %d %d ",id,type,start,end);
printf("%lf %lf %lf\n",R,X,C);
}
}LineData[N_Line];
class fushu //自定义复数类
{//定义复数 re+j im
public:
double gre,bim;
fushu () //构造函数
{gre=0,bim=0;}
~fushu () {} //析构函数
void Out ()
{
printf("%8lf+J%8lf\n",gre,bim);
}
//以下通过重载运算符的方式定义复数运算规则
fushu operator + (const fushu a) const
{//加法
struct fushu ans;
ans.gre=gre+a.gre;
ans.bim=bim+a.bim;
return ans;
}
fushu operator - (const fushu a) const
{//减法
struct fushu ans;
ans.gre=gre-a.gre;
ans.bim=bim-a.bim;
return ans;
}
fushu operator * (const double x) const
{//复数乘以实数
struct fushu ans;
ans.gre=gre*x;
ans.bim=bim*x;
return ans;
}
fushu operator * (const fushu x) const
{//两个复数相乘
struct fushu ans;
ans.gre=gre*x.gre-bim*x.bim;
ans.bim=gre*x.bim+bim*x.gre;
return ans;
}
fushu operator ^ (const fushu x) const
{//两个复数各自取共轭,然后相乘
double tbim=-bim,xbim=-x.bim;
struct fushu ans;
ans.gre=gre*x.gre-tbim*xbim;
ans.bim=gre*xbim+tbim*x.gre;
return ans;
}
};
fushu yData[N_Node][N_Node]; //节点导纳矩阵
struct Node_Ans
{//节点的解
int id,type;
fushu u; //电压的实部虚部
double r,alpha; //电压的幅值和相角 (角度制)
double Pin,Qin; //节点的注入功率
double B; //节点的对地导纳
void Out ()
{
printf("%d %d %lf %lf %lf %lf %lf %lf %lf\n",id,type,u.gre,u.bim,r,alpha,Pin,Qin,B);
}
}NodeAns[N_Node];
double CAS[N_Line][N_Line]; //连接支路的端点在该支路的对地导纳
double e[N_Node],f[N_Node]; //节点电压 e+jf
double delta[2*(N_Node-1)]; //功率及电压不平衡量
fushu flow [N_Node][N_Node]; //节点间功率流动
fushu fee[N_Line]; //网损
/*
对于PQ节点,先分别计算 H,N,J,L
对于PV节点,先分别计算 H,N,R,S
其中 R与J,L与S 保存在一起
最后将上述矩阵整合到一起 JMatrix
*/
double H[N_Node-1][N_Node-1],N[N_Node-1][N_Node-1];
double J[N_Node-1][N_Node-1],L[N_Node-1][N_Node-1];
double JMatrix[2*(N_Node-1)][2*(N_Node-1)]; //雅克比矩阵
/*
计算电压的幅值和相角
传入参数为：电压的实部虚部 a+jb,用于保存结果的r,alpha
其中alpha 为角度
*/
void Cal_Vol (double a,double b,double &r,double &alpha)
{
r=sqrt(a*a+b*b);
alpha=atan(b/a);
alpha/=PI;
alpha*=180;
}
bool Input () //读入节点和支路的函数,成功返回true,失败返回false
{
FILE *fp; //以下为读入节点信息
fp=fopen(NodeLoadPath,"r");
if (fp==NULL)
{
printf("Node open fail!");
return false;
}
int i;
for (i=0;i<N_Node;i++)
NodeData[i].Get(fp);
fclose(fp); //以下为支路读入信息
fp=fopen(LineLoadPath,"r");
if (fp==NULL)
{
printf("Line open fail!");
return false;
}
for (i=0;i<N_Line;i++)
LineData[i].Get(fp);
fclose(fp);
return true;
}
/*
计算节点导纳矩阵
*/
void Caculate_YMatrix ()
{
memset(yData,0,sizeof(yData));
int i;
for (i=0;i<N_Line;i++)
{
//由阻抗计算导纳
double tmp=LineData[i].R*LineData[i].R+LineData[i].X*LineData[i].X;
double LG=LineData[i].R/tmp;
double LB=(-1.0)*LineData[i].X/tmp;
int start=LineData[i].start-1;
int end=LineData[i].end-1;
if (LineData[i].type==1)
{//数据为线路支路时
double B=LineData[i].C;
yData[start][start].gre+=LG;
yData[start][start].bim+=LB+B;
yData[end][end].gre+=LG;
yData[end][end].bim+=LB+B;
yData[start][end].gre+=(-1)*LG;
yData[start][end].bim+=(-1)*LB;
yData[end][start].gre+=(-1)*LG;
yData[end][start].bim+=(-1)*LB;
CAS[start][end]+=B;
CAS[end][start]+=B;
}
else if (LineData[i].type==2)
{//数据为变压器支路时
double k=LineData[i].C; //变比
yData[start][start].gre+=LG/(k*k);
yData[start][start].bim+=LB/(k*k);
yData[end][end].gre+=LG;
yData[end][end].bim+=LB;
yData[start][end].gre+=(-1)*LG/k;
yData[start][end].bim+=(-1)*LB/k;
yData[end][start].gre+=(-1)*LG/k;
yData[end][start].bim+=(-1)*LB/k;
CAS[start][end]=(k-1)/k*LB;
CAS[end][start]=(1-k)/k/k*LB;
}
}
}
/*
节点导纳矩阵调试函数
输出节点线路信息，向文件及屏幕输出节点导纳矩阵
分别输出实部和虚部
*/
bool Debug_YMatrix ()
{
int i;
for (i=0;i<N_Node;i++)
NodeData[i].Out();
for (i=0;i<N_Line;i++)
LineData[i].Out();
FILE *fp;
fp=fopen("Ymatrix.txt","w");
if (fp==NULL)
{
printf("Ymatrix open fail!\n");
return false;
}
for (i=0;i<N_Node;i++)
{
for (int j=0;j<N_Node;j++)
{
printf("%8lf ",yData[i][j].gre);
fprintf(fp,"%8lf ",yData[i][j].gre);
}
printf("\n");
fprintf(fp,"\n");
}
for (i=0;i<N_Node;i++)
{
for (int j=0;j<N_Node;j++)
{
printf("%8lf ",yData[i][j].bim);
fprintf(fp,"%8lf ",yData[i][j].bim);
}
printf("\n");
fprintf(fp,"\n");
}
fclose (fp);
return true;
}
/*
共有2*(N_Node-1)个不平衡量
对于PQ节点,delta对应2个值 有功功率不平衡del_P,无功功率不平衡del_Q
对于PV节点,delta对应2个值 有功功率不平衡del_P,电压不平衡量del_V
对于平衡节点,电压幅值和相角已知,不参加电压迭代
*/
void Init_UPQ () //初始化节点电压
{
int i;
for (i=0;i<N_Node;i++)
{
switch (NodeData[i].type)
{
case 1:e[i]=1,f[i]=0;break;
case 2:e[i]=NodeData[i].second,f[i]=0;break;
case 0:e[i]=NodeData[i].first*cos(NodeData[i].second);
f[i]=NodeData[i].first*sin(NodeData[i].second);
break;
default:printf("没有这种节点!\n");
}
}
memset(delta,0,sizeof(delta));
}
void UPQ () //计算电压U,有功P,无功Q各自的不平衡量
{
int i,j;
for (i=0;i<N_Node-1;i++)
{
double tmp; //临时变量
if (NodeData[i].type==1)
{
delta[2*i]=NodeData[i].first;
delta[2*i+1]=NodeData[i].second;
for (j=0;j<N_Node;j++)
{
tmp=e[i]*(yData[i][j].gre*e[j]-yData[i][j].bim*f[j]);
tmp+=f[i]*(yData[i][j].gre*f[j]+yData[i][j].bim*e[j]);
delta[2*i]-=tmp; //有功不平衡量
tmp=f[i]*(yData[i][j].gre*e[j]-yData[i][j].bim*f[j]);
tmp-=e[i]*(yData[i][j].gre*f[j]+yData[i][j].bim*e[j]);
delta[2*i+1]-=tmp; //无功不平衡量
}
}
else if (NodeData[i].type==2)
{
delta[2*i]=NodeData[i].first;
delta[2*i+1]=NodeData[i].second*NodeData[i].second; //U^2
for (j=0;j<N_Node;j++)
{
tmp=e[i]*(yData[i][j].gre*e[j]-yData[i][j].bim*f[j]);
tmp+=f[i]*(yData[i][j].gre*f[j]+yData[i][j].bim*e[j]);
delta[2*i]-=tmp; //有功不平衡量
}
delta[2*i+1]-=(e[i]*e[i]+f[i]*f[i]); //电压^2不平衡量
}
}
}
void Cal_JacobianMatrix () //计算雅克比矩阵
{
int i,j;//初始化四个矩阵为0
memset(H,0,sizeof(H));
memset(N,0,sizeof(N));
memset(J,0,sizeof(J));
memset(L,0,sizeof(L));
//计算H/N/J/L四个矩阵
for (i=0;i<N_Node-1;i++) for (j=0;j<N_Node-1;j++)
{
if (i!=j)
{
H[i][j]=yData[i][j].bim*e[i]-yData[i][j].gre*f[i];
N[i][j]=-yData[i][j].gre*e[i]-yData[i][j].bim*f[i];
if (NodeData[i].type==1) //是PQ节点
{
J[i][j]=-N[i][j];//yData[i][j].bim*f[i]+yData[i][j].gre*e[i];
L[i][j]=H[i][j];//-yData[i][j].gre*f[i]+yData[i][j].bim*e[i];
}
else if (NodeData[i].type==2) //是PV节点
{
J[i][j]=0;
L[i][j]=0;
}
else
{
printf("在计算HNJL时出现未知节点\n");
}
}
else //i==j
{
double tmph=0,tmpn=0;
int k;
for (k=0;k<N_Node;k++) //注意循环长度
{
tmph+=yData[i][k].gre*f[k]+yData[i][k].bim*e[k];
tmpn+=yData[i][k].gre*e[k]-yData[i][k].bim*f[k];
}
H[i][i]=-tmph-yData[i][i].gre*f[i]+yData[i][i].bim*e[i];
N[i][i]=-tmpn-yData[i][i].gre*e[i]-yData[i][i].bim*f[i];
if (NodeData[i].type==1) //是PQ节点
{
J[i][j]=-tmpn+yData[i][i].gre*e[i]+yData[i][i].bim*f[i];
L[i][j]=tmph-yData[i][i].gre*f[i]+yData[i][i].bim*e[i];
}
else if (NodeData[i].type==2) //是PV节点
{
J[i][j]=-2*f[i];
L[i][j]=-2*e[i];
}
else
printf("在计算HNJL时出现未知节点\n");
}
}
for (i=0;i<N_Node-1;i++) for (j=0;j<N_Node-1;j++)
{
JMatrix[2*i][2*j]=-H[i][j];
JMatrix[2*i][2*j+1]=-N[i][j];
JMatrix[2*i+1][2*j]=-J[i][j];
JMatrix[2*i+1][2*j+1]=-L[i][j];
}
}
/*
功率不平衡量及HNJL四个矩阵的调试函数
调用后将会向屏幕输出当前的功率不平衡量及HNJL四个矩阵
正式运行时不建议调用
*/
void Debug_unbalance_HNJL ()
{
int i,j;
printf("有功不平衡量:\n");
for (i=0;i<(N_Node-1);i++)
printf("%8lf ",delta[2*i]);
printf("\n无功及电压不平衡量:\n");
for (i=0;i<(N_Node-1);i++)
printf("%8lf ",delta[2*i+1]);
printf("\n\nH矩阵:");
for (i=0;i<N_Node-1;i++)
for (j=0;j<N_Node-1;j++)
printf(j==N_Node-2?"%lf\n":"%lf ",H[i][j]);
printf("\n\nN矩阵:");
for (i=0;i<N_Node-1;i++)
for (j=0;j<N_Node-1;j++)
printf(j==N_Node-2?"%lf\n":"%lf ",N[i][j]);
printf("\n\nJ矩阵:");
for (i=0;i<N_Node-1;i++)
for (j=0;j<N_Node-1;j++)
printf(j==N_Node-2?"%lf\n":"%lf ",J[i][j]);
printf("\n\nL矩阵:");
for (i=0;i<N_Node-1;i++)
for (j=0;j<N_Node-1;j++)
printf(j==N_Node-2?"%lf\n":"%lf ",L[i][j]);
}
/*
向文件输出每次迭代时的雅克比矩阵,功率不平衡量,当前的电压情况
*/
bool Debug_JMatrix_PQ (int id) //传入参数为当前正在进行的迭代次数
{
int i,j;
FILE *fp;
if (id==1) fp=fopen(JMatrixPath,"w");
else fp=fopen(JMatrixPath,"a");
if (fp==NULL)
{
printf("JMatrix Open fail !\n");
return false;
}
fprintf(fp,"\n\n第%d次迭代时各参数情况如下:\n\n",id);
fprintf(fp,"功率不平衡\t有功\t\t无功\n");
for (i=0;i<N_Node-1;i++)
{
if (NodeData[i].type==1)
fprintf(fp,"节点%d\t\t%lf\t%lf\n",i+1,delta[2*i],delta[2*i+1]);
else if (NodeData[i].type==2)
fprintf(fp,"节点%d\t\t%lf\n",i+1,delta[2*i]);
else
fprintf(fp,"平衡节点");
}
fprintf(fp,"\n电压\te\t\tf\n");
for (i=0;i<N_Node;i++)
{
fprintf(fp,"节点%d\t%lf\t%lf\n",i+1,e[i],f[i]);
}
fprintf(fp,"\n电压\t幅值\t\t相角\n");
for (i=0;i<N_Node;i++)
{
double r,alpha; //电压的幅值和相角
Cal_Vol(e[i],f[i],r,alpha);
fprintf(fp,"节点%d\t%lf\t%lf\n",i+1,r,alpha);
}
fprintf(fp,"\n雅克比矩阵\n");
for (i=0;i<2*(N_Node-1);i++)
{
for (j=0;j<2*(N_Node-1);j++)
fprintf(fp,"%lf ",JMatrix[i][j]);
fprintf(fp,"\n");
}
fclose(fp);
return true;
}
bool Function () //处理修正量,判断迭代精度
{//精度满足要求,返回true,否则返回false
int i;
double MaxDelta=0;
for (i=0;i<N_Node-1;i++)
{
e[i]+=delta[2*i+1];
f[i]+=delta[2*i];
Upmax(MaxDelta,fabs(delta[2*i]));
Upmax(MaxDelta,fabs(delta[2*i+1]));
}
if (MaxDelta>STD) //精度不满足要求
return false;
return true;
}
void SolveEquation (unsigned int Dimension,double FactorMatrix[16][16], double ConstVector[])
{//解方程函数
unsigned int i,j,k;
//消元过程
for(i=0;i<Dimension;i++)
{
//规格化过程:(每次消元前,先规格化,以便以下各行消元时,消元系数直接取待消
//列元素的值即可,也便于回代过程,而运算量并不增加)
for( j = i+1; j < Dimension; j++ )
{
FactorMatrix[i][j] = FactorMatrix[i][j] / FactorMatrix[i][i];
}
ConstVector[i] = ConstVector[i]/FactorMatrix[i][i];
for( j = i+1; j < Dimension; j++ ) //消去第i列(从i+1行到最后一行)
{
if( FactorMatrix[j][i] != 0 ) //如果第j行第i列元素本就是0,则不需本列对应的消元过程
{
for( k = i + 1; k < Dimension; k++ ) //当FactorMatrix[i][k]=0,a[j][k]值不变,可省去运算
if( FactorMatrix[i][k] != 0 )
FactorMatrix[j][k] = FactorMatrix[j][k] - FactorMatrix[j][i] * FactorMatrix[i][k];
}
//常数项的消元
ConstVector[j] = ConstVector[j] - FactorMatrix[j][i]* ConstVector[i];
}
}
//回代过程
for( i = Dimension-1; i > 0; i-- ) //Dimension-2:最后一个变量可直接获得,从n-1个变量求起
{
for(j = Dimension-1; j > i-1 ; j-- )
{
if( FactorMatrix[i-1][j] != 0 )
ConstVector[i-1] = ConstVector[i-1] - FactorMatrix[i-1][j] * ConstVector[j]; //a[i][k]=0时可以不算
}
}
return;
}
/*计算节点的解*/
void Caculate_Node_Ans ()
{
int i,j;
for (i=0;i<N_Node;i++)
{//记录节点的电压等基本信息
NodeAns[i].id=NodeData[i].id;
NodeAns[i].type=NodeData[i].type;
NodeAns[i].u.gre=e[i];
NodeAns[i].u.bim=f[i];
Cal_Vol(e[i],f[i],NodeAns[i].r,NodeAns[i].alpha);
//计算节点的功率信息
if (NodeAns[i].type==1) // PQ节点
{
NodeAns[i].Pin=NodeData[i].first;
NodeAns[i].Qin=NodeData[i].second;
}
else if (NodeAns[i].type==2) //PV节点
{
NodeAns[i].Pin=NodeData[i].first;
double tmp=0; //临时变量
for (j=0;j<N_Node;j++)
{
tmp+=f[i]*(yData[i][j].gre*e[j]-yData[i][j].bim*f[j]);
tmp-=e[i]*(yData[i][j].gre*f[j]+yData[i][j].bim*e[j]);
}
NodeAns[i].Qin=tmp;
}
else if (NodeAns[i].type==0) //平衡节点
{
fushu tmp; //复数型的临时变量
for (j=0;j<N_Node;j++)
{//在自定义的复数运算规则中有 ‘^’运算的实现
tmp=tmp+(yData[i][j]^NodeAns[j].u);
}
fushu s=NodeAns[i].u*tmp;
NodeAns[i].Pin=s.gre;
NodeAns[i].Qin=s.bim;
}
else
{
printf("不存在这种节点,节点类型出错");
break;
}
}
}
//计算节点间功率流动及线路网损
void Caculate_Flow ()
{
memset(flow,0,sizeof(flow));
memset(fee,0,sizeof(fee));
for (int i=0;i<N_Line;i++)
{
int start=LineData[i].start-1;
int end=LineData[i].end-1;
double a,b;
fushu tmp=NodeAns[start].u-NodeAns[end].u;
tmp=tmp^yData[start][end];
flow[start][end]=NodeAns[start].u*tmp;
flow[start][end].bim+=NodeAns[start].r*NodeAns[start].r*CAS[start][end];
flow[start][end].gre=-flow[start][end].gre; //取相反数,参考方向相反
flow[start][end].bim=-flow[start][end].bim;
a=NodeAns[start].r*NodeAns[start].r*CAS[start][end];
start=LineData[i].end-1;
end=LineData[i].start-1;
tmp=NodeAns[start].u-NodeAns[end].u;
tmp=tmp^yData[start][end];
flow[start][end]=NodeAns[start].u*tmp;
flow[start][end].bim+=NodeAns[start].r*NodeAns[start].r*CAS[start][end];
flow[start][end].gre=-flow[start][end].gre; //取相反数,参考方向相反
flow[start][end].bim=-flow[start][end].bim;
b=NodeAns[start].r*NodeAns[start].r*CAS[start][end];
fee[i]=flow[start][end]+flow[end][start]; //网损
}
}
//向文件输出计算结果
bool Output_Ans ()
{
FILE *fp;
fp=fopen("电压和相角.txt","w");
if (fp==NULL)
{
printf("电压和相角.txt打开失败!\n");
return false;
}
int i,j;
fprintf(fp,"\t节点\t电压\t相角\n");
for (i=0;i<N_Node;i++)
fprintf(fp,"\t%d\t%.3lf\t%.3lf\n",i+1,NodeAns[i].r,NodeAns[i].alpha);
fclose(fp);
fp=fopen("节点输入的有功和无功.txt","w");
if (fp==NULL)
{
printf("节点输入的有功和无功.txt打开失败!\n");
return false;
}
fprintf(fp,"\t节点\t有功\t无功\n");
for (i=0;i<N_Node;i++)
fprintf(fp,"\t%d\t%.3lf\t%.3lf\n",i+1,NodeAns[i].Pin,NodeAns[i].Qin);
fclose(fp);
fp=fopen("线路潮流.txt","w");
if (fp==NULL)
{
printf("线路潮流.txt打开失败!\n");
return false;
}
fprintf(fp,"有功\t1\t2\t3\t4\t5\t6\t7\t8\t9\n");
for (i=0;i<N_Node;i++)
{
fprintf(fp,"%d\t",i+1);
for (j=0;j<N_Node;j++)
fprintf(fp,"%.3lf\t",flow[i][j].gre);
fprintf(fp,"\n");
}
fprintf(fp,"\n\n无功\t1\t2\t3\t4\t5\t6\t7\t8\t9\n");
for (i=0;i<N_Node;i++)
{
fprintf(fp,"%d\t",i+1);
for (j=0;j<N_Node;j++)
fprintf(fp,"%.3lf\t",flow[i][j].bim);
fprintf(fp,"\n");
}
fclose(fp);
fp=fopen("线路网损.txt","w");
if (fp==NULL)
{
printf("线路网损.txt打开失败!\n");
return false;
}
fprintf(fp,"\t线路号\t有功\t无功\n");
for (i=0;i<N_Line;i++)
fprintf(fp,"\t%d\t%.3lf\t%.3lf\n",i+1,fee[i].gre,fee[i].bim);
fclose(fp);
return true;
}
int main ()
{
if (Input ()==false) return 0; //读入数据
printf("读入数据成功,正在初始化\n");
Caculate_YMatrix (); //计算节点导纳矩阵
// Debug_YMatrix (); //节点导纳矩阵调试函数
Init_UPQ (); //按节点类型赋初值
printf("开始迭代.......\n");
int cnt=1; //标记迭代次数
do
{
UPQ (); //计算功率电压不平衡量
Cal_JacobianMatrix (); //计算雅克比矩阵
// Debug_unbalance_HNJL (); //调试函数正式运行时不建议调用
Debug_JMatrix_PQ (cnt); //向文件输出迭代过程的中间量
SolveEquation (16,JMatrix,delta); //解方程
cnt++;
}while(cnt<=30 && Function ()==false); //判断本次迭代结果
if (cnt>30)
{
printf("迭代不收敛,程序已退出\n");
return 0;
}
printf("迭代收敛,达到目标精度 %lf 共进行了 %d 次迭代\n",STD,cnt-1);
Caculate_Node_Ans (); //计算各节点的电压,注入功率
Caculate_Flow (); //计算线路的功率流动及网损
if (Output_Ans ()) //向文件输出结果
printf("运算完毕,请查看相关结果文件\n");
else
printf("向文件输出结果失败!\n");
return 0;
}
/*
node.txt
1 1 -1.25 -0.5
2 1 0.000000 0.000000
3 1 -0.9 -0.3
4 1 0.000000 0.000000
5 1 -1.000000 -0.35
6 1 0.000000 0.000000
7 2 1.63 1.025
8 2 0.85 1.025
9 0 1.04 0.000000
line.txt
1 1 1 2 0.01 0.085000001 0.088
2 1 1 6 0.032000002 0.160999998 0.152999997
3 1 2 3 0.017000001 0.092 0.079000004
4 2 2 9 0 0.057999998 1
5 1 3 4 0.039000001 0.170000002 0.179000005
6 1 4 5 0.012 0.101000004 0.104999997
7 2 4 8 0 0.059 1
8 1 5 6 0.018999999 0.071999997 0.075000003
9 2 6 7 0 0.063000001 1
*/
扫一扫
952
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
