期末总结-电力系统潮流上机计算

拼搏现实的明天。 2022-09-20 11:14 229阅读 0赞

随着敲下明天要交的报告的最后一个句号，本学期的在校学习任务就算正式结束了。

已经快三周没有好好休息了吧，期末考完之后关于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
    */