📄 cpp1.cpp

📁 最优化解法的运行程序大家帮忙支持下我自己做的哦！！谢谢了
💻 CPP
字号:
#include<iostream.h>
#include<math.h>
#include<time.h>
 
#define N1 40
//cF54                             DIXMAANE                              
//*                             
//*                    Initial point x0=[2.,2.,2...,2.].
//*

void CALF(int n,double x[N1],double &f){
        double alpha,beta,gamma,delta;                      
        int  k1,k2,k3,k4;
    
        alpha=1.0;
        beta=0.0;
        gamma=0.1250;
        delta=0.1250;
                      
        k1=1;
        k2=0;
        k3=0;
        k4=1;
           
        int m = n/3;
        
        f=1.0;
        for(int i=0;i<n;i++)
          f=f+alpha*x[i]*x[i]*pow((float(i+1)/float(n)),k1);
   
    
        for(i=0;i<n-1;i++)
          f=f+beta*x[i]*x[i]*(x[i+1]+x[i+1]*x[i+1])*pow((float(i+1)/float(n)),k2);
   
        
        for(i=0;i<2*m;i++)
          f=f+gamma*x[i]*x[i]*(pow(x[i+m],4))*pow((float(i+1)/float(n)),k3);
    
        
       for(i=0;i<m;i++)
          f=f+delta*x[i]*x[i+2*m]*pow((float(i+1)/float(n)),k4);
}
//*-- Gradient    
void CALG(int n,double x[N1],double g[N1]){
        double alpha,beta,gamma,delta;                       
        int  k1,k2,k3,k4;
    
        alpha=1.0;
        beta=0.0;
        gamma=0.1250;
        delta=0.1250;
                      
        k1=1;
        k2=0;
        k3=0;
        k4=1;
           
        int m=n/3;

        g[0]=2.0*alpha*x[0]*pow((float(1)/float(n)),k1)+2.0*beta*x[0]*pow((x[0]+x[1]*x[1]),2)*pow((float(1)/float(n)),k2)+2.0*gamma*x[0]*pow(x[m],4)*pow((float(1)/float(n)),k3)+delta*x[2*m]*pow((float(1)/float(n)),k4);
       for(int i=1;i<m;i++)
          g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i+1)/float(n)),k2)+2.0*gamma*x[i]*pow(x[i+m],4)*pow((float(i+1)/float(n)),k3)+delta*x[i+2*m]*pow((float(i+1)/float(n)),k4);
        
       for(i=m;i<2*m;i++)
          g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i+1)/float(n)),k2)+2.0*gamma*x[i]*pow(x[i+m],4)*pow((float(i+1)/float(n)),k3)+4.0*gamma*pow(x[i-m],2)*pow(x[i],3)*pow((float(i-m+1)/float(n)),k3);
        
       for(i=2*m;i<n-1;i++)
          g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i)/float(n)),k2)+4.0*gamma*pow(x[i-m],2)*pow(x[i],3)* pow((float(i-m+1)/float(n)),k3)+delta*x[i-2*m]*pow((float(i-2*m+1)/float(n)),k4);
        
        
        g[n-1]=2.0*alpha*x[n-1]+2.0*beta*x[n-2]*x[n-2]*(x[n-1]+x[n-1]*x[n-1])*(1.0+2.0*x[n-1])*pow((float(n-1)/float(n)),k2)+4.0*gamma*x[2*m-1]*x[2*m-1]*x[n-1]*x[n-1]*x[n-1]*pow((float(2*m)/float(n)),k3)+delta*x[m-1]*pow((float(m)/float(n)),k4);
}




int conclude1(double G[N1])
{double a=0.00;
 for(int i=0;i<N1;i++)
	{a+=G[i]*G[i];}
 if(sqrt(a)>=0.0001)    {return 1;}
 else return 0;
}

double DDOT1(double p[N1],double Q[N1])
{double t=0.00;
 for(int i=0;i<N1;i++)
	{t+=p[i]*Q[i];}
 return t;
}



void main()
{clock_t start,end;
 start=clock();
 double k=1,kesi=0.0001;
double beta,m=0.00,n=0.00;
double X[N1],G[N1],WORK[N1],D[N1],P[N1],G1[N1]  ; 
double F,ALPHA,IPRINT=1,C1=0.01,C2=0.05; 
int NLSMAX ;
	  double FO ;
      int I; 
      double DGO; 
      double A1; 
      double F1; 
      double DG1; 
      double A2; 
      double A21; 
      double RATIO1; 
      double RATIO2; 
      int NCALFN,NCALG,NCALF; 
      double DD,XX;
	  double TRY,F2,DG,DG2;
	  double T1,T2,T3,T4;
      double  INFO=1; 
for(int i=0;i<N1;i++)
 {X[i]=2;}
CALG(N1,X,G);
for(i=0;i<N1;i++)
	{G1[i]=G[i];} 

for(i=0;i<N1;i++)
	{P[i]=-G[i];} 
//................................................................................
step1: CALG(N1,X,G);
	   //cout<<x[1]<<'\t';
	   for(i=0;i<N1;i++)
	   {D[i]=-G[i];}
step2: if(conclude1(G)){
       if(k==1)
	   {beta=0;
		for(i=0;i<N1;i++)
		{ G1[i]=G[i];} 
	   }
	   else {for(i=0;i<N1;i++)
		     {m+=G[i]*G[i];
	          n+=G1[i]*G1[i];
	   }
	   beta=m/n;
	   for(i=0;i<N1;i++)
	   { G1[i]=G[i];}
	   
	   }
	 
  for(i=0;i<N1;i++)
  {P[i]=-G[i]+beta*P[i];
  		
  }
goto step3; 
	   }
    else goto step7;
 
//=========   step 1  设初始值    ===========/
step3: NLSMAX=50;
	   FO=F;
      for( I=0;I<N1;I++)
		  WORK[I]=X[I];
       DGO=DDOT1(D,G);
       A1=0.0;
       F1=F;
       DG1=DGO;
       A2=-1.0;
       A21=0;
       RATIO1=1.0;
       RATIO2=0.1;
       NCALFN=0,NCALG=1,NCALF=1;
       ALPHA=1.0;
	   //INFO;
       if(DGO>=0.0) goto enter_1100;
       if(IPRINT!=0) goto enter_40;
	   cout<<"ENTER LINE SEARCH"<<endl;
////////////////////////////////////////////////////
//==========   step 2  计算新的试验点  ===========//
      enter_40 : for(I=0;I<N1;I++)
					 X[I]=WORK[I]+ALPHA*D[I];
      CALF(N1,X,F);//目标函数值
      NCALFN=NCALFN+1;
      if(NCALFN>NLSMAX) goto enter_1300;
      if(F<=FO+C1*ALPHA*DGO) goto enter_200;

//==========  Step 3   如果第一个等式不成立========// 
//==========	  通过二次修正降低步长   ==========//
      RATIO1=RATIO1*0.1;
      if(DGO/(1.0+fabs(FO))<-1.0e-16)  goto enter_80;
      DD=DDOT1(D,D);
      XX=DDOT1(X,X);
      if(ALPHA*sqrt(DD)/(1.0+sqrt(XX))<=1.0e-12) goto enter_1400;
      enter_80 : A2=ALPHA;
      F2=F;
      A21=0;
      enter_90 : TRY=((F2-F1)/(A2-A1)-DG1)/(A2-A1);
      TRY=-DG1/TRY/2.0;
      goto enter_240;

//=========   Step 4 如果第一个满足测试第二个不等式=========//
      enter_200 : CALG(N1,X,G);//函数的梯度
      NCALG=NCALG+1;
      DG=DDOT1(D,G);
      if(fabs(DG)<=-C2*DGO) goto enter_1500;
      if(DG>0.0) goto enter_210;
      if(A2<0.0) goto enter_260;
      RATIO2=RATIO2*0.1;
      A1=ALPHA;
      F1=F;
      DG1=DG;
      RATIO1=0.1;
      if(!A21) goto enter_90;
      goto enter_220;

//==========   Step 5 使用三次插值寻找===========// 
//                 在[A1,A2]内的新的点           //
      enter_210 : RATIO1=0.1;
      A2=ALPHA;
      F2=F;
      DG2=DG;
      A21=1;
      enter_220 : 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 enter_230;
      if(sqrt(TRY)<=T3) goto enter_230;
      TRY=-DG1*(A2-A1)/(sqrt(TRY)+T3);
      goto enter_240;
      enter_230 : TRY=((F2-F1)/(A2-A1)-DG1)/(A2-A1);
      TRY=-DG1/TRY/2.0;

//========   如果必要去掉这步==========//
      enter_240 : if(TRY<RATIO1*(A2-A1)) TRY=RATIO1*(A2-A1);
      if(TRY>(1.0-RATIO2)*(A2-A1)) TRY=(1.0-RATIO2)*(A2-A1);
              ALPHA=TRY+A1;
      goto enter_40;

//=========   Step 6   ALPHA是最大的可行点,通过使用三次插值找到另一个可行的值使STEP>ALPHA=======//
      enter_260 : 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 enter_275;
      if(sqrt(TRY)<=T3) goto enter_275;
      TRY=-DG*(ALPHA-A1)/(sqrt(TRY)-T3);
      if(TRY<=1.0*10*(ALPHA-A1)) goto enter_280;
      enter_275 : TRY=1.0*10*(ALPHA-A1);
      enter_280 : if(A2<0.0) goto enter_290;
      if(TRY<0.1*(A2-A1)) TRY=0.1*(A2-A1);
      if(TRY>0.9*(A2-A1)) TRY=0.9*(A2-A1);
      enter_290 : A1=ALPHA;
      F1=F;
      DG1=DG;
      ALPHA=A1+TRY;
	  goto step6;
      goto enter_40;


//==========   使用上升方向搜索或者参量是不正确的   ===========//
      enter_1100 : //cout<<"AN UPHILL SEARCH DIRECTION IS USED!"<<endl;
      INFO=1;
      goto enter_1500;
      enter_1300 : //cout<<"LINE SEARCH FAILS AFTER NLSMAX TRIAL STEPSIZES!"<<endl;
      INFO=2;
      goto enter_1500;
      enter_1400 : CALG(N1,X,G);
      NCALG=NCALG+1;
      //cout<<"STEP IS TOO SMALL!"<<endl;
      INFO=3;
      enter_1500 : NCALF=NCALF+NCALFN;
   //return;
step6:
 k=k+1;
 goto step1;
      
step7: CALF(N1,X,F);
	   cout<<"问题的规模:"<<N1<<endl<<"迭代次数:"<<k<<endl<<"精度为:"<<kesi<<endl<<"F="<<F<<endl;
	   end=clock();
	   cout<<"计算用时:"<<(double)(end-start)/CLOCKS_PER_SEC<<"ms"<<endl;
}
💿 文件大小 4 K
👤 上传用户 lidaniel
📂 所属分类数据结构
🏷️ 相关标签

#运行程序 #家
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -