wolfe.cpp

最优化实验中的 WOLF 算法 黄金分割法 平分法程序

#include<iostream.h>
#include<string.h>
#include<math.h>
const int N=2;
int Ncalf;
int Ncalg;
int iprint;
int info;
double F;
double work[N];
double G[N];
double alpha;
double c1=0.1,c2=0.5;
int Ncalfn;
bool A21;

double DDOT1(double *x, double *y)
{
	long double ddot1=0.0;
	for(int i=0;i<N;i++)
		ddot1+=x[i]*y[i];
	
	return ddot1;
}
void CALF( double *x)
{
	F=0.0;
	for(int i=0;i<N-1;i++)
		F+=0.01*(x[i+1]-x[i]*x[i])*(x[i+1]-x[i]*x[i])+(1-x[i])*(1-x[i]);
}
void CALG(double*x )
{
	
	G[0]=-400.0*x[0]*(x[1]-x[0]*x[0])-2.0*(1.0-x[0]);
	G[N-1]=200.0*(x[N-1]-x[N-2]*x[N-2]);
	
}
int inexls ( double *x,double *D)
{
	 double T1,T2,T3,T4;
	 double dd,xx;
	 double DGO,DG,DG2;
	 double F1,F2;
	 double DG1;
     double A1,A2;
	 double rati01,rati02;
	 double TRY;
	int i,s=0;
	long double F0;
	while(1)
	{
		int Nlsmax=50;
		info=1;
		DGO=DDOT1(D,G);
		if(DGO>=0.0)goto a1100;	
		F0=F;
		for(i=0;i<N;i++)work[i]=x[i];
		A1=0.0;
		F1=F;
		DG1=DGO;
		A2=-0.1;
		A21=false;
		rati01=1.0;
		rati02=0.1;
		alpha=1.0;
        if(iprint!=0) goto a40;
		else goto a2010;
a40:
	for( i=0;i<N;i++)
		x[i]=work[i]+alpha*D[i];
	CALF(x);
	Ncalfn+=1;
	if(Ncalfn>Nlsmax) goto a1300;
	if(F<=F0+c1*alpha*DGO)goto a200;
	rati01*=0.1;
	if(DGO/(1.0+fabs(F0))<-1e-16)goto a80;

	dd=DDOT1(D,D);
	xx=DDOT1(x,x);
	if(alpha*sqrt(dd)/(1.0+sqrt(xx))<=1.0e-14)goto a1400;
a80:
	A2=alpha;
	F2=F;
	A21=false;
a90:
	
	TRY=((F2-F1)/(A2-A1)-DG1)/(A2-A1);
	TRY=-DG1/TRY/2.0;
	goto a240;

a200:CALG(x);
	 Ncalg+=1;
	 DG=DDOT1(D,G);
	 if(fabs(DG)<=-c2*DGO)goto a1500;
	 if(DG>0.0)goto a210;
	 if(A2<0.0)goto a260;
	 rati02*=0.1;
	 A1=alpha;
	 DG1=DG;
	 rati01=0.1;
	 if(!A21)goto a90;
	 goto a220;
a210:rati01=0.1;
	 A2=alpha;
	 F2=F;
	 DG2=DG;
	 A21=true;
a220:
	 T1=(F2-F1)/(A2-A1)-DG1;
	 T2=DG2-DG1;
	 T3=3.0*T1-T2;
	 T4=T2-T1-T1;
	 TRY=T3*T3-3.0*T4*DG1;
	 if(TRY<0.0)goto a230;
	 if(sqrt(TRY)<=T3)goto a230;
	 TRY=-DG1*(A2-A1/sqrt(TRY)+T3);
	 goto a240;
a230:
	TRY=((F2-F1)/(A2-A1)-DG1)/(A2-A1);
	TRY=-DG1/TRY/2.0;
a240:
	if(TRY<rati01*(A2-A1))TRY=rati01*(A2-A1);
	if(TRY>(1.0-rati02)*(A2-A1))TRY=(1.0-rati02)*(A2-A1);
	alpha=TRY+A1;
	goto a40;
a260:T1=(F-F1)/(alpha-A1)-DG;
	 T2=DG1-DG;
	 T3=3.0*T1-T2;
	 T4=T2-T1-T1;
	 TRY=T3*T3-3.0*T4*DG;
	 if(TRY<0.0)goto a275;
	 if(sqrt(TRY)<=T3)goto a275;
	 TRY=-DG*(alpha-A1)/(sqrt(TRY)-T3);
	 if(TRY<=0.1*(alpha-A1))goto a280;
a275:TRY=0.1*(alpha-A1);
a280:if(A2<0.0)goto a290;
	 if(TRY<0.1*(A2-A1))TRY=0.1*(A2-A1);
	 if(TRY<0.9*(A2-A1))TRY=0.9*(A2-A1);
a290:A1=alpha;
	 F1=F;
	 DG1=DG;
	 alpha=A1+TRY;
	 goto a40;
a1100:info=1;
	  goto a1500;
a1300:info=2;
	  goto a1500;
a1400:CALG(x);
	  Ncalg+=1;
	  goto a2080;
a1500:Ncalf+=Ncalfn;
	  return 1;
a2010:cout<<"'enter line search"<<endl;
	  cout<<"F="<<F<<endl;
	  cout<<"DGO="<<DGO<<endl;
	  goto end;
a2080:cout<<"!!step is too small"<<endl;
end:;
	s++;
	if(s==5)break;
}
return 1;
}

int main()
{

	double X[N],G[N],D[N],WORK[N];
	for(int i=1;i<=N;i++)
	{
		if(i%2==0)X[i-1]=-1.2;
		else X[i-1]=1.0;
	}
	
	 CALF(X);
	 Ncalf=1;
	
	
	CALG(X);
    Ncalg=1;

	for(i=0;i<N;i++)
		D[i]=-G[i];
	cout<<" F= "<<F<<endl;
	for(i=0;i<N;i++)
		cout<<" G["<<i<<"] "<<G[i]<<'\t';
	cout<<endl;

	inexls (X,D);
	cout<<"alpha="<<alpha<<endl;
	for(i=0;i<N;i++)
		cout<<" G["<<i<<"] "<<G[i]<<'\t';
	cout<<endl;
	cout<<"new F="<<F<<endl;
	cout<<"Ncalf="<<Ncalf<<endl;
	cout<<"Ncalg="<<Ncalg<<endl;
	cout<<"!!press enter enter key to end!!"<<endl;
	cin.get();
	return 0;
}