模拟退火算法解决路径优化 的源代码.docx
《模拟退火算法解决路径优化 的源代码.docx》由会员分享,可在线阅读,更多相关《模拟退火算法解决路径优化 的源代码.docx(9页珍藏版)》请在冰豆网上搜索。
模拟退火算法解决路径优化的源代码
ÎÒÓÐsimulatedannealingwithmetropolies(MonteCarlo)×öµÄÒ»¸öÏîÄ¿µÄ´úÂ룬ÄãÒª¿´¿´Ã´£¿
voidanneal(intnparam,intnstep,intnstep_per_block,doublet0,
constdouble*param_in,
doublecost_in,double*params_out,double*cost_out){
intnblock;
intstep;
intblock;
intnactive;
intrank;
intn_accepted=0;
inti,j,n;
doublecost_current,cost_trial;
int*param_index;
double*param_current;
double*param_trial;
double*Q;
double*S;
double*u;
double*dp;
double*A;
FILE*fp_log_file;
charfname[FILENAME_MAX];
doubletemp=t0;
doubletempmax=temp;
doubleebar,evar,emin,eta,specific_heat;
doubledelta;
doublechi=0.8;//Annealingschedule
doublechi_s=3.0;//Vanderbilt/Louie'growthfactor'
doublerm;
doubleroot3=sqrt(3.0);
doublep=0.02/sqrt(3.0);//maxsizeofannealingstep
param_current=newdouble[nparam];
param_trial=newdouble[nparam];
cost_current=cost_in;
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
sprintf(fname,"a_%4.4d.log",rank);
fp_log_file=fopen(fname,"a");
if(fp_log_file==(FILE*)NULL)errorMessage("fopen(log)failed\n");
//Workoutthenumberofactiveparameters,andsetupthe
//indextableoftheactiveparameters.
//Notethatthecompletearrayofparameters(param_trial)must
//beusedtoevaluatethecostfunction.
nactive=0;
for(n=0;nparam_current[n]=param_in[n];
param_trial[n]=param_in[n];
if(P.is_active[n])nactive++;
}
param_index=newint[nactive];
i=0;
for(n=0;nif(P.is_active[n])param_index[i++]=n;
}
//InitialisethestepdistributionmatrixQ_ij
Q=newdouble[nactive*nactive];
S=newdouble[nactive*nactive];
u=newdouble[nactive];
dp=newdouble[nactive];
A=newdouble[nactive];
double*Qtmp;
Qtmp=newdouble[nactive*nactive];
for(i=0;ifor(j=0;jdelta=(i==j);
Q[i*nactive+j]=p*delta*param_current[param_index[j]];
}
}
//carryoutannealingpoints
nblock=nstep/nstep_per_block;
rm=1.0/(double)nstep_per_block;
for(block=0;block
//Setthescheduleforthisblock,andinitialiseblockwisequantities.
//Wealsoensurethestepdistributionmatrixisdiagonal.
temp=chi*temp;
for(i=0;iA[i]=0.0;
for(j=0;jS[i*nactive+j]=0.0;
delta=(i==j);
Q[i*nactive+j]*=delta;
}
}
ebar=0.0;
evar=0.0;
emin=cost_current;
for(i=0;iprintf("Step:
%d%g\n",i,Q[i*nactive+i]);
}
for(step=0;step
//Settherandomvectoru,andcomputethestepsizedp
for(i=0;iu[i]=root3*(r_uniform()*2.0-1.0);
}
for(i=0;idp[i]=0.0;
for(j=0;jdp[i]+=Q[i*nactive+j]*u[j];
}
}
for(i=0;i
n=param_index[i];
param_trial[n]=param_current[n]+dp[i];
if(param_trial[n]
if(param_trial[n]>P.max[n])param_trial[n]=P.max[n];
}
//calculatenewcostfunctionscore
p_model->setParameters(param_trial);
cost_trial=p_costWild->getCost();
cost_trial+=p_costLHY->getCost();
cost_trial+=p_costTOC1->getCost();
cost_trial+=p_costAPRR->getCost();
//Metropolis
delta=cost_trial-cost_current;
if(delta<0.0||r_uniform()
for(n=0;nparam_current[n]=param_trial[n];
}
cost_current=cost_trial;
++n_accepted;
}
//'Energy'statistics
ebar+=cost_current;
evar+=cost_current*cost_current;
if(cost_current
//Pertimesteplog
fprintf(fp_log_file,"%6d%6d%10.4f%10.4f%10.4f%10.4f\n",
block,step,temp,
cost_current,cost_trial,
(float)n_accepted/(float)(block*nstep_per_block+step));
//Accumulateaverage,measuredcovariance
for(i=0;iA[i]+=param_current[param_index[i]];
for(j=0;jS[i*nactive+j]
+=param_current[param_index[i]]*param_current[param_index[j]];
}
}
/*Nextstep*/
}
//Setthepreviousblockaverageandmeasuredcovariance
for(i=0;iA[i]=rm*A[i];
}
for(i=0;ifor(j=0;jS[i*nactive+j]=rm*S[i*nactive+j]-A[i]*A[j];
if(i==j)printf("Average:
%d%g%g\n",i,A[i],S[i*nactive+j]);
//Settheconvarienceforthenextiterations=6chi_sS/M
S[i*nactive+j]=6.0*chi_s*rm*S[i*nactive+j];
}
}
//Resetthestepdistributionmatrixforthenextblock
i=do_cholesky(nactive,S,Q);
j=test_cholesky(nactive,S,Q);
printf("Cholesky%d%d\n",i,j);
//Blockstatistics
ebar=rm*ebar;
evar=rm*evar;
specific_heat=(evar-ebar*ebar)/temp*temp;
eta=(ebar-emin)/ebar;
fprintf(fp_log_file,"%d%d%f%f%f%f%f%f\n",
block,nstep_per_block,temp,ebar,evar,emin,
specific_heat,eta);
/*Nextblock*/
}
*cost_out=cost_current;
for(n=0;nparams_out[n]=param_current[n];
}
fclose(fp_log_file);
deleteparam_index;
deleteparam_current;
deleteparam_trial;
deleteQ;
deleteu;
deletedp;
deleteS;
deleteA;
return;
}