c++gongetidufa.c

可以通过目标行数进行

#include "stdio.h"
#include "math.h"


//***目标函数***//
double func1(double x[])
{
    return(5*x[0]*x[0]+4*x[1]*x[1]+2*x[2]*x[2]+11*x[3]*x[3]+x[4]*x[4]);
}
double func(double t,double x[],double p[])
{
    double xx[5];
    int i=0;
    for(i=0;i<5;i++)
        xx[i]=x[i]+t*p[i];
    return(func1(xx));
}
//***目标函数***//


//***目标函数的梯度***//
double g1(double x[])
{
    return(10*x[0]);
}
double g2(double x[])
{
    return(8*x[1]);
}
double g3(double x[])
{
    return(4*x[2]);
}
double g4(double x[])
{
    return(22*x[3]);
}
double g(double x[])
{
    return(2*x[4]);
}
//***目标函数的梯度***//


//***搜索区间确定***//
void sosuoqujian(double *a,double *b,double x[],double p[])
{
    double t0=0.0;
    double h=1.0;
    double s,t1,t2,t3,v,u,w,a,b,c;
    a=func(t0,x,p);
    b=func(t0+h,x,p);
    if(b>a)
    {   
        s=-h;
        c=func(t0+s,x,p);
        if(c>a)
        {
            t1=t0+s;
            t2=t0;
            t3=t0-s;
            *a=t1;
            *b=t3;
            return;
        }
        else
            b=c;
    }
    else 
    { 
        s=h;
    }
    u=t0;
    v=t0+s;
    while(1)
    {  
        s=2.0*s;
        w=u+s; 
        c=func(w,x,p);  
        if(c<=b)  
        {     
                
            a=b;
            u=w; 
            b=c;          
        }  
        else 
            break; 
    }
    if(u>w)
    {
        *a=w; 
        *b=u;
    }
    else 
    {
        *a=u;
        *b=w; 
    }
}  
//***搜索区间确定***//


//***黄金分割法直线搜索***// 
double goldencut(double a,double b,double x[],double p[])
{
    double beta,fa1,fa2,t1,t2,t;
    double epxl=1.0e-5;
    beta=(sqrt(5.0)-1.0)/2.0; 
    t2=a+beta*(b-a); 
    fa2=func(t2,x,p);
    while(1)
    {
        t1=a+b-t2; 
        fa1=func(t1,x,p);
        while(1)
        {
            if(fabs(t1-t2)<epxl)
        	{
                t=(t1+t2)/2.0;
                return(t);
        	}
            else 
        	{
                if(fa1<=fa2)
                {  
                    b=t2;
                    t2=t1;
                    fa2=fa1;
                    break;
            	} 
                else 
            	{
                    a=t1;     
                    t1=t2;
                    fa1=fa2;  	
                    t2=a+beta*(b-a);  
                    fa2=func(t2,x,p);
            	}
        	}
        }
    }
}
//***黄金分割法直线搜索***// 


//***H终止准则***//
int hlimmelblau(double x1[],double x2[],double fk,double fl,double g[])
{
    double epxl1=1.0e-5;
    double epxl2=1.0e-5;
    double epxl3=1.0e-4;
    double row;
    if(sqrt(g[0]*g[0]+g[1]*g[1]+g[2]*g[2]+g[3]*g[3]+g[4]*g[4])>=epxl3)
        return(1);
    else 
    {
        row=sqrt(x1[0]*x1[0]+x1[1]*x1[1]+x1[2]*x1[2]+x1[3]*x1[3]+x1[4]*x1[4]);
        if(row<epxl2)
            row=1.0;
        else 
            row=row*epxl1;
        if(((x2[0]-x1[0])*(x2[0]-x1[0])+(x2[1]-x1[1])*(x2[1]-x1[1])+(x2[2]-x1[2])*(x2[2]-x1[2])+(x2[3]-x1[3])*(x2[3]-x1[3])+(x2[4]-x1[4])*(x2[4]-x1[4]))>=row)
            return(1);
        else
            row=fabs(fk);
        if(row<epxl2)
            row=1.0;
        else
            row=epxl1*row;
        if(fabs(fl-fk)>=row)
            return(1);
        else
            return(0);      	
    }
}
//***H终止准则***//


//***共轭梯度法函数***//
void main()
{
    double t,a,b,arfa,x0[3],f0,f,g0[5],p0[5],p[5],g[5],x[5],epxl;
    double *point1,*point2;
    point1=&a;
    point2=&b;
    int k,i,j;
    j=k=0;
    epxl=1.0e-5;
    cout<<"请输入初始点（五维）：\n";
    for(i=0;i<5;i++)
         cin>>x0[i];    	 
    for(i=0;i<5;i++)
     x[i]=x0[i];
    f=func1(x0);
    g[0]=g1(x0);
    g[1]=g2(x0);
    g[2]=g3(x0);
        g[3]=g4(x0);
        g[4]=g5(x0);
    while(1)
    {
        for(i=0;i<=4;i++)  
            p0[i]=-g[i];
        while(1)    
        {	
            for (i=0;i<5;i++)   	
        	{
                x0[i]=x[i];
                g0[i]=g[i];
        	}
        
            f0=f;   
            sosuoqujian(point1,point2,x0,p0);   
            t=goldencut(a,b,x0,p0);
            for(i=0;i<=4;i++)      	
                x[i]=x0[i]+t*p0[i]; 
            printf("%f %f %f \n",x[0],x[1],x[2],x[3],x[4]);
            f=func1(x); 
            printf("搜索点x:\n%f %f %f %f %f\n 函数值：f= %f\n",x[0],x[1],x[2],x[3],x[4],f);
            j++;
            g[0]=g1(x);
            g[1]=g2(x);
            g[2]=g3(x);
                    g[3]=g4(x);
                        g[4]=g5(x);
            if(hlimmelblau(x0,x,f0,f,g))
            {		
                if(k==3)    	
                {  		
                    k=0;      			
                    break;    		
                }		
                else	
                {  			
                    arlfa=(g[0]*g[0]+g[1]*g[1]+g[2]*g[2]+g[3]*g[3]+g[4]*g[4])/(g0[0]*g0[0]+g0[1]*g0[1]+g0[2]*g0[2]+g0[3]*g0[3]+g0[4]*g0[4]);    	
                    for(i=0;i<5;i++)           
                        p[i]=-g[i]+arlfa*p0[i];     	
                    if(fabs(p[0]*g[0]+p[1]*g[1]+p[2]*g[2]+p[3]*g[3]+p[4]*g[4])<=epxl)   		
                    {				
                        k=0;				
                        break;				
                    }			
                    else			
                    {					
                        if((p[0]*g[0]+p[1]*g[1]+p[2]*g[2]+p[3]*g[3]+p[4]*g[4])<-epxl)       		
                        {					
                            for(i=0;i<5;i++)        				
                                p0[i]=p[i];					
                            k=k+1;					
                        }     					
                        else      					
                        { 				
                            for(i=0;i<5;i++)        					
                                p0[i]=p[i];					
                            k=k+1;						
                        }   						
                    }	  			
                }			
            }			
            else    	
                goto loop1;     
    
        }
    }
loop1:   cout<<"最优点："<<x[0]<<"    "<<x[1]<<"  "<<x[2]<<" "<<x[3]<<"  "<<x[4]<<endl<<"最优值："<<f<<endl<<"计算步数："<<j<<endl;


}