?? 分水系統模型.cpp
字號:
#include <stdio.h>
#include <math.h>
#define n 6
FILE *fp1,*fp2,*fp3;
double min(double a,double b)
{ if(a>b) return b; else return a; }
void eq(float A[][n],float B[][n]) //A[i][j]<--B[i][j]
{ int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++) A[i][j]=B[i][j]; }
void uni_eq(float C[][n],float D[][n]) //C[i][j]<--D[i][j]
{ int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++) C[i][j]=(-1)*D[i][j];}
float product(float E[][n],float L[][n]) //return E[][]*L[][]
{ float sum=0;int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++) sum=sum+E[i][j]*L[i][j];
return sum; }
void add(float c1,float M[][n],float c2,float N[][n],float R[][n]) //R[][]=c1*M[][]+c2*N[][]
{ int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++) R[i][j]=c1*M[i][j]+c2*N[i][j] ; }
void Add_line(float Line_a[],float Line_x[][n],float Line_c[][n])
{ int i,j;
for(i=0;i<n;i++) Line_a[i]=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
Line_a[i]=Line_a[i]+Line_x[i][j]*Line_c[i][j];
} //Line_x[i]為左行c*X和
void Add_x_line(float Line_a[],float Line_x[][n])
{ int i,j;
for(i=0;i<n;i++)
Line_a[i]=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
Line_a[i]=Line_a[i]+Line_x[i][j];
} //Line_x[i]為左行X和
void ceshihanghe(float Line_a[],float Line_x[][n])
{ int i,j;
for(i=0;i<n;i++)
Line_a[i]=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
Line_a[i]=Line_a[i]+Line_x[i][j];
fp3=fopen("ls_共厄梯度解目標.txt","a");
for(i=0;i<n;i++)
fprintf(fp3,"\t行和:%f\t",Line_a[i]);
fprintf(fp3,"over\n");
fclose(fp3); } //Line_x[i]為左行X和--------------------測試
void Add_list(float List_a[],float List_x[][n],float List_c[][n])
{int i,j;
for(j=0;j<n;j++)
{ List_a[j]=0;
for(i=0;i<n;i++)
List_a[j]=List_a[j]+List_x[i][j]*List_c[i][j]; } } //List_x[i]為列c*X和
void Add_x_list(float List_a[],float List_x[][n])
{int i,j;
for(j=0;j<n;j++)
{ List_a[j]=0;
for(i=0;i<n;i++)
List_a[j]=List_a[j]+List_x[i][j]; } } //List_x[i]為列X和
void ceshiliehe(float List_a[],float List_x[][n])
{ int i,j;
for(j=0;j<n;j++)
{ List_a[j]=0;
for(i=0;i<n;i++)
List_a[j]=List_a[j]+List_x[i][j]; }
fp3=fopen("ls_共厄梯度解目標.txt","a");
for(i=0;i<n;i++)
fprintf(fp3,"\t列和:%f\t",List_a[i]);
fprintf(fp3,"over\n");
fclose(fp3); } //Line_x[i]為列c*X和--------------------測試
double F_min(float F_C[][n],float F_X[][n],float F_A[],float F_B[],float F_mu) //F_min()無約束目標函數 , F_A[]=A[i] //f(x)
{ int i,j,k;
double f=0;
float a[n],b[n];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{if(F_X[i][j]>=0) k=0;
else k=1;
f=f+k*pow(F_X[i][j],2); //解決階躍函數,Xij平方和，解決x為正
}
Add_x_line(a,F_X); //行和a[i]
Add_x_list(b,F_X); //列和b[i]
for(i=0;i<n;i++) f=f+pow(a[i]-F_A[i],2)+pow(b[i]-F_B[i],2); //f追加等式約束
f=f*F_mu; //引入罰因子F_mu
for(i=0;i<n;i++)
for(j=0;j<n;j++) f=f+F_C[i][j]*F_X[i][j]; //追加原約束問題的目標函數c*x
return f;
}
double rf(float F_C[][n],float F_X[][n],float F_A[],float F_B[]) //F_min()無約束目標函數 , F_A[]=A[i] //f(x)
{ int i,j,k;
double f=0;
float a[n],b[n];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{if(F_X[i][j]>=0) k=0;
else k=1;
f=f+k*pow(F_X[i][j],2); //解決階躍函數,Xij平方和，解決x為正
}
Add_x_line(a,F_X); //行和a[i]
Add_x_list(b,F_X); //列和b[i]
for(i=0;i<n;i++) f=f+pow(a[i]-F_A[i],2)+pow(b[i]-F_B[i],2);//f追加等式約束
fp3=fopen("ls_共厄梯度解目標.txt","a");
fprintf(fp3,"\n+++++++++++++++++++++++++++++++++++++++++++++++rf(a)=%-17.5e+++++++++++++++++++++++++++++++++++++\n",f);
fclose(fp3); //測試罰函數
return f;
}
void F_DF(float F_g[][n],float F_C[][n],float F_X[][n],float F_A[],float F_B[],float F_mu) //無約束目標函數F梯度 F_g[][]存放 , F_B[]=B[i]
{ int i,j,k;
float f=0,a[n],b[n];
Add_x_line(a,F_X); //行和a[i]
Add_x_list(b,F_X); //列和b[i]
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{ if(F_X[i][j]>=0) k=0;
else k=1; //解決階躍函數
F_g[i][j]=F_C[i][j]+2*F_mu*(a[i]-F_A[i]+b[j]-F_B[j]+k*F_X[i][j]); } //梯度分量分配結束
}
void ceshitidu(float F_g[][n],float F_C[][n],float F_X[][n],float F_A[],float F_B[],float F_mu) //無約束目標函數F梯度 F_g[][]存放 , F_B[]=B[i]
{ int i,j,k;
float f=0,a[n],b[n];
Add_x_line(a,F_X); //行和a[i]
Add_x_list(b,F_X); //列和b[i]
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{ if(F_X[i][j]>=-1e-9) k=0;
else k=1; //解決階躍函數
F_g[i][j]=F_C[i][j]+2*F_mu*(a[i]-F_A[i]+b[j]-F_B[j]+k*F_X[i][j]); } //梯度分量分配結束----------------------測試
fp3=fopen("ls_共厄梯度解目標.txt","a");
for(i=0;i<n;i++)
{for(j=0;j<n;j++)
fprintf(fp3,"DF%-d%-d=%f",i,j,F_g[i][j]);
fprintf(fp3,"over\n");}
fclose(fp3);
}
void ls(float NP[][n],float F_X0[][n],float F_C[][n],float F_X[][n],float F_A[],float F_B[],float F_mu) //RX[][]=ls(MX[][]F_X,NP[][]_-F_g,RX[][]_NEW_X_Point)
{ float p=0.1,B=2,v=0.4,gg1[n][n],gg2[n][n],g1=0,g2=0;
int i,k=0,j;
float f1,f2;
double t=1,a=0,b=1e+6;
f1=F_min(F_C,F_X0,F_A,F_B,F_mu); //f1 為無約束函數值
F_DF(gg1,F_C,F_X0,F_A,F_B,F_mu);
g1=product(gg1,NP);
for(k=0;k<500;k++)
{ add(1,F_X0,t,NP,F_X); //MX has been changed
f2=F_min(F_C,F_X,F_A,F_B,F_mu);
F_DF(gg2,F_C,F_X,F_A,F_B,F_mu); //g=DF(x)
g2=product(gg2,NP);
//fp3=fopen("ls_共厄梯度解目標.txt","a");
// for(i=0;i<n;i++)
// {for(j=0;j<n;j++)
// fprintf(fp3,"F_X(%-d%-d)=%f",i,j,F_X[i][j]);
// fprintf(fp3,"over\n");}
// fprintf(fp3,"over\n");
// fclose(fp3);
//fp3=fopen("ls_共厄梯度解目標.txt","a");
// fprintf(fp3,"g2=%-17.20eg1=%-17.5ef1=%-17.15ef2=%-17.15et=%.30f",g2,g1,f1,f2,t);
// fprintf(fp3,"\n");
//fclose(fp3);
if(g2-v*g1<-1e-10&&fabs(a-b)>1e-20)
{ a=t;t=min((a+b)/2.0,2.0*a);
//fp3=fopen("ls_共厄梯度解目標.txt","a");
//fprintf(fp3,"條件一a=%-17.20e\tb=%-17.20e\tt=%-17.25e\tf2=%-17.15e",a,b,t,f2);
//fprintf(fp3,"\n");
// fclose(fp3);
continue;}
else if(f1-f2<(-1.0)*p*t*g1&&fabs(a-b)>1e-20)
{b=t;
t=0.5*(a+b);
//fp3=fopen("ls_共厄梯度解目標.txt","a");
// fprintf(fp3,"條件二a=%-17.20e\tb=%-17.20e\tt=%-17.25e\tf2=%-17.15e",a,b,t,f2);
// fprintf(fp3,"\n"); fclose(fp3);
continue; }
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -