sumt.cpp

基于单纯形法的惩罚函数外点法的类。 在优化设计中使用。

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

	
	double fx(double *p)  //目标函数
	{ 
	     return(4*(*p-5)*(*p-5)+(*(p+1)-6)*(*(p+1)-6));
	     //return (*p+*(p+1)); 
	}

    double gg(double *p,int r)  //约束函数,p为优化变量,r为约束方程的标识。
	{
         double *g = new double [3];
         double some;
	     *(g+0) = *p**p-*(p+1);
	     *(g+1) = -*p;
	     *(g+2) = *p**p+20;

		 some = *(g+r);
         delete g;
	     return (some);
	}

/*double g1(double *p)
{
	
	return(*p**p-*(p+1));
}

double g2(double *p)
{

	return(-*p);
}  



typedef double (*pp)(double *p);

pp gg[2] = { g1, g2};  */



class sumt    //外点法惩罚函数优化方法类
{
private:
	int n,ng;   // n是维数，ng是优化变量的个数。
	double e,jing,M; // e,jing 是精度，M系数。
	double r;
	double xl,xh,xg,xc;  //单纯形法用到的：最好点，最坏点，次坏点。
	int f_xl,f_xh,f_xg,k;  //单纯形法用到的最好点，最坏点，次坏点的位置。
public:
	double daxiao(double a,double b)  //比较大小。
	{
		return (a>=b?a:b);
	}
	double fxy(double *p)   //构造外惩罚函数的无约束函数。
	{
		int i;
		double sum;
		double *m = new double [ng];
		for(i=0;i<ng;i++)
			*(m+i) = daxiao(gg(p,i),0);
		sum = 0;
		for(i=0;i<ng;i++)
			sum += *(m+i)**(m+i);
        delete m;
		return(fx(p)+M*sum);
	}

	void paixu(double *m)  //单纯形法中：排序求最好点，最坏点和次坏点。
	{
		int i;
		xl = *m;
		xh = *m;
		xg = *m;
		f_xl = 0;
		f_xh = 0;
		f_xg = 0;
		for(i=0;i<n+1;i++)
		{
			if(*(m+i)<xl)
			{
				xl = *(m+i);
				f_xl = i;
			}
			if(*(m+i)>xh)
			{
				xh = *(m+i);
				f_xh = i;
			}
		}
		xg = xl;
		f_xg = f_xl;
		for(i=0;i<n+1;i++)
		{
			if(xg<*(m+i)&&i!=f_xh)
			{
				xg = *(m+i);
				f_xg = i;
			}
		}
	}

	void xingxin(double *mm,double *cc) //单纯形法中求n个点的形心。
	{
		int i,j;
		for(i=0;i<n;i++)
			*(cc+i) = 0;
		for(i=0;i<n;i++)
		{
			for(j=0;j<n+1;j++)
				if(j!=f_xh)
					*(cc+i) += *(mm+i*(n+1)+j);
		}
		for(i=0;i<n;i++)
			*(cc+i) = *(cc+i)/n;
		xc = fxy(cc);
	}

	void dcx(int nn,double jing,double *z)  //单纯形法求无约束最优点。
	{
		int i,j;
		double h = 1.6;
		double ass;
		n = nn+1;   //设置构造单纯形的点的个数。
		k=0;
		double *p = new double [n*(n+1)];   //所有分量值
		double *e = new double [n*n];   //单位分量

		
  
     //  double temp[6];
	  // double temp2[3];

	   

		
		for(i=0;i<(n*n);i++) //构造单位分量。
		{
			if(i/n==i%n)
				*(e+i) = 1.0;
			else
				*(e+i) = 0.0;
		}

		

		for(i=0;i<n*(n+1);i+=n+1)   //设置初值为0
			*(p+i) = 0.0;

		
		for(i=0;i<n*(n+1);i++)
			if(i%(n+1)!=0)
				*(p+i) = *(p+i-i%(n+1))+h**(e+i-1-i/(n+1));
		delete e;

		
		
     ass = 2*(n+1)*jing*jing;  //保证进入循环

	 xl = ass;
	 

        double *w = new double [n];   //暂存点的坐标分量
		double *g = new double [n+1];  //各个点对应的函数值


     while(k++,sqrt(1.0/(n+1)*ass)>jing)   //  ?? sqrt(1.0/(n+1)*ass)
	 {

		

		for(j=0;j<n+1;j++)       //求每个点的函数值放入 g
		{
			for(i=0;i<n;i++)
				*(w+i) = *(p+i*(n+1)+j);
			*(g+j) = fxy(w);
		}

     /*  for(i=0;i<3;i++)
	         temp2[i] = *(g+i);*/



		paixu(g);    //对各个点进行排序，记录最优，最差，次差点。

		double *c = new double [n];  //形心点
		xingxin(p,c);
		xc = fxy(c);

		

		double *r = new double [n];  //反射点

		for(i=0;i<n;i++)
		{
			*(r+i) = *(c+i)+1.0*(*(c+i)-*(p+i*(n+1)+f_xh));
		}

        if(fxy(r)<xl)
		{
			double *ee = new double [n];
			for(i=0;i<n;i++)
			{
				*(ee+i) = *(c+i) + 2.0*(*(r+i)-*(c+i));
			}

			if(fxy(ee)<xl)
				for(i=0;i<n;i++)
					*(p+i*(n+1)+f_xh) = *(ee+i);
			else
				for(i=0;i<n;i++)
					*(p+i*(n+1)+f_xh) = *(r+i);
			delete ee;
			delete r;
			delete c;
		}
		else if(fxy(r)<xg)
		{
            for(i=0;i<n;i++)
					*(p+i*(n+1)+f_xh) = *(r+i);
			delete r;
			delete c;
		}
		else
		{
			if(fxy(r)<xh)
			{
			
			for(i=0;i<n;i++)
				*(p+i*(n+1)+f_xh) = *(c+i);
			}
		    
			
			double *s = new double[n];
			for(i=0;i<n;i++)
				*(s+i) = *(c+i)+0.5*(*(p+i*(n+1)+f_xh)-*(c+i));
			
			if(fxy(s)>xh)
			{
				for(i=0;i<n*(n+1);i++)
				{
					*(p+i) = *(p+(i/(n+1))+f_xl)+0.5*(*(p+i)-*(p+(i/(n+1))+f_xl));
				}
				delete s;

			}
			else
			{
                for(i=0;i<n;i++) 
				*(p+i*(n+1)+f_xh) = *(s+i);
				delete s;
			}
			delete r;
			delete c;
		}

      /* for(i=0;i<6;i++)
		   temp[i] = *(p+i);*/


         ass = 0;
		 for(j=0;j<n+1;j++)       //求每个点的函数值放入 g
		{
			for(i=0;i<n;i++)
				*(w+i) = *(p+i*(n+1)+j);
			*(g+j) = fxy(w);
		}
		for(i=0;i<n+1;i++)
		{
	       ass+=(*(g+i)-xc)*(*(g+i)-xc);
		}

  /*for(i=0;i<3;i++)
	  temp2[i] = *(g+i);*/

        
	

		
	 }
       

	 delete g;
	 delete w;
	 
     //cout<<endl;
	 //cout<<f_xl<<"   "<<f_xh<<"    "<<f_xg;
	 //cout<<endl;cout<<endl;
	 /*for(i=0;i<n-1;i++)
		 cout<<"X"<<i+1<<" = "<<*(p+i*(n+1)+f_xl)<<'\t'<<endl;
	 cout<<"叠代次数为"<<k<<endl;
*/
	 for(i=0;i<n-1;i++)
		 *(z+i) = *(p+i*(n+1)+f_xl);
	 delete p;
	 n--;

	}

	void action(int nn,double jingdu,int n2,double *jg)  //nn为维数，jingdu为精度，n2为约束方程个数，jg为优化结果。
	{
		int i,j;
		double R,ass,QQ;


        ng = n2;
		n = nn;
		j = 0;
		jing = jingdu;
		M = 1.0;
		R = 100;
		e = 0.0001;

		double *a = new double [n];  //a,b为迭代用的两个空间。
		double *b = new double [n];

	

		for(i=0;i<n;i++)
			*(a+i) = 1.0;

		while (++j)  
		{
			dcx(n,jing,b);

			QQ = gg(b,0);

			for (i=0;i<ng;i++)
			{ 
                if(gg(b,i)>QQ)
					QQ = gg(b,i);
			}
			if(QQ<=e)   
				break;

            if(M>R)
			{
				ass = 0;
				for(i=0;i<n;i++)
			        ass +=  (*(a+i)-*(b+i))*(*(a+i)-*(b+i));
				if(sqrt(ass)<jing)
					break;
			}
			
			M = M*8;
			for(i=0;i<n;i++)
				*(a+i) = *(b+i);

			if(j%10000==0)
				cout<<"迭代"<<j<<"次了"<<endl;
			if(j>200000)
				break;
		}
		for(i=0;i<n;i++)
			*(jg+i) = *(b+i);
		
		delete a;
		delete b;

		cout<<"迭代次数  "<<j<<endl;
	}


};

void main()
{   
	int i,n;
	n = 2;
    double *a = new double [n];
	sumt eg;
	
	eg.action(n,0.00000000001,3,a);  //第一参数为维数，第二参数为精度，第三参数为约束方程个数,a存储优化结果。
	
	for(i=0;i<2;i++)
			cout<<*(a+i)<<endl;

	delete a;
	
}