zyh-1.c

最优化理论与方法

                // #endif
                
                // #ifdef lsno
               //t=lsno_accurate(a_xk,a_gk,a_pk,n,a_v,a_lumda,a_miu,n1,n2,a_xk_last);
			    //printf("The lsno_accurate's t=%f(FR)\n",t);
				//getch();
                //#endif
        		
                                   
                
                    for(i=0;i<n;i++)
                      {
                       a_xk[i]=a_xk[i]+t*a_pk[i];   /*Xk+1*/
                       printf("the (Xk+1): xk%d=%f gk%d=%f pk%d=%f\n",i,a_xk[i],i,a_gk[i],i,a_pk[i]);
                      }            
            	                     printf("\n");
                        for(i=0;i<n;i++)
                           a_gk[i]=d_fun(a_xk,i,n,a_v,a_lumda,a_miu,n1,n2);  /*gk+1*/
                     
                        for(i=0;i<n;i++)
                           printf("xk%d=%f gk%d=%f ",i,a_xk[i],i,a_gk[i]);
                                    //getch();
                 
                         result=Him(a_xk,a_xk_last,a_gk,a_gk_last,n,a_v,a_lumda,a_miu,n1,n2);
                    
                    if(result==1)
                        {
                         result1=function(a_xk,control_df,n,a_v,a_lumda,a_miu,n1,n2);
                         printf("result=%f",result1);
						 printf("\n************          The FR end              *************\n");
						 //getch();
                         return (result1);

                        };
                  printf("time_k=%d\n",++time_k);
                    if(time_k==m)
                        {
                         time_k=0;
                         for(i=0;i<n;i++)
							{a_pk[i]=-a_gk[i];printf("pk%d=%f ",i,a_pk[i]);}printf("\n");
						    // getch();
						}
                      else
					  {         
						        
                          t2=fanshu(a_gk,n)*fanshu(a_gk,n);
                               
                          t1=fanshu(a_gk_last,n)*fanshu(a_gk_last,n);
                          alpha=t2/t1;
                          printf("alpha=%f t2=%f %f t1=%f %f\n",alpha,t2,pow(fanshu(a_gk,n),2),t1,pow(fanshu(a_gk_last,n),2));
                                //getch();
                          for(i=0;i<n;i++)
			    printf("pk%d=%f ",i,a_pk[i]);
			    printf("\n");


                          for(i=0;i<n;i++)
                            a_pk[i]=-a_gk[i]+alpha*a_pk_last[i];
                               
                          for(i=0;i<n;i++)
			    printf("pk%d=%f ",i,a_pk[i]);
			    printf("\n");

                          sum1=0;

                          for(i=0;i<n;i++)
                            sum1=sum1+a_pk[i]*a_gk[i];
                            printf("sum1=%f\n",sum1);
                                //getch();
                          if(fabs(sum1)>0.0001)
                            {
                              if(sum1>-0.0001)
                                  {
                                    for(i=0;i<n;i++)
                                     a_pk[i]=-a_pk[i];
                                   }
                              ++time_k;
                             }
                           else
                            {
                               k=0;
                               for(i=0;i<n;i++)
                               a_pk[i]=-a_gk[i];
                             }
                                
                        }

            keeplastp(a_xk,a_xk_last,a_gk,a_gk_last,a_pk,a_pk_last,n);

          }
	  for(i=0;i<n;i++)
		  printf("xk%d=%f ",i,a_xk[i]);
          printf("the FR's M=%d\n",M);
}

double BFGS(double a_xk[],double a_xk_last[],double a_gk[],double a_gk_last[],double a_pk[],double a_pk_last[],double a_H[][N],double a_H_last[][N],int n,double a_v[],double a_lumda[],double a_miu,int n1,int n2)
{
 int i,j,a,m=15,M=0,k=0,result1;
 double yk[N],sk[N],t,result,a1,w[N],H[N][N],H1[N][N],H2[N][N],H3[N][N],sum;
  printf("************      The BFGS       *************\n");
                                //getch();

  for(i=0;i<n;i++)
	{ 
	  printf("xk%d=%f gk%d=%f pk%d=%f\n",i,a_xk[i],i,a_gk[i],i,a_pk[i]);
          printf("xk_%d=%f gk_%d=%f pk_%d=%f\n",i,a_xk_last[i],i,a_gk_last[i],i,a_pk_last[i]);
	} 

while(++M)
   {
        //#ifdef ls0 
         t=ls0618(a_xk,a_pk,n,a_v,a_lumda,a_miu,n1,n2,a_xk_last);printf("the ls0618's t is %f(BFGS)\n",t);
        //#endif
        

         //#ifdef lsno
        //t=lsno_accurate(a_xk,a_gk,a_pk,n,a_v,a_lumda,a_miu,n1,n2,a_xk_last);printf("the lsno_accurate's t(BFGS)  is %f\n",t);
		//getch();
         //#endif
          for(i=0;i<n;i++)
			  printf("xk%d=%f xk_%d=%f gk%d=%f gk_%d=%f\n",i,a_xk[i],i,a_xk_last[i],i,a_gk[i],i,a_gk_last[i]);

		  for(i=0;i<n;i++)
             a_xk[i]=a_xk[i]+t*a_pk[i];/*xk+1*/
          for(i=0;i<n;i++) 
			  a_gk[i]=d_fun(a_xk,i,n,a_v,a_lumda,a_miu,n1,n2);/*gk+1*/
		  
		  for(i=0;i<n;i++)
			  printf("xk%d=%f xk_%d=%f gk%d=%f gk_%d=%f\n",i,a_xk[i],i,a_xk_last[i],i,a_gk[i],i,a_gk_last[i]);
		  
		  result1=Him(a_xk,a_xk_last,a_gk,a_gk_last,n,a_v,a_lumda,a_miu,n1,n2);
          
                            

           if(result1==1)
             {
               result=function(a_xk,control_df,n,a_v,a_lumda,a_miu,n1,n2);
               printf("************      The BFGS End!!       *************\n");
			   return (result);
               break;
             }
           if(k==m)
             {
               k=0;
               for(i=0;i<n;i++)
                a_pk[i]=-a_gk[i];
		for(i=0;i<n;i++)
                   for(j=0;j<n;j++)
                       {
                        if(i==j)a_H[i][j]=1;
                        if(i!=j)a_H[i][j]=0;  
 		       }

	   //printf("the k is equal to m and start from next_point\n");
			   //getch();
             }
             else
			 { 
				  for(i=0;i<n;i++)
					  printf("xk%d=%f xk_%d=%f gk%d=%f gk_%d=%f pk%d=%f pk_%d=%f\n",i,a_xk[i],i,a_xk_last[i],i,a_gk[i],i,a_gk_last[i],i,a_pk[i],i,a_pk_last[i]);
                 
				 for(i=0;i<n;i++)
                   {
                    yk[i]=a_gk[i]-a_gk_last[i];
                    sk[i]=a_xk[i]-a_xk_last[i];
                   }

                            for(i=0;i<n;i++)
								{
								 printf("yk%d=%f sk%d=%f \n",i,yk[i],i,sk[i]);
								}
				            //getch();
				 
				 a1=0;
                  for(i=0;i<n;i++)
                    a1=a1+sk[i]*yk[i];    /*sk'*yk*/
				            
				            printf("sk'*yk=%f\n",a1);
					        //getch();

                  for(i=0;i<n;i++)
                    {
					  for(j=0;j<n;j++)
				    	H[i][j]=sk[i]*sk[j];
					}  /*sk*sk'*/
				            
				             printf("sk*sk'\n");
				             for(i=0;i<n;i++)
							   for(j=0;j<n;j++)
									{printf("H[%d][%d]=%f ",i,j,H[i][j]);
									 if(j==n-1)printf("\n");}

                  for(i=0;i<n;i++)
                    {for(j=0;j<n;j++)
                      {
                        H1[i][j]=yk[i]*sk[j]; /*yk*sk',sk*yk'*/
                        H2[i][j]=sk[i]*yk[j];
                      }
                     }

                             
				              printf("yk*sk'=H1[]\n");
				              for(i=0;i<n;i++)
								for(j=0;j<n;j++)
									{
									  printf("H1[%d][%d]=%f",i,j,H1[i][j]);
									  if(j==n-1)printf("\n");
									}
							    printf("sk*yk'=H2[]\n");  
                                                               for(i=0;i<n;i++)
								  for(j=0;j<n;j++)
								  {
								   printf( "H2[%d][%d]=%f ",i,j,H2[i][j]);
								   if(j==n-1)printf("\n");
								  }
				  
				  
				  for(i=0;i<n;i++)
                     {
                       w[i]=0;
                       for(j=0;j<n;j++)
                         w[i]=w[i]+yk[j]*a_H[j][i];   /*yk'*Hk*/
                      }

				                     printf("yk'*Hk\n");
                                     for(i=0;i<n;i++)
										 printf("w%d=%f ",i,w[i]);
                                       printf("\n");
                   sum=0;
                   for(i=0;i<n;i++)
                     sum=sum+w[i]*yk[i];   /*1+(yk'*Hk*yk)/(sk'*yk)*/
                     sum=sum/a1+1;

                                       printf("1+(yk'*Hk*yk)/(sk'*yk)=%f\n",sum); 
					 
				   for(i=0;i<n;i++)
                      for(j=0;j<n;j++)
                         H[i][j]=sum*H[i][j]; /*(1+(yk'*Hk*yk)/(sk'*yk))*sk*sk'*/

					              printf("(1+(yk'*Hk*yk)/(sk'*yk))*sk*sk'\n");
								  for(i=0;i<n;i++)
									for(j=0;j<n;j++)
										{printf("H[%d][%d]=%f ",i,j,H[i][j]);
									    if(j==n-1)printf("\n");}

                   sum=0;
                 for(i=0;i<n;i++)
				   {  
					   for(j=0;j<n;j++)
					   {
						   sum=0;
						   for(a=0;a<n;a++)
						   sum=sum+a_H[i][a]*H1[a][j];  /*Hk*yk*sk'*/
				           H3[i][j]=sum;
					   }
				   }     /*Hk*yk*sk'*/  

					              printf("Hk*yk*sk'\n");
								  for(i=0;i<n;i++)
									for(j=0;j<n;j++)
										{
									     printf("H3[i][j]=%f ",i,j,H3[i][j]);
										 if(j==n-1)printf("\n");
										}

                   sum=0;
                  for(i=0;i<n;i++)   /*sk*yk'*Hk*/
				   {
					   for(j=0;j<n;j++)
						{
					       sum=0;
						   for(a=0;a<n;a++)
						    sum=sum+H2[i][a]*a_H[a][j];
					        H1[i][j]=sum;	
					   } 
				   } /*sk*yk'*Hk*/

					              printf("sk*yk'*Hk\n");
								  for(i=0;i<n;i++)
									for(j=0;j<n;j++)
										{
									     printf("H1[%d][%d]=%f ",i,j,H1[i][j]);
										 if(j==n-1)printf("\n");
										}

                   for(i=0;i<n;i++)
                     for(j=0;j<n;j++)
                       {
                        a_H[i][j]=a_H[i][j]+(1/a1)*(H[i][j]-H3[i][j]-H1[i][j]);/*Hk+1*/
                        a_H_last[i][j]=a_H[i][j];
						printf("H[%d][%d]=%f ",i,j,a_H[i][j]);
						if(j==n-1)printf("\n");
                       }
					            

                   for(i=0;i<n;i++)
                     {
                      sum=0;
                       for(j=0;j<n;j++)
                         sum=sum+a_H[i][j]*a_gk[j];
                         a_pk[i]=-sum;
                     }
                ++k;
              for(i=0;i<n;i++)
				  printf("xk%d=%f pk%d=%f",i,a_xk[i],i,a_pk[i]);
			   printf("\n");

           }
       keeplastp(a_xk,a_xk_last,a_gk,a_gk_last,a_pk,a_pk_last,n);
	   
	   
	   printf("the time(BFGS) is %d\n",k);
               //getch();
    }

}


