求讲解

陀飞轮

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. }
     

复制代码