?? 分水系統模型.cpp
字號:
break; } //break必須有,達到前兩條件后調處,k不再自加
}
void main()
{
int i=0,j=0,k,H=0,I=1,J=0,II;
float C[n][n],A[n],B[n],h1,h2,h3,FDC=10;
float mu=100;
float X0[n][n],X[n][n];
float a[n],b[n];
double f,f0;
float t,sum=0,p[n][n],p0[n][n],g[n][n],g0[n][n],rfa,e=1e-3,e1,e2,e3;
FILE *fp1,*fp2,*fp3;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
X0[i][j]=1; // X 被初始化 均為1
fp1=fopen("參數庫.txt","w+");
fprintf(fp1,"2.5 1.60 4.55 2.18 2.18 10000 589\n"); //題目參量矩陣
fprintf(fp1,"2.08 0 2.56 1.40 0.25 10000 49\n");
fprintf(fp1,"0.95 0 1.53 0.15 0.12 10000 295\n");
fprintf(fp1,"2.46 0 2.94 1.66 1.28 10000 48\n");
fprintf(fp1,"1.11 0.099 1.87 0.42 0.42 10000 5.4\n");
fprintf(fp1,"0.80 0 0.74 0 0 10000 18\n");
//fprintf(fp1,"10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 0\n");
//fprintf(fp1,"10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 0\n");
//fprintf(fp1,"10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 0\n");
//fprintf(fp1,"10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 0\n");
fprintf(fp1,"486 97 90 295 36.4 0\n");//第11行:B[j]
//fprintf(fp1,"1 2 3 4 5 6 7 8 9 0 10\n"); //測試矩陣
//fprintf(fp1,"11 12 13 14 15 16 17 18 19 0 10\n");
//fprintf(fp1,"21 22 23 24 25 26 27 28 29 0 10\n");
//fprintf(fp1,"31 32 33 34 35 36 37 38 39 0 10\n");
//fprintf(fp1,"41 42 43 44 45 46 47 48 49 0 10\n");
//fprintf(fp1,"51 52 53 54 55 56 57 58 59 0 10\n");
//fprintf(fp1,"61 62 63 64 65 66 67 68 69 0 10\n");
//fprintf(fp1,"71 72 73 74 75 76 77 78 79 0 10\n");
//fprintf(fp1,"81 82 83 84 85 86 87 88 89 0 10\n");
//fprintf(fp1,"91 92 93 94 95 96 97 98 99 0 10\n");
//fprintf(fp1,"10 10 10 10 10 10 10 10 10 10\n");//第11行:B[j]
fclose(fp1);
fp1=fopen("參數庫.txt","r");
for(i=0;i<n;i++)
{ for(j=0;j<n;j++) fscanf(fp1,"%f",&C[i][j]);//輸入C[i][j]
fscanf(fp1,"%f",&A[i]); } //輸入A[i][j]——行和
for(i=0;i<n;i++) fscanf(fp1,"%f",&B[i]); //輸入B[i][j]——列和
fclose(fp1);
fp2=fopen("結果庫.txt","w+");//輸出文件,格式(規整)化系數矩陣
fprintf(fp2,"A====================================A\n");
for(i=0;i<n;i++)
{ for(j=0;j<n;j++)
fprintf(fp2,"%5.3f\t",C[i][j]);
fprintf(fp2,"%-5.3f\t",A[i]);
fprintf(fp2,"\n"); }
fprintf(fp2,"\n");
for(i=0;i<n;i++)
fprintf(fp2,"%8.3f",B[i]);
fprintf(fp2,"\nA===========A===========A============A\n");
fclose(fp2);
for(II=0;II<20;II++) //最外層_起始
{
k=1; J=0; j=0; //進入外罰法======================================
fp3=fopen("ls_共厄梯度解目標.txt","w+");
fprintf(fp3,"START=========================================================================WRITE\n",I++);
fclose(fp3); ++I;
ceshihanghe(a,X0);
ceshiliehe(a,X0);
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\t=====mu=%-17.e===========\n",mu);
fclose(fp3);
eq(X,X0);
f=F_min(C,X0,A,B,mu);
F_DF(g,C,X0,A,B,mu) ;
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\nFx=%-17.5e\n",F_min(C,X,A,B,mu));
rf(C,X,A,B);
fclose(fp3);
do
{ k=0; uni_eq(p0,g);
for(;H!=1&&k<300;k++)
{ eq(X0,X);
f0=f;
eq(g0,g);
rf(C,X,A,B);
ls(p0,X0,C,X,A,B,mu);
//x0 has been put into x
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"ls%d=====================================================================ls%d\n",J,J);
rf(C,X,A,B);
for(i=0;i<n;i++)
{ for(j=0;j<n;j++)
fprintf(fp3,"x%-d%d=%-17.5e",i,j,X[i][j]);
fprintf(fp3,"\n"); }
fclose(fp3);
++J;
f=F_min(C,X,A,B,mu);
F_DF(g,C,X,A,B,mu) ;
ceshihanghe(a,X0);
ceshiliehe(a,X0);
e1=e2=1e-5;e3=1e-4; //H order
h1=h2=h3=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{ h1=h1+fabs(g[i][j]);
h2=h2+fabs(X[i][j]-X0[i][j]);
h3=h3+fabs(X0[i][j]); }
h2=h2/(h3+1);
double(h3)=fabs(F_min(C,X,A,B,mu)-F_min(C,X0,A,B,mu))/(fabs(F_min(C,X0,A,B,mu))+1);
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\n\nh1=%-17.eh2=%-17.eh3=%-17.e",h1,h2,h3);
fprintf(fp3,"\ne3=%-17.ee1=%-17.ee2=%-17.e\n",e3,e1,e2);
fclose(fp3);
if(h1<e3&&h2<=e1&&h3<=e2)
{ H=1;
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\n=======================H order has been passed.=======================\n");
for(i=0;i<n;i++)
for(j=0;j<n;j++)
fprintf(fp3,"x%-d%-d=%-17.5e",i,j,X[i][j]);
fprintf(fp3,"\nH=%-18dFx=%-17.5ek=%-18dJ=%-17d\n",H,F_min(C,X,A,B,mu),k,J);
fprintf(fp3,"\nh1=%-17.eh2=%-17.eh3=%-17.10e\n",h1,h2,h3);
fprintf(fp3,"\ne3=%-17.ee1=%-17.ee2=%-17.e\n",e3,e1,e2);
fclose(fp3);
break; }
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\n");
fprintf(fp3,"H=%-18dFx=%-17.5ek=%-18dJ=%-17d\n",H,F_min(C,X,A,B,mu),k,J);
fprintf(fp3,"\t\tJump over the H order and go down.\n");
fclose(fp3);
if(k==n*n) { fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\n======================從k=n跳出====================");
fclose(fp3);
break;} //break from 'for' and continue 'do while'
rfa=product(g,g)/product(g0,g0);
add(-1,g,rfa,p0,p); //p=-g+a*p0
if(fabs(product(p,g))<=e) {fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\n======================從||p*g||<=e跳出====================");
fclose(fp3);
break; }
else if(product(p,g)<-1*e)
eq(p0,p);
else uni_eq(p0,p);
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\n======================從k=k+1回for循環====================\n\n");
fclose(fp3); } //end for
J++;
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\n======================從while返回,繼續while循環====================\n\n");
fclose(fp3);
}while(H==0&&J<500);
//output on the square-eye
for(i=0;i<n;i++)
{ for(j=0;j<n;j++)
printf("x%-d%-d=%-17.5e",i,j,X[i][j]);
printf("\nFx=%.5f\tk=%d\tj=%d\n",F_min(C,X,A,B,mu),k,j); }
if(log10(mu)+log10(rf(C,X,A,B))<-4) {fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\t=============計算終止=============mu=%-17.e======計算終止==========II=%d==============\n",mu,II);
fclose(fp3); break;}
else { fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\t===計算沒有終止==mu=%-17.e===========\n",mu);
fclose(fp3);
mu=mu*10; eq(X0,X); printf("mu=%17.e\n",mu);J=0;}
} //最外層_終止
//改F_MIN為double
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -