牛顿拉夫逊算法进行潮流计算，数据处理程序，方便处理实验数据，并取得数学表达式

#include<iostream> #include<complex> #include<fstream> using namespace std; #include<cmath> #include<iomanip> #define pi 3.1415926 complex<float> conj(complex<float>); float Q(complex<float> [15][1],complex<float> [15][15],int); float AM(complex<float>); void main() { float realY,imagY; float PL[15][1],PG[15][1],QL[15][1],QG[15][1],value; complex<float> Y[15][15]; int i,j; ifstream YDATA1("YDATAC1_10.txt",ios::in); if(!YDATA1){ cout<<"Adata cannot be opened!"<<endl; exit(1); } else{ for(i=1;i<15;i++) for(j=1;j<=10;j++){ YDATA1>>realY>>imagY; Y[i][j].real(realY); Y[i][j].imag(imagY); } } ifstream YDATA2("YDATAC11_14.txt",ios::in); if(!YDATA2){ cout<<"Adata cannot be opened!"<<endl; exit(1); } else{ for(i=1;i<=14;i++) for(j=11;j<=14;j++){ YDATA2>>realY>>imagY; Y[i][j].real(realY); Y[i][j].imag(imagY); } } ifstream SDATA("SDATA.txt",ios::in); if(!SDATA){ cout<<"Sdata cannot be opened!"<<endl; exit(1); } else{ for(i=1;i<=14;i++){ SDATA>>value; PG[i][0]=value; } for(i=1;i<=14;i++){ SDATA>>value; QG[i][0]=value; } for(i=1;i<=14;i++){ SDATA>>value; PL[i][0]=value; } for(i=1;i<=14;i++){ SDATA>>value; QL[i][0]=value; } } float Psp[15][1],Qsp[15][1]; for(i=1;i<=14;i++){ Psp[i][0]=PG[i][0]-PL[i][0]; Qsp[i][0]=QG[i][0]-QL[i][0]; } complex<float> U[15][1]; U[1][0].real(1.06); U[1][0].imag(0); for(i=2;i<=14;i++){ U[i][0].real(1); U[i][0].imag(0); } float dU=0; complex<float> temp; ofstream IEEE14("IEEE14.txt",ios::out); int k=0; do{ IEEE14<<left<<setw(10)<<k; for(i=1;i<=14;i++){ IEEE14<<left<<setw(10)<<AM(U[i][0])<<left<<setw(10)<<atan(imag(U[i][0])/real(U[i][0]))/(2*pi)*360; } IEEE14<<endl; dU=0; for(i=2;i<=14;i++){ if(i==2){ temp=U[i][0]; complex<float> Ssp(Psp[i][0],-Qsp[i][0]); U[i][0]=Ssp/conj(U[i][0])/Y[i][i]; for(j=1;j<=14;j++){ if(j==i) continue; else U[i][0]-=Y[i][j]/Y[i][i]*U[j][0]; } U[i][0]*=1.045/AM(U[i][0]); temp=U[i][0]-temp; dU+=AM(temp); if(Q(U,Y,i)>0.5) Qsp[i][0]=0.5; else if (Q(U,Y,i)<-0.4) Qsp[i][0]=-0.4; else Qsp[i][0]=Q(U,Y,i); } else if(i==3){ temp=U[i][0]; complex<float> Ssp(Psp[i][0],-Qsp[i][0]); U[i][0]=Ssp/conj(U[i][0])/Y[i][i]; for(j=1;j<=14;j++){ if(j==i) continue; else U[i][0]-=Y[i][j]/Y[i][i]*U[j][0]; } U[i][0]*=1.01/AM(U[i][0]); temp=U[i][0]-temp; dU+=AM(temp); if(Q(U,Y,i)>0.4) Qsp[i][0]=0.4; else if (Q(U,Y,i)<0) Qsp[i][0]=0; else Qsp[i][0]=Q(U,Y,i); } else if(i==6){ temp=U[i][0]; complex<float> Ssp(Psp[i][0],-Qsp[i][0]); U[i][0]=Ssp/conj(U[i][0])/Y[i][i]; for(j=1;j<=14;j++){ if(j==i) continue; else U[i][0]-=Y[i][j]/Y[i][i]*U[j][0]; } U[i][0]*=1.07/AM(U[i][0]); temp=U[i][0]-temp; dU+=AM(temp); if(Q(U,Y,i)>0.24) Qsp[i][0]=0.24; else if (Q(U,Y,i)<-0.06) Qsp[i][0]=-0.06; else Qsp[i][0]=Q(U,Y,i); } else if(i==8){ temp=U[i][0]; complex<float> Ssp(Psp[i][0],-Qsp[i][0]); U[i][0]=Ssp/conj(U[i][0])/Y[i][i]; for(j=1;j<=14;j++){ if(j==i) continue; else U[i][0]-=Y[i][j]/Y[i][i]*U[j][0]; } U[i][0]*=1.09/AM(U[i][0]); temp=U[i][0]-temp; dU+=AM(temp); if(Q(U,Y,i)>0.24) Qsp[i][0]=0.24; else if (Q(U,Y,i)<-0.06) Qsp[i][0]=-0.06; else Qsp[i][0]=Q(U,Y,i); } else{ temp=U[i][0]; complex<float> Ssp(Psp[i][0],-Qsp[i][0]); U[i][0]=Ssp/conj(U[i][0])/Y[i][i]; for(j=1;j<=14;j++){ if(j==i) continue; else U[i][0]-=Y[i][j]/Y[i][i]*U[j][0]; } temp=U[i][0]-temp; dU+=AM(temp); } } k++; }while(dU>pow(10,-5)); } complex<float> conj(complex<float> x) { x.real(real(x)); x.imag(-imag(x)); return x; } float Q(complex<float> U[15][1],complex<float> Y[15][15],int p){ complex<float> S(0,0); S=conj(U[p][0])*Y[p][p]*U[p][0]; for(int J=1;J<=14;J++){ if(J==p) continue; else S+=conj(U[p][0])*Y[p][J]*U[J][0]; } return -imag(S); } float AM(complex<float> x){ return sqrt(pow(real(x),2)+pow(imag(x),2)); }