powell.cpp

来自「优化设计中powell法的类」· C++ 代码 · 共 501 行

#include <iostream.h>
#include <math.h>

double fxy1(double *p)
{
	return (10*(*p+*(p+1)-5)*(*p+*(p+1)-5)+(*p-*(p+1))*(*p-*(p+1)));
}


double fxy2(double *p)
{
	return ((*p-2)*(*p-2)+(*(p+1)-3)*(*(p+1)-3)+(*(p+2)-10)*(*(p+2)-10));
}

double fxy(double *p)
{ 
	return ((*p-2)*(*p-2)+(*(p+1)-3)*(*(p+1)-3)+1.0/(*p)+1.0/(*(p+1))); 
}
double fxy4(double *p)
{
	return(*p);
}
class powell
{
private:
	int n;
	double jing;
	double h;

	double kkk[6];



public:
	void tansuo(double *b,double *e1,int ii)
	{
	   int j; 

       double h0;
	   double yy2;
	   double x1,x2,x3,y1,y2,y3,xx,yy;
	 
	   h0 = h;
	 
	 
	   double *pm = new double [n];  //临时存放前一分量
	   double *pt = new double [n];  //临时存放探索分量
	   double *ee = new double [n];  //临时存放方向向量

	   double *p1 = new double [n];  //临时存放区间分量
	   double *p2 = new double [n];  //临时存放区间分量
	   double *p3 = new double [n];  //临时存放区间分量

	   
       double aa,bb,cc;  //用于探索区间
		
	   
	   for(j=0;j<n;j++)
		{
            *(pm+j) = *(b+j*(n+1)+ii-1);  //把初值传给pm  
		}
		aa = fxy(pm);   //把pm的函数值给aa

	
		for(j=0;j<n;j++)
		{
			*(ee+j) = *(e1+j*n+ii-1);  //设置ee为当前方向
		}


		

		for(j=0;j<n;j++)
		{
			*(pt+j) = *(pm+j)+*(ee+j)*h;   //向前一个步长
		}
        bb = fxy(pt);

			
		if(aa>bb)   
		{   
			for(j=0;j<n;j++)
				*(p2+j) = *(pt+j);  //pt当前值给p2
			
			cc = bb-10;  
			while(bb>cc)
			{
				h = h*2;
				for(j=0;j<n;j++)
					*(pt+j) = *(pt+j)+*(ee+j)*h;

				cc = fxy(pt);
			}

			for(j=0;j<n;j++)
			{
				*(p1+j) = *(pm+j);
				*(p3+j) = *(pt+j);
			}
			
			x1 = 0;
			x2 = h0;
			x3 = h0 + h;
			
		}
		else
		{
			cc = aa-10;
			h = -h/16;
			for(j=0;j<n;j++)
			{
				*(p3+j) = *(pt+j);
			}

			for(j=0;j<n;j++)
			{
				*(pt+j) = *(pm+j);
			}
				
			while(cc<aa)
			{
				h = 2*h;
				for(j=0;j<n;j++)
				{
					*(pt+j) = *(pt+j)+*(ee+j)*h;
				}
					
				cc = fxy(pt);
			}

		
			for(j=0;j<n;j++)
			{
				*(p1+j) = *(pt+j);
				*(p2+j) = *(pm+j);
			}

			x1 = h;
			x2 = 0;
			x3 = h0;
		}
       
		

		aa = fxy(p1);
		bb = fxy(p2);
		cc = fxy(p3);

		
            y1 = aa;
			y2 = bb;
			y3 = cc;
		
			double *zx = new double [n];
						
			yy = jing*10+y2+10;
			double xx2;
			xx = 0;
			while((yy-yy2)*(yy-yy2)>jing*100)
			{
				xx2 = xx;
				
				xx = 0.5*(x1+x3-((y3-y1)/(x3-x1)*(x2-x3))/((y2-y1)/(x2-x1)-(y3-y1)/(x3-x1)));  //二次插值公式
				yy2 = y2;
						
				for(j=0;j<n;j++)
				
					*(zx+j) = *(pm+j) + *(ee+j)*xx;
				                 				
               	yy = fxy(zx);

				if(xx<x2)
				{
					x3 = x2;
					x2 = xx;
					y3 = y2;
					y2 = yy;
                	
				}
				else
				{   
				
					x1 = x2;
				    x2 = xx;				
					y1 = y2;
					y2 = yy;
				}

				if(yy>1e+010||yy<-1e+010||yy2>1e+010||yy2<-1e+010)
				{
					xx = xx2;
					break;
				}

				
			}
		

		  delete zx;

		for(j=0;j<n;j++)	
		  *(b+j*(n+1)+ii) = *(pm+j) + *(ee+j)*xx;  //一维探索得到一个新点

		

		delete pt;
		
		delete pm;
		delete ee;

		delete p1;
		delete p2;
		delete p3;	

		

	}

