能源行业(电力石油煤炭)

开发平台：
C/C++

main.cpp：源码内容
							// main.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#define n 4
#define N (2*n)
/*************************求逆*********************/
float *inv(float *a,int x)
{ 
int *is,*js,i,j,k,l,u,v;
    float d,p;
    is=(int*)malloc(x*sizeof(int));
    js=(int*)malloc(x*sizeof(int));
    for (k=0; k<=x-1; k++)
      { d=0.0;
        for (i=k; i<=x-1; i++)
        for (j=k; j<=x-1; j++)
          { l=i*x+j; p=fabs(a[l]);
            if (p>d) { d=p; is[k]=i; js[k]=j;}
          }
        if (d+1.0==1.0)
          { free(is); free(js); printf("err**not invn");
            return NULL;
          }
        if (is[k]!=k)
          for (j=0; j<=x-1; j++)
            { u=k*x+j; v=is[k]*x+j;
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (js[k]!=k)
          for (i=0; i<=x-1; i++)
            { u=i*x+k; v=i*x+js[k];
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        l=k*x+k;
        a[l]=1.0/a[l];
        for (j=0; j<=x-1; j++)
          if (j!=k)
            { u=k*x+j; a[u]=a[u]*a[l];}
        for (i=0; i<=x-1; i++)
          if (i!=k)
            for (j=0; j<=x-1; j++)
              if (j!=k)
                { u=i*x+j;
                  a[u]=a[u]-a[i*x+k]*a[k*x+j];
                }
        for (i=0; i<=x-1; i++)
          if (i!=k)
            { u=i*x+k; a[u]=-a[u]*a[l];}
      }
    for (k=x-1; k>=0; k--)
      { if (js[k]!=k)
          for (j=0; j<=x-1; j++)
            { u=k*x+j; v=js[k]*x+j;
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
        if (is[k]!=k)
          for (i=0; i<=x-1; i++)
            { u=i*x+k; v=i*x+is[k];
              p=a[u]; a[u]=a[v]; a[v]=p;
            }
      }
    free(is); free(js);
    return a;
}
void display(float M[],int x) /*显示结果*/
{ 
int i,j;
float *result=NULL,*A=NULL;
/*显示原矩阵*/
printf("nn转置前矩阵为:n");
     for(i=0;i<x;i++)
{
   for(j=0;j<x;j++)
   {   
    printf("%10.4ft", M[i*x+j]);
    
   }
      printf("n");
}
     
A=(float*)malloc(x*x*sizeof(float));
for(i=0;i<x*x;i++) A[i]=M[i];
result=inv(A,N);
if(result!=NULL)
{
   printf("nn转置后矩阵为:n");
   for(i=0;i<x;i++)
   {
   for(j=0;j<x;j++)
   {   
	   M[i*x+j]=result[i*x+j];
    printf("%10.4ft", M[i*x+j]);
    
   }
      printf("n");
   }
}
else
{
   printf("The Matrix is singular ! Can't be inverse!n");
}
free(A);
}
#include <conio.h>
/* Test */
/**********************************************************/
float absu(float x)
{
	if(x>=0)
		return x;
	else 
		return (0-x);
}
int main(int argc, char* argv[])
{
    float P[n]={0},p[n],p1[n];
	float Q[n]={0},q[n],q1[n];
	float im[n],re[n];
	float e[n+1]={0};
    float f[n+1]={0};
     int i,j,k;
	  float sum1=0.0;
	float sum2=0.0;
	int s=1;
	 float G[n+1][n+1],B[n+1][n+1];
     float a,b,sump,sumq;
	 float Jac[2*n][2*n];
	 float result[2*n],res[N];
printf("n步骤一:录入导纳！nn平衡节点默认为最后节点，不进行录入!!nn");
printf("输入电导，行优先:n");
for(i=0;i<n+1;i++)
for(j=i;j<n+1;j++)
{
scanf("%f",&G[i][j]);
G[j][i]=G[i][j];
}
for(i=0;i<n+1;i++)
{
for(j=0;j<n+1;j++)
{
printf("%10.4ft",G[i][j]);
//printf("  ");
}
printf("n");
}
printf("n输入电纳:n");
for(i=0;i<n+1;i++)
for(j=i;j<n+1;j++)
{
scanf("%f",&B[i][j]);
B[j][i]=B[i][j];
}
for(i=0;i<n+1;i++)
{
for(j=0;j<n+1;j++)
{
printf("%10.4ft",B[i][j]);
}
printf("n");
}
printf("n节点导纳矩阵为：n");
for(i=0;i<n+1;i++)
{
for(j=0;j<n+1;j++)
{
	printf("%f+j%f",G[i][j],B[i][j]);
	printf("  ");
}
printf("n");
}	
printf("nn");
printf("步骤二:计算deta P,deta Q:nn");
printf("请输入P:n");
for(i=0;i<n;i++)
scanf("%f",&p1[i]);
printf("请输入Q:n");
for(i=0;i<n;i++)
scanf("%f",&q1[i]);
printf("请输入e:n");
for(i=0;i<n+1;i++)
scanf("%f",&e[i]);
printf("请输入f:n");
for(i=0;i<n+1;i++)
scanf("%f",&f[i]);
while(1)
{
	printf("nn第%d次迭代nn",s);
   for(i=0;i<n;i++)
	{sump=0;
	sumq=0;
	a=1;
	b=1;
	for(j=0;j<n+1;j++)
	{   
		a=G[i][j]*e[j]-B[i][j]*f[j];
        b=G[i][j]*f[j]+B[i][j]*e[j];
		sump=sump+e[i]*a+f[i]*b;
		sumq=sumq+f[i]*a-e[i]*b;
	}
		p[i]=sump;
		q[i]=sumq;
		P[i]=p1[i]-sump;
		Q[i]=q1[i]-sumq;
	}
	printf("deP,deQ分别为：n");
	for(i=0;i<n;i++)
	{
		printf("P[%d]=%f",i,P[i]);
	printf("  ");
	}
	printf("n");
   	for(i=0;i<n;i++)
	{
		printf("Q[%d]=%f",i,Q[i]);
	printf("  ");
	}
	printf("nP,Q分别为：n");
	for(i=0;i<n;i++)
	{
		printf("p[%d]=%f",i,p[i]);
	printf("  ");
	}
	printf("n");
   	for(i=0;i<n;i++)
	{
		printf("q[%d]=%f",i,q[i]);
	printf("  ");
	}
                 
		for(k=0;k<n;k++)
		{   
				sum1=e[k]*e[k]+f[k]*f[k];
				re[k]=(p[k]*e[k]-q[k]*f[k])/sum1;
				im[k]=(p[k]*f[k]+q[k]*e[k])/sum1;
						
		}
		printf("n电流为:nn");
		for(k=0;k<n;k++)
			printf("I[%d]=%f+j%f ",k,re[k],im[k]);
	printf("n步骤三：计算Jaccobi矩阵:nn");
     for(i=0;i<n;i++)
	{
	
			for(j=0;j<n;j++)
			{
				if(j!=i)
				{
					float a=0.0,b=0.0;
					a=(G[i][j]*e[i]+B[i][j]*f[i]);
					b=(B[i][j]*e[i]-G[i][j]*f[i]);
					Jac[2*i][2*j]=a;
					Jac[2*i][2*j+1]=0-b;
					Jac[2*i+1][2*j]=0-b;
					Jac[2*i+1][2*j+1]=0-a;
				}
				else
				{
					
                   Jac[2*i][2*j]=re[i]+G[i][i]*e[i]+B[i][i]*f[i];
				   Jac[2*i][2*j+1]=0.0-B[i][i]*e[i]+G[i][i]*f[i]-im[i];
				    Jac[2*i+1][2*j]=0.0-B[i][i]*e[i]+G[i][i]*f[i]+im[i];
				   Jac[2*i+1][2*j+1]=0.0-G[i][i]*e[i]-B[i][i]*f[i]+re[i];
				}
			}
	}
              printf("nJaccobi矩阵是:n");
             for(i=0;i<2*n;i++)
			 {
              for(j=0;j<2*n;j++)
			  {
            	printf("%10.4ft",Jac[i][j]);
	            //printf("      ");
			  }
                printf("n");
			 }	
     printf("对Jaccobi矩阵处理nn");
	  for(i=0;i<n;i++)
	  {
		  res[2*i]=P[i];
		  res[2*i+1]=Q[i];
	  }
     for(i=0;i<2*n;i++)
			 {
              
            	printf("%f",res[i]);
	            printf("  ");
			  }
               
			
     printf("n第四步:计算Jaccobi阵列nn");
     display(Jac[0],N);
    for(i=0;i<N;i++)
	{
		sum2=0.0;
		for(j=0;j<N;j++)
		sum2=sum2+Jac[i][j]*res[j];
		result[i]=sum2;
	}
	 printf("n电压值为:nn");
	 for(i=0;i<n;i++)
		 {e[i]=e[i]+result[2*i];
		 f[i]=f[i]+result[2*i+1];
		 printf("%f+j%f  ",e[i],f[i]);
	 }
     if(absu(result[0])<0.00001||(getchar()=='o'))
		 break;
	  
	 s++;
}
	return 0;
}

						