啊哈磊_编程从这里起步
标题:
求讲解
[打印本页]
作者:
陀飞轮
时间:
2021-1-8 21:09
标题:
求讲解
#include<iostream>
#include<fstream>
#include<math.h>
#include<iomanip>
#include<complex>
using namespace std;
#define KM 50
#define PI 3.1415926
int jzqn(void*p,void*q,int nn);
int main()
{
ifstream in("input.txt");
ofstream out("output.txt");
typedef complex<double> cp;
int i,j,k,l,m,n,o;
char c,s[120];
double r,x,bb,kb,a,d,e,fm,f;
struct jiedian
{
int t;
double u,a1,p,q,a,dp,dq,du,da,e,x;
};
struct zhilu
{
int t,i,j;
cp Z,I,K;
};
in.getline(s,100);
in>>n>>o>>f;
double g[n+1][n+1]={},b[n+1][n+1]={},b1[n][n]={},c1[n][n]={};
struct jiedian jd[n+1]={};
struct zhilu zl[o]={};
for(i=1;i<=3;i++)
in.getline(s,100);
for(k=0;k<o;k++)
{
in>>zl[k].t;
in>>i>>j>>r>>x;
zl[k].i=i;
zl[k].j=j;
real(zl[k].Z)=r;
imag(zl[k].Z)=x;
a=r/(r*r+x*x);
d=x/(r*r+x*x);
if(zl[k].t==1)
{
in>>c>>kb;
g[i][i]+=a;
b[i][i]-=d;
g[j][j]+=a/(kb*kb);
b[j][j]-=d/(kb*kb);
g[i][j]-=a/kb;
g[j][i]-=a/kb;
b[i][j]+=d/kb;
b[j][i]+=d/kb;
real(zl[k].K)=kb;
continue;
}
if(zl[k].t==2)
{
in>>bb>>c;
g[i][i]+=a;
g[j][j]+=a;
g[i][j]-=a;
g[j][i]-=a;
b[i][i]+=bb-d;
b[j][j]+=bb-d;
b[i][j]+=d;
b[j][i]+=d;
real(zl[k].K)=1;
continue;
}
if(zl[k].t==3)
{
in>>c>>c;
g[i][i]+=a;
b[i][i]-=d;
}
}
for(i=1;i<=3;i++)
in.getline(s,120);
e=0;
for(i=1;i<=n;i++)
{
in>>c>>k>>jd[i].t;
if(k==1)
{
in>>c>>c>>c>>c>>jd[i].u>>jd[i].a1;
jd[i].a=jd[i].a1*PI/180;
}
if(k==2)
{
in>>a>>c>>d>>c>>jd[i].u>>c;
jd[i].p=a-d;
}
if(k==3)
{
if(e==0)
{ e=1;
l=i; }
in>>jd[i].p>>jd[i].q>>a>>d>>c>>c;
jd[i].p-=a;
jd[i].q-=d;
jd[i].u=1;
}
if(jd[i].t==1) in>>jd[i].x>>jd[i].e;
else in>>c>>c;
}
m=n+1-l;
double b2[m+1][m+1]={},c2[m+1][m+1]={};
for(i=2;i<=n;i++)
for(j=2;j<=n;j++)
b1[i-1][j-1]=-b[i][j];
for(i=l;i<=n;i++)
for(j=l;j<=n;j++)
b2[i-l+1][j-l+1]=-b[i][j];
jzqn(b1,c1,n-1);
jzqn(b2,c2,m);
k=0;
while(1)
{
fm=0;
for(i=2;i<=(l-1);i++)
{
e=0;
for(j=1;j<=n;j++)
{
a=jd[i].a-jd[j].a;
e+=jd[j].u*(g[i][j]*cos(a)+b[i][j]*sin(a));
}
jd[i].dp=jd[i].p-jd[i].u*e;
if(abs(jd[i].dp)>fm) fm=abs(jd[i].dp);
}
for(i=l;i<=n;i++)
{
e=0;
r=0;
for(j=1;j<=n;j++)
{
x=jd[i].a-jd[j].a;
a=cos(x);
d=sin(x);
e+=jd[j].u*(g[i][j]*a+b[i][j]*d);
r+=jd[j].u*(g[i][j]*d-b[i][j]*a);
}
jd[i].dp=jd[i].p-jd[i].u*e;
jd[i].dq=jd[i].q-jd[i].u*r;
if(abs(jd[i].dp)>fm) fm=abs(jd[i].dp);
if(abs(jd[i].dq)>fm) fm=abs(jd[i].dq);
}
if(fm<f) break;
for(i=2;i<=n;i++)
{
r=0;
for(j=2;j<=n;j++)
r+=c1[i-1][j-1]*jd[j].dp/jd[j].u;
jd[i].da=r;
}
for(i=l;i<=n;i++)
{
r=0;
for(j=l;j<=n;j++)
r+=c2[i-l+1][j-l+1]*jd[j].dq/jd[j].u;
jd[i].du=r;
}
for(i=2;i<=n;i++)
jd[i].a+=jd[i].da;
for(i=l;i<=n;i++)
jd[i].u+=jd[i].du;
k++;
if(k>KM) break;
}
if(k>KM)
{
out<<"潮流不收敛";
return 0;
}
for(i=1;i<=1;i++)
{
e=0;
r=0;
for(j=1;j<=n;j++)
{
x=jd[i].a-jd[j].a;
a=cos(x);
d=sin(x);
e+=jd[j].u*(g[i][j]*a+b[i][j]*d);
r+=jd[j].u*(g[i][j]*d-b[i][j]*a);
}
jd[i].p=jd[i].u*e;
jd[i].q=jd[i].u*r;
}
for(i=2;i<=(l-1);i++)
{
e=0;
for(j=1;j<=n;j++)
{
a=jd[i].a-jd[j].a;
e+=jd[j].u*(g[i][j]*sin(a)-b[i][j]*cos(a));
}
jd[i].q=jd[i].u*e;
}
for(i=1;i<=n;i++)
jd[i].a1=jd[i].a*180/PI;
out<<"一、YN矩阵"<<endl;
out<<"实部"<<endl;
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
out<<left<<setw(8)<<g[i][j]<<" ";
out<<endl;
}
out<<"虚部"<<endl;
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
out<<left<<setw(8)<<b[i][j]<<" ";
out<<endl;
}
out<<"二、潮流计算结果"<<endl;
out<<left<<setw(13)<<"节点编号";
out<<left<<setw(19)<<"电压幅值(p.u.)";
out<<left<<setw(18)<<"电压相角(度)";
out<<left<<setw(19)<<"节点有功(p.u.)";
out<<left<<setw(19)<<"节点无功(p.u.)";
out<<endl;
for(i=1;i<=n;i++)
{
out<<" "<<left<<setw(7)<<i;
out<<left<<setw(15)<<jd[i].u;
out<<left<<setw(13)<<jd[i].a1;
out<<left<<setw(15)<<jd[i].p;
out<<left<<setw(15)<<jd[i].q;
out<<endl;
}
int p;
cp ZF;
for(i=1;i<=2;i++)
in.getline(s,100);
in>>p>>real(ZF)>>imag(ZF);
for(i=1;i<=n;i++)
{
if(jd[i].t==1) b[i][i]-=1/jd[i].x;
if(jd[i].t==2)
{
g[i][i]-=jd[i].p/jd[i].u/jd[i].u;
b[i][i]+=jd[i].q/jd[i].u/jd[i].u;
}
}
cp Y1[n+1];
for(i=1;i<=n;i++)
{
real(Y1[i])=g[i][i];
imag(Y1[i])=b[i][i];
}
double gb[2*n+1][2*n+1]={}, rx[2*n+1][2*n+1]={};
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
{
gb[i][j]=g[i][j];
gb[i][n+j]=-b[i][j];
gb[n+i][j]=b[i][j];
gb[n+i][n+j]=g[i][j];
}
jzqn(gb,rx,2*n);
cp Z1[n+1];
for(i=1;i<=n;i++)
{
real(Z1[i])=rx[i][p];
imag(Z1[i])=rx[n+i][p];
}
cp U0[n+1];
for(i=1;i<=n;i++)
{
real(U0[i])=jd[i].u*cos(jd[i].a);
imag(U0[i])=jd[i].u*sin(jd[i].a);
}
cp IF1,U1[n+1];
IF1=U0[p]/(Z1[p]+ZF);
for(i=1;i<=n;i++)
U1[i]=U0[i]-Z1[i]*IF1;
if(abs(U1[p])<0.0000001)
{
real(U1[p])=0;
imag(U1[p])=0;
}
for(k=0;k<o;k++)
if(zl[k].t!=3)
zl[k].I=(U1[zl[k].i]-U1[zl[k].j]/zl[k].K)/zl[k].Z;
for(i=1;i<=n;i++)
if(jd[i].t==2)
{
g[i][i]+=jd[i].p/jd[i].u/jd[i].u;
b[i][i]-=jd[i].q/jd[i].u/jd[i].u;
}
cp Y2[n+1];
for(i=1;i<=n;i++)
{
real(Y2[i])=g[i][i];
imag(Y2[i])=b[i][i];
}
for(i=1;i<=n;i++)
{
gb[i][i]=g[i][i];
gb[i][n+i]=-b[i][i];
gb[n+i][i]=b[i][i];
gb[n+i][n+i]=g[i][i];
}
float lyh;
for(i=0;i<=2*n;i++)
for(j=0;j<=2*n;j++)
{
lyh=gb[i][j];
gb[i][j]=lyh;
}
jzqn(gb,rx,2*n);
cp Z2[n+1];
for(i=1;i<=n;i++)
{
real(Z2[i])=rx[i][p];
imag(Z2[i])=rx[n+i][p];
}
cp IF2,U2[n+1],A(1,0);
IF2=A/(Z2[p]+ZF);
for(i=1;i<=n;i++)
U2[i]=A-Z2[i]*IF2;
out<<"三、节点的自导纳"<<endl;
out<<"节点序号 精确 近似"<<endl;
for(i=1;i<=n;i++)
{
out<<" "<<left<<setw(5)<<i;
out<<left<<setw(20)<<Y1[i];
out<<Y2[i]<<endl;
}
out<<"四、节点阻抗矩阵的第6列"<<endl;
out<<"序号 精确 近似"<<endl;
for(i=1;i<=n;i++)
{
out<<" "<<left<<setw(4)<<i;
out<<left<<setw(23)<<Z1[i];
out<<Z2[i]<<endl;
}
out<<"五、短路电流If"<<endl;
double ifm1,ifm2,mwc,ifj1,ifj2,jwc;
ifm1=abs(IF1);
ifm2=abs(IF2);
mwc=abs(ifm1-ifm2);
ifj1=arg(IF1)*180/PI;
ifj2=arg(IF2)*180/PI;
jwc=abs(ifj1-ifj2);
out<<" 精确 近似 绝对误差"<<endl;
out<<"模值 ";
out<<left<<setw(12)<<ifm1;
out<<left<<setw(12)<<ifm2;
out<<mwc<<endl;
out<<"角度 ";
out<<left<<setw(12)<<ifj1;
out<<left<<setw(12)<<ifj2;
out<<jwc<<endl;
out<<"六、各节点电压的模值"<<endl;
out<<"节点序号 精确 近似 绝对误差"<<endl;
for(i=1;i<=n;i++)
{
out<<" "<<left<<setw(7)<<i;
out<<left<<setw(12)<<abs(U1[i]);
out<<left<<setw(12)<<abs(U2[i]);
out<<abs(abs(U2[i])-abs(U1[i]))<<endl;
}
out<<"七、各支路电流"<<endl;
out<<"首节点序号 尾节点序号 支路电流"<<endl;
for(k=0;k<o;k++)
if(zl[k].t!=3)
{
out<<" "<<left<<setw(11)<<zl[k].i;
out<<left<<setw(7)<<zl[k].j;
out<<zl[k].I<<endl;
}
out<<"八、短路电流容量"<<endl;
out<<"1.46KA"<<endl;
return 0;
}
int jzqn(void*p,void*q,int nn)
{
double(*a)[nn+1]=(double (*)[nn+1])p;
double(*c)[nn+1]=(double (*)[nn+1])q;
double b[nn+1][nn+1]={};
int i,j,k;
double e;
for(i=1;i<=nn;i++)
for(j=1;j<=nn;j++)
b[i][j]=a[i][j];
for(i=0;i<=nn;i++)
for(j=0;j<=nn;j++)
c[i][j]=0;
for(i=1;i<=nn;i++)
c[i][i]=1;
for(j=1;j<=nn;j++)
{
for(i=j;i<=nn;i++)
if(b[i][j]!=0) break;
if(i!=j)
for(k=1;k<=nn;k++)
{
e=b[i][k];
b[i][k]=b[j][k];
b[j][k]=e;
e=c[i][k];
c[i][k]=c[j][k];
c[j][k]=e;
}
e=b[j][j];
if(e!=1)
for(k=1;k<=nn;k++)
{
b[j][k]=b[j][k]/e;
c[j][k]=c[j][k]/e;
}
for(i=1;i<=nn;i++)
{
if(i==j || b[i][j]==0) continue;
e=b[i][j];
for(k=1;k<=nn;k++)
{
b[i][k]=b[i][k]-e*b[j][k];
c[i][k]=c[i][k]-e*c[j][k];
}
}
}
return 0;
}
复制代码
欢迎光临 啊哈磊_编程从这里起步 (https://bbs.codeaha.com/)
Powered by Discuz! X3.2