搜索
查看: 99|回复: 0

求讲解

[复制链接]
 楼主| 发表于 2021-1-8 21:09:50 | 显示全部楼层 |阅读模式
15啊哈币
  1. #include<iostream>
  2. #include<fstream>
  3. #include<math.h>
  4. #include<iomanip>
  5. #include<complex>
  6. using namespace std;
  7. #define KM 50
  8. #define PI 3.1415926
  9. int jzqn(void*p,void*q,int nn);
  10. int main()
  11. {
  12. ifstream in("input.txt");
  13. ofstream out("output.txt");
  14. typedef complex<double> cp;
  15. int i,j,k,l,m,n,o;
  16. char c,s[120];
  17. double r,x,bb,kb,a,d,e,fm,f;
  18. struct jiedian
  19. {
  20.       int t;
  21.       double u,a1,p,q,a,dp,dq,du,da,e,x;
  22. };
  23. struct zhilu
  24. {
  25.       int t,i,j;
  26.       cp Z,I,K;
  27. };
  28. in.getline(s,100);
  29. in>>n>>o>>f;
  30. double g[n+1][n+1]={},b[n+1][n+1]={},b1[n][n]={},c1[n][n]={};
  31. struct jiedian jd[n+1]={};
  32. struct zhilu zl[o]={};
  33. for(i=1;i<=3;i++)
  34.      in.getline(s,100);
  35. for(k=0;k<o;k++)
  36. {
  37.       in>>zl[k].t;
  38.       in>>i>>j>>r>>x;
  39.       zl[k].i=i;
  40.       zl[k].j=j;
  41.       real(zl[k].Z)=r;
  42.       imag(zl[k].Z)=x;
  43.       a=r/(r*r+x*x);
  44.       d=x/(r*r+x*x);
  45.       if(zl[k].t==1)
  46.       {
  47.             in>>c>>kb;
  48.             g[i][i]+=a;
  49.             b[i][i]-=d;
  50.             g[j][j]+=a/(kb*kb);
  51.             b[j][j]-=d/(kb*kb);
  52.             g[i][j]-=a/kb;
  53.             g[j][i]-=a/kb;
  54.             b[i][j]+=d/kb;
  55.             b[j][i]+=d/kb;
  56.             real(zl[k].K)=kb;
  57.             continue;
  58.       }
  59.       if(zl[k].t==2)
  60.       {
  61.             in>>bb>>c;
  62.             g[i][i]+=a;
  63.             g[j][j]+=a;
  64.             g[i][j]-=a;
  65.             g[j][i]-=a;
  66.             b[i][i]+=bb-d;
  67.             b[j][j]+=bb-d;
  68.             b[i][j]+=d;
  69.             b[j][i]+=d;
  70.             real(zl[k].K)=1;
  71.             continue;
  72.       }
  73.       if(zl[k].t==3)
  74.       {
  75.             in>>c>>c;
  76.             g[i][i]+=a;
  77.             b[i][i]-=d;
  78.       }
  79. }
  80. for(i=1;i<=3;i++)
  81.      in.getline(s,120);
  82. e=0;
  83. for(i=1;i<=n;i++)
  84. {
  85.       in>>c>>k>>jd[i].t;
  86.       if(k==1)
  87.       {
  88.           in>>c>>c>>c>>c>>jd[i].u>>jd[i].a1;
  89.           jd[i].a=jd[i].a1*PI/180;
  90.       }
  91.       if(k==2)
  92.       {
  93.            in>>a>>c>>d>>c>>jd[i].u>>c;
  94.            jd[i].p=a-d;
  95.       }
  96.       if(k==3)
  97.       {
  98.             if(e==0)
  99.             {    e=1;
  100.                   l=i;   }
  101.             in>>jd[i].p>>jd[i].q>>a>>d>>c>>c;
  102.             jd[i].p-=a;
  103.             jd[i].q-=d;
  104.             jd[i].u=1;
  105.       }
  106.       if(jd[i].t==1)  in>>jd[i].x>>jd[i].e;
  107.       else in>>c>>c;
  108. }
  109. m=n+1-l;
  110. double b2[m+1][m+1]={},c2[m+1][m+1]={};
  111. for(i=2;i<=n;i++)
  112.       for(j=2;j<=n;j++)
  113.           b1[i-1][j-1]=-b[i][j];
  114. for(i=l;i<=n;i++)
  115.       for(j=l;j<=n;j++)
  116.           b2[i-l+1][j-l+1]=-b[i][j];
  117. jzqn(b1,c1,n-1);
  118. jzqn(b2,c2,m);
  119. k=0;
  120. while(1)
  121. {
  122.       fm=0;
  123.       for(i=2;i<=(l-1);i++)
  124.       {
  125.             e=0;
  126.             for(j=1;j<=n;j++)
  127.             {
  128.                   a=jd[i].a-jd[j].a;
  129.                   e+=jd[j].u*(g[i][j]*cos(a)+b[i][j]*sin(a));
  130.             }
  131.             jd[i].dp=jd[i].p-jd[i].u*e;
  132.             if(abs(jd[i].dp)>fm) fm=abs(jd[i].dp);
  133.       }
  134.       for(i=l;i<=n;i++)
  135.       {
  136.             e=0;
  137.             r=0;
  138.             for(j=1;j<=n;j++)
  139.             {
  140.                   x=jd[i].a-jd[j].a;
  141.                   a=cos(x);
  142.                   d=sin(x);
  143.                   e+=jd[j].u*(g[i][j]*a+b[i][j]*d);
  144.                   r+=jd[j].u*(g[i][j]*d-b[i][j]*a);
  145.             }
  146.             jd[i].dp=jd[i].p-jd[i].u*e;
  147.             jd[i].dq=jd[i].q-jd[i].u*r;
  148.             if(abs(jd[i].dp)>fm) fm=abs(jd[i].dp);
  149.             if(abs(jd[i].dq)>fm) fm=abs(jd[i].dq);
  150.       }
  151.      if(fm<f) break;
  152.      for(i=2;i<=n;i++)
  153.      {
  154.            r=0;
  155.            for(j=2;j<=n;j++)
  156.                 r+=c1[i-1][j-1]*jd[j].dp/jd[j].u;
  157.            jd[i].da=r;
  158.      }
  159.      for(i=l;i<=n;i++)
  160.      {
  161.            r=0;
  162.            for(j=l;j<=n;j++)
  163.                 r+=c2[i-l+1][j-l+1]*jd[j].dq/jd[j].u;
  164.            jd[i].du=r;
  165.      }
  166.      for(i=2;i<=n;i++)
  167.           jd[i].a+=jd[i].da;
  168.      for(i=l;i<=n;i++)
  169.           jd[i].u+=jd[i].du;
  170.      k++;
  171.      if(k>KM) break;
  172. }
  173. if(k>KM)
  174. {
  175.       out<<"潮流不收敛";
  176.       return 0;
  177. }
  178. for(i=1;i<=1;i++)
  179. {
  180.       e=0;
  181.       r=0;
  182.       for(j=1;j<=n;j++)
  183.       {
  184.             x=jd[i].a-jd[j].a;
  185.             a=cos(x);
  186.             d=sin(x);
  187.             e+=jd[j].u*(g[i][j]*a+b[i][j]*d);
  188.             r+=jd[j].u*(g[i][j]*d-b[i][j]*a);
  189.       }
  190.       jd[i].p=jd[i].u*e;
  191.       jd[i].q=jd[i].u*r;
  192. }
  193. for(i=2;i<=(l-1);i++)
  194. {
  195.       e=0;
  196.       for(j=1;j<=n;j++)
  197.       {
  198.             a=jd[i].a-jd[j].a;
  199.             e+=jd[j].u*(g[i][j]*sin(a)-b[i][j]*cos(a));
  200.       }
  201.       jd[i].q=jd[i].u*e;
  202. }
  203. for(i=1;i<=n;i++)
  204.       jd[i].a1=jd[i].a*180/PI;
  205. out<<"一、YN矩阵"<<endl;
  206. out<<"实部"<<endl;
  207. for(i=1;i<=n;i++)
  208.    {
  209.       for(j=1;j<=n;j++)
  210.              out<<left<<setw(8)<<g[i][j]<<"   ";
  211.       out<<endl;
  212.    }
  213. out<<"虚部"<<endl;
  214. for(i=1;i<=n;i++)
  215.    {
  216.       for(j=1;j<=n;j++)
  217.              out<<left<<setw(8)<<b[i][j]<<"   ";
  218.       out<<endl;
  219.    }
  220. out<<"二、潮流计算结果"<<endl;
  221. out<<left<<setw(13)<<"节点编号";
  222. out<<left<<setw(19)<<"电压幅值(p.u.)";
  223. out<<left<<setw(18)<<"电压相角(度)";
  224. out<<left<<setw(19)<<"节点有功(p.u.)";
  225. out<<left<<setw(19)<<"节点无功(p.u.)";
  226. out<<endl;
  227. for(i=1;i<=n;i++)
  228. {
  229.       out<<"   "<<left<<setw(7)<<i;
  230.       out<<left<<setw(15)<<jd[i].u;
  231.       out<<left<<setw(13)<<jd[i].a1;
  232.       out<<left<<setw(15)<<jd[i].p;
  233.       out<<left<<setw(15)<<jd[i].q;
  234.       out<<endl;
  235. }
  236. int p;
  237. cp ZF;
  238. for(i=1;i<=2;i++)
  239.      in.getline(s,100);
  240. in>>p>>real(ZF)>>imag(ZF);
  241. for(i=1;i<=n;i++)
  242. {
  243.       if(jd[i].t==1) b[i][i]-=1/jd[i].x;
  244.       if(jd[i].t==2)
  245.     {
  246.          g[i][i]-=jd[i].p/jd[i].u/jd[i].u;
  247.          b[i][i]+=jd[i].q/jd[i].u/jd[i].u;
  248.      }
  249. }
  250. cp Y1[n+1];
  251. for(i=1;i<=n;i++)
  252. {
  253.       real(Y1[i])=g[i][i];
  254.       imag(Y1[i])=b[i][i];
  255. }
  256. double gb[2*n+1][2*n+1]={}, rx[2*n+1][2*n+1]={};
  257. for(i=1;i<=n;i++)
  258.       for(j=1;j<=n;j++)
  259.       {
  260.            gb[i][j]=g[i][j];
  261.            gb[i][n+j]=-b[i][j];
  262.            gb[n+i][j]=b[i][j];
  263.            gb[n+i][n+j]=g[i][j];
  264.       }
  265. jzqn(gb,rx,2*n);
  266. cp Z1[n+1];
  267. for(i=1;i<=n;i++)
  268.     {
  269.       real(Z1[i])=rx[i][p];
  270.       imag(Z1[i])=rx[n+i][p];
  271.     }
  272. cp U0[n+1];
  273. for(i=1;i<=n;i++)
  274.     {
  275.       real(U0[i])=jd[i].u*cos(jd[i].a);
  276.       imag(U0[i])=jd[i].u*sin(jd[i].a);
  277.     }
  278. cp IF1,U1[n+1];
  279. IF1=U0[p]/(Z1[p]+ZF);
  280. for(i=1;i<=n;i++)
  281.      U1[i]=U0[i]-Z1[i]*IF1;
  282. if(abs(U1[p])<0.0000001)
  283. {
  284.       real(U1[p])=0;
  285.       imag(U1[p])=0;
  286. }
  287. for(k=0;k<o;k++)
  288.      if(zl[k].t!=3)
  289.           zl[k].I=(U1[zl[k].i]-U1[zl[k].j]/zl[k].K)/zl[k].Z;
  290. for(i=1;i<=n;i++)
  291.       if(jd[i].t==2)
  292.       {
  293.             g[i][i]+=jd[i].p/jd[i].u/jd[i].u;
  294.             b[i][i]-=jd[i].q/jd[i].u/jd[i].u;
  295.       }
  296. cp Y2[n+1];
  297. for(i=1;i<=n;i++)
  298. {
  299.       real(Y2[i])=g[i][i];
  300.       imag(Y2[i])=b[i][i];
  301. }
  302. for(i=1;i<=n;i++)
  303. {
  304.       gb[i][i]=g[i][i];
  305.       gb[i][n+i]=-b[i][i];
  306.       gb[n+i][i]=b[i][i];
  307.       gb[n+i][n+i]=g[i][i];
  308. }
  309. float lyh;
  310. for(i=0;i<=2*n;i++)
  311.     for(j=0;j<=2*n;j++)
  312.     {
  313.          lyh=gb[i][j];
  314.          gb[i][j]=lyh;
  315.     }
  316. jzqn(gb,rx,2*n);
  317. cp Z2[n+1];
  318. for(i=1;i<=n;i++)
  319.     {
  320.       real(Z2[i])=rx[i][p];
  321.       imag(Z2[i])=rx[n+i][p];
  322.     }
  323. cp IF2,U2[n+1],A(1,0);
  324. IF2=A/(Z2[p]+ZF);
  325. for(i=1;i<=n;i++)
  326.      U2[i]=A-Z2[i]*IF2;
  327. out<<"三、节点的自导纳"<<endl;
  328. out<<"节点序号     精确                近似"<<endl;
  329. for(i=1;i<=n;i++)
  330. {
  331.       out<<"    "<<left<<setw(5)<<i;
  332.       out<<left<<setw(20)<<Y1[i];
  333.       out<<Y2[i]<<endl;
  334. }
  335. out<<"四、节点阻抗矩阵的第6列"<<endl;
  336. out<<"序号       精确                   近似"<<endl;
  337. for(i=1;i<=n;i++)
  338. {
  339.       out<<"  "<<left<<setw(4)<<i;
  340.       out<<left<<setw(23)<<Z1[i];
  341.       out<<Z2[i]<<endl;
  342. }
  343. out<<"五、短路电流If"<<endl;
  344. double ifm1,ifm2,mwc,ifj1,ifj2,jwc;
  345. ifm1=abs(IF1);
  346. ifm2=abs(IF2);
  347. mwc=abs(ifm1-ifm2);
  348. ifj1=arg(IF1)*180/PI;
  349. ifj2=arg(IF2)*180/PI;
  350. jwc=abs(ifj1-ifj2);
  351. out<<"        精确         近似      绝对误差"<<endl;
  352. out<<"模值   ";
  353. out<<left<<setw(12)<<ifm1;
  354. out<<left<<setw(12)<<ifm2;
  355. out<<mwc<<endl;
  356. out<<"角度   ";
  357. out<<left<<setw(12)<<ifj1;
  358. out<<left<<setw(12)<<ifj2;
  359. out<<jwc<<endl;
  360. out<<"六、各节点电压的模值"<<endl;
  361. out<<"节点序号    精确       近似       绝对误差"<<endl;
  362. for(i=1;i<=n;i++)
  363. {
  364.       out<<"   "<<left<<setw(7)<<i;
  365.       out<<left<<setw(12)<<abs(U1[i]);
  366.       out<<left<<setw(12)<<abs(U2[i]);
  367.       out<<abs(abs(U2[i])-abs(U1[i]))<<endl;
  368. }
  369. out<<"七、各支路电流"<<endl;
  370. out<<"首节点序号 尾节点序号    支路电流"<<endl;
  371. for(k=0;k<o;k++)
  372.      if(zl[k].t!=3)
  373.      {
  374.            out<<"    "<<left<<setw(11)<<zl[k].i;
  375.            out<<left<<setw(7)<<zl[k].j;
  376.            out<<zl[k].I<<endl;
  377.      }
  378. out<<"八、短路电流容量"<<endl;
  379. out<<"1.46KA"<<endl;
  380. return 0;
  381. }
  382. int jzqn(void*p,void*q,int nn)
  383. {
  384.      double(*a)[nn+1]=(double (*)[nn+1])p;
  385.      double(*c)[nn+1]=(double (*)[nn+1])q;
  386.      double b[nn+1][nn+1]={};
  387.      int i,j,k;
  388.      double e;
  389.      for(i=1;i<=nn;i++)
  390.           for(j=1;j<=nn;j++)
  391.                b[i][j]=a[i][j];
  392.      for(i=0;i<=nn;i++)
  393.           for(j=0;j<=nn;j++)
  394.                c[i][j]=0;
  395.      for(i=1;i<=nn;i++)
  396.             c[i][i]=1;
  397.      for(j=1;j<=nn;j++)
  398.      {
  399.            for(i=j;i<=nn;i++)
  400.                 if(b[i][j]!=0) break;
  401.             if(i!=j)
  402.                   for(k=1;k<=nn;k++)
  403.                    {
  404.                         e=b[i][k];
  405.                         b[i][k]=b[j][k];
  406.                         b[j][k]=e;
  407.                         e=c[i][k];
  408.                         c[i][k]=c[j][k];
  409.                         c[j][k]=e;
  410.                    }
  411.             e=b[j][j];
  412.             if(e!=1)
  413.                  for(k=1;k<=nn;k++)
  414.                   {
  415.                     b[j][k]=b[j][k]/e;
  416.                     c[j][k]=c[j][k]/e;
  417.                   }
  418.             for(i=1;i<=nn;i++)
  419.             {
  420.                   if(i==j || b[i][j]==0) continue;
  421.                   e=b[i][j];
  422.                   for(k=1;k<=nn;k++)
  423.                         {
  424.                               b[i][k]=b[i][k]-e*b[j][k];
  425.                               c[i][k]=c[i][k]-e*c[j][k];
  426.                         }
  427.             }
  428.      }
  429.      return 0;
  430. }
复制代码

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

广播台
特别关注
快速回复 返回顶部 返回列表