optimumt1.c

这是一个高通滤波器的最优化算法。有约束问题采用的是惩罚函数法

//This program is used to optimize a high Pass Filter.
//Optimization can be implemented through External_Point_Punishment Method + Nelder and Mead's Simplex method

#define N 4	//the number of variables in Simplex method
#define pi 3.1415927
#include "stdio.h"
#include "math.h"

double r=2;	//Default Value of the Punishment Coefficient
int k=1;    //Times of iteration

//----------------------Function Declaration----------------------
double get_Q(double x[]);                   //Get quality factor Q
double target_function(double x[]);	      //Target Function F(x)
double punishment_function(double x[]);   //Punishment Function //Used in External_Point_Punishment Method 
double punishment_detect(double x[]);         //To See if the Optimization can go to the end
void simplex(double x0[],double x_def[],double t,double alpha,double beta,double gamma,double eps);                       //Simplex Method
void Freq_Res(double x[]);	               //Print Frequency Response
void display_x(double x[]);                  //print Values of X
//----------------------Function Declaration----------------------

void main()
{

	double x_def[N]={4,2.0,30,15};
	double eps=9E-30;

	printf("Default Value is:\n");
    display_x(x_def);
	printf("Frequency Response Before Optimization:\n");
	Freq_Res(x_def);	
	printf("Punishment(x)=%e\n",punishment_function(x_def));
	printf("Quanlity Factor is:  Q=%e\n",get_Q(x_def));
	getchar();//Display cover by cover

	do{	simplex(x_def,x_def,0.3,1,0.5,2,0.1);
		r=r*5;
		k=k+1;
	}while (punishment_detect(x_def)>=eps);

	printf("Values After Optimazation is:\n");
    display_x(x_def);
	printf("Frequency Response After Optimization:\n");
	Freq_Res(x_def);	
	printf("punishment(x)=%e\n",punishment_function(x_def));
	printf("Quanlity Factor is:  Q=%e\n",get_Q(x_def));
	printf("Times of iteration is: %d\n",k);
	getchar();//Display cover by cover
}
double get_Q(double x[])
{	double t;
        t=1/(2-x[3]/x[2]);
	//t=sqrt(x[0]*x[1]*x[2])/((x[0]+x[1])*sqrt(x[3]));
	return t;
	printf("\n");
}
double get_sqr(double x)
{return x*x;}

double target_function(double x[])
{	
	double wk[18]={4.0,3.5,3.2,3.0,2.7,2.5,2.0,1.90,1.70,1.60,1.50,1.40,1.30,1.20,1.10,1.00,1.00,1.00};	
        double f1[18]={15000,15050,15100,15200,15400,15800,16000,17000,
                       19000,21000,23000,25000,27000,29000,31000,33000,36000,40000};
	double wj[10]={3.5,3.0,2.7,2.5,2.0,1.6,1.3,1,1,1};	
	double f2[10]={14950,14900,14800,14400,14000,13000,
                       11000,9000,6000,1000};
	double h;
	double t;
	int i;

	t=0;
	for(i=0;i<18;i++)
	{		
		h=((1+x[3]/x[2])*get_sqr(2*pi*f1[i]*x[0]*x[1]*1E-6))/sqrt(get_sqr(1-get_sqr(2*pi*f1[i]*x[0]*x[1]*1E-6))+get_sqr(2*pi*f1[i]*(2-x[3]/x[2])*x[0]*x[1]*1E-6));
		//h=1/sqrt(get_sqr(1-get_sqr(2*pi*f1[i])*x[0]*x[1]*x[2]*x[3]*1E-12)+get_sqr(2*pi*f1[i]*(x[0]+x[1])*(1E-6)*x[3]));
		t=t+wk[i]*get_sqr(h-1.4);
	}
	for(i=0;i<10;i++)
	{	
		h=((1+x[3]/x[2])*get_sqr(2*pi*f2[i]*x[0]*x[1]*1E-6))/sqrt(get_sqr(1-get_sqr(2*pi*f2[i]*x[0]*x[1]*1E-6))+get_sqr(2*pi*f2[i]*(2-x[3]/x[2])*x[0]*x[1]*1E-6));	
		//h=1/sqrt(get_sqr(1-get_sqr(2*pi*f2[i])*x[0]*x[1]*x[2]*x[3]*1E-12)+get_sqr(2*pi*f2[i]*(x[0]+x[1])*(1E-6)*x[3]));
		t=t+wj[i]*get_sqr(h);
	}
	return t;
}

double punishment_function(double x[4])
{
	double t;
	
	t=target_function(x);
	if (x[0]-10<0)		t=t+r*get_sqr(x[0]-10);
	if (x[1]-0.1<0)		t=t+r*get_sqr(x[1]-0.1);
	if (x[2]-10<0)		t=t+r*get_sqr(x[2]-10);
	t=t+r*get_sqr(((1+x[3]/x[2])*get_sqr(2*pi*15000*x[0]*x[1]*1E-6))/sqrt(get_sqr(1-get_sqr(2*pi*15000*x[0]*x[1]*1E-6))+get_sqr(2*pi*15000*(2-x[3]/x[2])*x[0]*x[1]*1E-6))-1.4*0.707);
	//t=t+r*get_sqr(sqrt(get_sqr(1-get_sqr(2*pi*10000)*x[0]*x[1]*x[2]*x[3]*1E-12)+get_sqr(2*pi*10000*(x[0]+x[1])*(1E-6)*x[3]))-1.414);
	return t;
}
double punishment_detect(double x[4])
{	
	double t;
	
	t=0;
	if (x[0]-10<0)		t=t+get_sqr(x[0]-10);
	if (x[1]-0.1<0)		t=t+get_sqr(x[1]-0.1);
	if (x[2]-10<0)		t=t+get_sqr(x[1]-10);
	t=t+get_sqr(((1+x[3]/x[2])*get_sqr(2*pi*15000*x[0]*x[1]*1E-6))/sqrt(get_sqr(1-get_sqr(2*pi*15000*x[0]*x[1]*1E-6))+get_sqr(2*pi*15000*(2-x[3]/x[2])*x[0]*x[1]*1E-6))-1.4*0.707);
	//t=t+get_sqr(sqrt(get_sqr(1-get_sqr(2*pi*10000)*x[0]*x[1]*x[2]*x[3]*1E-12)+get_sqr(2*pi*10000*(x[0]+x[1])*(1E-6)*x[3]))-1.414);
	return t;
}
void Freq_Res(double x[])	
{
	double f[28]={1000,6000,9000,11000,13000,14000,14400,14800,14900,14950,15000,15050,15100,15200,15400,
	15800,16000,17000,19000,21000,23000,25000,27000,29000,31000,33000,36000,40000};
	double h;
	int i;
	for(i=0;i<28;i++)
	{	
		h=((1+x[3]/x[2])*get_sqr(2*pi*f[i]*x[0]*x[1]*1E-6))/sqrt(get_sqr(1-get_sqr(2*pi*f[i]*x[0]*x[1]*1E-6))+get_sqr(2*pi*f[i]*(2-x[3]/x[2])*x[0]*x[1]*1E-6));	
		//h=1/sqrt(get_sqr(1-get_sqr(2*pi*f[i])*x[0]*x[1]*x[2]*x[3]*1E-12)+get_sqr(2*pi*f[i]*(x[0]+x[1])*(1E-6)*x[3]));
		printf("A(%f)=%f    ",f[i],h);
		if((i+1)%3==0)
			printf("\r");
	}
	printf("\n");
}
void display_x(double x[])
{
	int i;
	for(i=0;i<4;i++)
		printf("%e ",x[i]);
	printf("\n");
}

void simplex(double x0[],double x_def[],double t,double alpha,double beta,double gamma,double eps)
//X0-initial point value, t-side length of simplex, alpha-reflection coefficient
//beta-shrink coefficient,gamma-expansion coefficient
{
	double x[N+1][N];	//N+1 vertexs of the simplex
	double d1,d2;
	double tem1,tem2,tem3,tem4,tem5;
	double xc[N],xr[N],xe[N],xt[N];//Centroid,Reflection,Expansion,shrink point
	int h,g,l;
	int i,j;

	//Get N+1 vertexs of the simplex
	d1=t/(N*sqrt(2))*(sqrt(N+1)+N-1);
	d2=t/(N*sqrt(2))*(sqrt(N+1)-1);
	for (j=0;j<N;j++)
		x[0][j]=x0[j];
	for (i=1;i<N+1;i++)
		for (j=0;j<N;j++)
			if ((i-1)==j)
				x[i][j]=x0[j]+d1;
			else
				x[i][j]=x0[j]+d2;
	do
	{
		//Get Xh,Xg and Xl
		tem1=tem2=tem5=punishment_function(x[0]);
		h=g=l=0;
		for(i=1;i<N+1;i++)
		{
			tem3=punishment_function(x[i]);
			if(tem1<tem3)	
			{
				tem5=tem1;g=h;
				tem1=tem3;h=i;
			}
			else if(tem5<tem3)	
			{
				tem5=tem3;g=i;
			}

			tem4=punishment_function(x[i]);
			if(tem2>tem4)	
			{
				tem2=tem4;
				l=i;
			}
		}
		//Reflection
		//Get Xc
		for(i=0;i<N;i++) xc[i]=0;
		for(i=0;i<N+1;i++)
			if(i!=h)
				for (j=0;j<N;j++)
					xc[j]=xc[j]+x[i][j];
		for(i=0;i<N;i++) xc[i]=xc[i]/N;
		for (i=0;i<N;i++)
			xr[i]=xc[i]+alpha*(xc[i]-x[h][i]);
		
		tem4=punishment_function(xr);
		if (tem4<tem2)	//punishment_function(Xr)<punishment_function(Xl)
		{
			//expansion
			for(i=0;i<N;i++)
				xe[i]=xc[i]+gamma*(xr[i]-xc[i]);
			if(punishment_function(xe)<tem2)	//punishment_function(Xe)<punishment_function(Xl)=>Xh=Xe
				for (i=0;i<N;i++) x[h][i]=xe[i];
			else	//Xh=Xr
				for (i=0;i<N;i++) x[h][i]=xr[i];
		}
		else
			if(tem4<punishment_function(x[g]))	//punishment_function(Xr)<punishment_function(Xg)=>Xh=Xr;
				for (i=0;i<N;i++) x[h][i]=xr[i];
			else //Contruction
			{
				if(tem4<=tem1)	//punishment_function(Xr)<punishment_function(Xh)
					for(i=0;i<N;i++)
						xt[i]=xc[i]+beta*(xr[i]-xc[i]);
				else
					for(i=0;i<N;i++)
						xt[i]=xc[i]+beta*(x[h][i]-xc[i]);
				if (punishment_function(xt)<=tem1)
					//Get new Simplex:Xh=Xt
					for (i=0;i<N;i++) x[h][i]=xt[i];
				else
					//Construction failed
					//Start to shrink whole simplex
					for(i=0;i<N+1;i++)
						if(i!=l)
							for(j=0;j<N;j++)
								x[i][j]=x[l][j]+0.5*(x[i][j]-x[l][j]);
			}
		tem1=0;tem2=punishment_function(xc);
		for(i=0;i<N+1;i++)
			tem1=tem1+(punishment_function(x[i])-tem2)*(punishment_function(x[i])-tem2);
		tem1=tem1/(N+1);
	}
	while(sqrt(tem1)>=eps);
	for(i=0;i<N;i++)
		x_def[i]=xc[i];
}