	void action(int nn,double jingdu)
	{
		int i,j;
		double F0,F2,F3;
		double maxfc,aa,bb,cc,h0,yy2;
		int m;
		double x1,x2,x3,y1,y2,y3,xx,yy;


		n = nn;
		jing = jingdu;


		h = 10;  //设定初始步长
		h0 = h;

		double *e = new double [n*n];  //  方向向量
		double *p = new double [n*(n+1)];  // 存储全部点的分量
        
		
		for(i=0;i<n*(n+1);i++)  
			*(p+i) = i+1;

	
       
		for(i=0;i<n*n;i++)  //设置方向向量
		{
			if((i/n)==(i%n))
			       *(e+i) = 1;
			else 
				*(e+i) = 0;
		}
   
		double ass;
		ass = jing*jing+10;

   while(sqrt(ass)>jing)
   {  
	  for(i=1;i<n+1;i++)  //进行一维探索
	  {
		  h = 10;
		  tansuo(p,e,i);
	  }

	  

	  double *px = new double [n];
	  double *p0 = new double [n];
	  double *pn = new double [n];
	  double *pk = new double [n];  //方向

	  
	  

	  for(i=0;i<n;i++)         //设置始点，终点，反射点
	  {
          *(p0+i) = *(p+i*(n+1));
		  *(pn+i) = *(p+i*(n+1)+n);
          *(px+i) = 2**(p+i*(n+1)+n) - *(p+i*(n+1));
	  }	  *(pk+i) = *(p+i*(n+1)+n)-*(p+i*(n+1));

       F0 = fxy(p0);
	   F2 = fxy(pn);
	   F3 = fxy(px);

	   double *fx = new double [n+1];   // 记录各中间点的函数值
	   for(i=0;i<n+1;i++)
	   {
		   double *temp = new double [n];  //临时存储分量
		   for(j=0;j<n;j++)
		   {
			   *(temp+j) = *(p+j*(n+1)+i);
		   }
		   *(fx+i) = fxy(temp);
		   delete temp;
	   }
	   

	   double *fc = new double [n];  //存储n个函数值之差
	   for(i=0;i<n;i++)
	   {
		   *(fc+i) = *(fx+i) - *(fx+i+1);
	   }

	   maxfc = *fc;   //最大函数值差
	   for(i=0;i<n;i++)
	   {
		   if(*(fc+i)>maxfc)
		   {
			   maxfc = *(fc+i);
			   m = i;
		   }
	   }
	   delete fc;
	 
	   delete fx;

	   int some;
	   some = (F3>F0)&&(((F0-2*F2+F3)*(F0-F2-maxfc)*(F0-F2-maxfc))<(0.5*maxfc*(F0-F3)*(F0-F3)));
	   
	   if(0)
	   {
           h = 10;
		   for(i=m;i<n-1;i++)
			   for(j=0;j<n;j++)
				   *(e+j*n+i) = *(e+j*n+i+1);

		   for(j=0;j<n;j++)
		        *(e+j*n+n-1) = *(pk+j);

		   double *k1 = new double [n];
		   double *k2 = new double [n];
		   double *k3 = new double [n];
		   double *kt = new double [n];   //临时存储探索分量

		   aa = fxy(p0);

		   for(i=0;i<n;i++)
			   *(kt+i) = *(p0+i) + *(pk+i)*h;

		   bb = fxy(kt);

		   if(aa>bb)
		   {
			   for(i=0;i<n;i++)
				   *(k2+i) = *(kt+i);
			   cc = bb-10;

			   while(bb>cc)
			   {
				   h = h*2;
				   for(i=0;i<n;i++)
					   *(kt+i) = *(kt+i) + *(pk+i)*h;
				   cc = fxy(kt);
			   }

			   for(i=0;i<n;i++)
			   {
				   *(k1+i) = *(p0+i);
				   *(k3+i) = *(kt+i);
			   }
			 
			   x1 = 0;
			   x2 = h0;
			   x3 = h0 + h;
		   }

		   else
		   {
			   cc = aa-10;
			   h = -h/16;

			   for(i=0;i<n;i++)
				   *(k3+i) = *(kt+i);
			   for(i=0;i<n;i++)
				   *(kt+i) = *(p0+i);

			   while(cc<aa)
			   {
				   h = h*2;
				   for(i=0;i<n;i++)
					   *(kt+i) = *(kt+i) + *(pk+i)*h;
				   cc = fxy(kt);
			   }
			   for(i=0;i<n;i++)
			   {
				 *(k1+i) = *(kt+i);
				 *(k2+i) = *(p0+i);
			   }
			 
			   x1 = h;
			   x2 = 0;
			   x3 = h0;
		   }

		   aa = fxy(k1);
		   bb = fxy(k2);
		   cc = fxy(k3);

		   
            y1 = aa;
			y2 = bb;
			y3 = cc;
	
			double *zx = new double [n];
              			
			yy = jing*10+y2+10;
			double xx2;
			xx = 0;
			while((yy-yy2)*(yy-yy2)>jing*jing)
			{
				xx2 = xx;
				
				xx = 0.5*(x1+x3-((y3-y1)/(x3-x1)*(x2-x3))/((y2-y1)/(x2-x1)-(y3-y1)/(x3-x1)));  //二次插值公式
                yy2 = y2;

			for(j=0;j<n;j++)
				*(zx+j) = *(p0+j) + *(pk+j)*xx;
				
					yy = fxy(zx);

				if(xx<x2)
				{
					x3 = x2;
					x2 = xx;
					y3 = y2;
					y2 = yy;
				}
				else
				{					
					x1 = x2;
					x2 = xx;
					y1 = y2;
					y2 = yy;
				}
				if(yy>1e+010||yy<-1e+010||yy2>1e+010||yy2<-1e+010)
				{
					xx = xx2;
					break;
				}
									
			}
	

		  delete zx;

		   for(i=0;i<n;i++)
			   *(kt+i) = *(p0+i) + *(pk+i)*xx;   // kt是下一轮始点
		   for(i=0;i<n;i++)
			   *(p+i*(n+1)) = *(kt+i);    //输入初始位置

		   delete k1;
		   delete k2;
		   delete k3;
		   delete kt;
		   
	   }
	   else
	   {
		   if(F2>F3)
			   for(i=0;i<n;i++)
				   *(p+i*(n+1)) = *(px+i);
		   else
			   for(i=0;i<n;i++)
				   *(p+i*(n+1)) = *(pn+i);
			   
	   }

	   ass = 0;
       for(i=0;i<n;i++)
		   ass += (*(pn+i)-*(p0+i))*(*(pn+i)-*(p0+i));

	   for(i=0;i<2;i++)
		  kkk[i] = *(pk+i);

	   delete p0;
	   delete pn;
	   delete px;
	  // delete pk;   //??为什么不能delete
	  
      int bbs;
	  bbs = 8;
   }  

   for(i=0;i<n;i++)
	   cout<<endl<<*(p+i*(n+1));
       cout<<endl;
   delete e;
   delete p;
   
   
	}



};

void main()
{
	powell eg;

	eg.action(2,0.00001);
}