main()
{
  int n=N,i,j,n1=N1,n2=N2,k=0;
  double pk[N],pk_last[N],xk[N]={2,4},xk_last[N],gk[N],gk_last[N],v[N1]={0},lumda[N2]={0},H[N][N],H_last[N][N];
  double result,c=10,m=0,faik=0,faik_last=0,t=0,miu=1;
      printf("*****************************************\n");
       for(i=0;i<n;i++)
         {
           gk[i]=d_fun(xk,i,n,v,lumda,miu,n1,n2);
           gk_last[i]=gk[i];                      /*x点，梯度，P方向初始化*/
           pk[i]=-gk[i];
           xk_last[i]=xk[i];
           pk_last[i]=pk[i];
         }
        for(i=0;i<n;i++)
         {
         printf("xk%d=%f xk_%d=%f gk%d=%f gk_%d=%f pk%d=%f pk_%d=%f\n",i,xk[i],i,xk_last[i],i,gk[i],i,gk_last[i],i,pk[i],i,pk_last[i]);
         }
        for(i=0;i<n1;i++)
			printf("v%d=%f ",i,v[i]);
		    printf("\n");
		for(i=0;i<n2;i++)
			printf("lumda%d=%f",i,lumda[i]);
		    printf("\n");
		
		 printf("the miu=%f\n",miu);
        for(i=0;i<n;i++)
         for(j=0;j<n;j++)
           {
            if(i==j){H[i][j]=1;H_last[i][j]=1;}
              else {H[i][j]=0;H_last[i][j]=0;}
           }
         for(i=0;i<n;i++)
           for(j=0;j<n;j++)
              {
               printf("%f ",H[i][j]);
               if(j==(n-1))printf("\n");
              }
         printf("\n*****************************************\n");
          faik=fai_k(xk,v,n1,n2,miu);
		  faik_last=faik;
		  printf("faik=%f\n",faik);
		  //getch();
		while(k<=500)
			{
				//#ifdef F
			   FR(xk,xk_last,gk,gk_last,pk,pk_last,n,v,lumda,miu,n1,n2);
				
				
		    	//	#else
                 //BFGS(xk,xk_last,gk,gk_last,pk,pk_last,H,H_last,n,v,lumda,miu,n1,n2);
				
			
				//#endif
				 faik=fai_k(xk,v,n1,n2,miu);  /*fai(k+1)*/
				 for(i=0;i<n;i++)
					 printf("xk%d=%f ",i,xk[i]);
				 printf("v0=%f lumda0=%f miu=%f faik=%f\n",v[0],lumda[0],miu,faik);
				    //getch();
				 if(faik<0.1)
					{
					 printf("the xk are:");
					 for(i=0;i<n;i++)
						 printf("xk%d=%f ",i,xk[i]);
						 printf("\n");
						// getch();
					   break;
					}
			    
				 if(faik/faik_last>0.5)
					{
				     miu=c*miu;
					 //printf("faik/faik_last>0.5 and the miu=%f\n",miu);
					 //getch();
					}
				 for(i=0;i<n1;i++)
					{
				      m=v[i]-2*miu*function_s(xk,i);
					  printf("v[%d]-2*%f*function_s(xk,%d)=%f\n",i,miu,i,m);
					 // getch();
					  v[i]=max(t,m);
					}
				 for(i=0;i<n2;i++)
				 {
				   lumda[i]=lumda[i]-2*miu*function_h(xk,i);
				   printf("lumda%d=%f \n",i,lumda[i]);
				 }
		      ++k;
			  for(i=0;i<n1;i++)
				  printf("v%d=%f ",i,v[i]);
			      printf("\n");
			  for(i=0;i<n2;i++)
				  printf("lumda%d=%f ",i,lumda[i]);
			      printf("\n");
			  printf("miu=%f faik_last=%f faik=%f\n",miu,faik_last,faik);
			  for(i=0;i<n;i++)
				{
				  printf("xk%d=%f gk%d=%f pk%d=%f\n",i,xk[i],i,gk[i],i,pk[i]);
			      printf("xk_%d=%f gk_%d=%f pk%d=%f \n",i,xk[i],i,gk[i],i,pk[i]);
				}
			  faik_last=faik;
              result=function(xk,control_df,n,v,lumda,miu,n1,n2);
			  printf("the result=%f\n",result);
			  getch();
			}
     
	 result=function(xk,control_df,n,v,lumda,miu,n1,n2);
	 for(i=0;i<n;i++)
		 printf("xk%d=%f ",i,xk[i]);
	     printf("\n");
	 printf("the langer times is %d\n",k);
	 printf("the function value is %f\n",result);
		   getch();
}