rosenbrock1.cpp

最优化算法

// RosenBrock1.cpp: implementation of the CRosenBrock class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "RosenBrock.h"
#include "RosenBrock1.h"
#include <math.h>
#include "StringResolution.h"


#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CRosenBrock::CRosenBrock(CString gongshi,CString x1,CString l1,CString &jieguo,float a1,float b1,float e1,int n1)
{
	j=1;
	k=1;
	a=a1;
	b=b1;
	e=e1;
	n=n1;
	d=new float*[n+1];
	p=new float*[n+1];
	q=new float*[n+1];
	l0=new float[n+1];
	l=new float[n+1];
	r=new float[n+1];
	resual=new float[n];
	jisuan.SetFormula(gongshi.GetBuffer(gongshi.GetLength()));


	for(int z=1;z<20;z++)//为x,y,l分配空间
	{
		x[z]=new float[n];
		y[z]=new float[n];		
	}
	string2(x1,x[1],n1);
	string2(l1,&l0[1],n1);

	for(int t=1;t<n+1;t++)//为d(n)分配空间//单位正交方向化
	{
		r[t]=0;
		p[t]=new float[n];
		q[t]=new float[n];
		d[t]=new float[n];	
		float * temp=d[t];
		for(int c=0;c<n;c++)//d的单位矩阵
			temp[c]=0;
		temp[t-1]=1;
		l[t]=l0[t];//l(j)=l0(j)
	}
    transform1(y[1],x[1],n);

	jieguo+="进行第一轮迭代:\r\n";
	themain(jieguo);
	jieguo+="最终结果:\r\n";
	for(int h=0;h<n;h++)
	{
	CString aa;
	aa.Format("x%d=%f\r\n",h,resual[h]);
	jieguo+=aa;
	}
	
}

CRosenBrock::~CRosenBrock()
{

}
// RosenBrock.cpp : Defines the entry point for the console application.



void CRosenBrock::themain(CString &jieguo)
{
	int leiji=0;//预防死循环,当leiji>=200时,跳出死循环
	do
	{
	    float *temp1=new float[n];//步骤2
	    transform1(temp1,y[j],d[j],l[j],n);
	    if(fun(temp1)<fun(y[j]))
	    {
		   transform1(y[j+1],temp1,n);
           r[j]+=l[j];//方向步长累加
		   l[j]=a*l[j];
		//////////////
	    }
	    else
	    {
		   transform1(y[j+1],y[j],n);
		   l[j]=b*l[j];
		/////////////////
	    }
        if(j<n) //步骤3
		{
			j=j+1;
			continue;
		}
	    if(fun(y[n+1])<fun(y[1])) //步骤4 (yes)
    	{
		   transform1(y[1],y[n+1],n);
		   j=1;
		   continue;
	    }
		 if(fun(y[n+1])==fun(y[1])) //步骤4 (no)
		 {
	         if(fun(y[n+1])<fun(x[k])) //步骤5 (yes)/////////////////////
			 {
		         transform1(x[k+1],y[n+1],n); //步骤6
				 float temp[2];
				 transform1(temp,x[k+1],x[k],-1,n);
				 if(fanshu(temp,n)<=e)
				 {
					 transform1(resual,x[k+1],n);
					 break;
				 }
				 else
				 {
				    trans(p,q,d,r,n);
					jieguo+="正交化d:\r\n";///////////显示中间信息 
					for(int w=1;w<=n;w++)
					{	CString aa;
	                    aa.Format("d(%d)=",w);
						float *temp=d[w];
						for(int y=0;y<n;y++)
						{
							CString bb;
	                         bb.Format("%f  ",temp[y]);
							 aa+=bb;
						}
						aa+="\r\n";
						jieguo+=aa;
					}////////////////////////////////
					transform1(y[1],x[k+1],n);
					k=k+1;
					CString cc;
					cc.Format("%d",k);
                    jieguo+="进行第"+cc+"轮迭代:\r\n";
					j=1;
					for(int t=1;t<n+1;t++)
					{
		               l[t]=1;
					   r[t]=0;
					}
					continue;
				 }
			 }
			 else                             //步骤5 (no)
			 {
				 int flag=1;///标记--判断是否|l[j]|<=e
				 for(int i=1;i<=n;i++)// |l[j]|<=e是否成立
				 {
					 if(fabsf(l[i])>e)
					 {
						 flag=0;
						 break;
					 }					 
				 }
				 if(flag==1)
				 {
					 transform1(resual,x[k],n);
					 break;
				 }
				 transform1(y[1],y[n+1],n);
				 j=1; 
				 continue;
			 }
		 }
	leiji++;
	}while(leiji<200);
}


void CRosenBrock::transform1(float *p,float *q,int count)//q数组的值赋值给q数组,如y(1)=x(1)
{
	for(int i=0;i<count;i++)
		p[i]=q[i];   
}
void CRosenBrock::transform1(float *temp,float *p,float *q,float l,int count)//公式:y(2)=y(1)+l*d(1)
{
	for(int i=0;i<count;i++)
		temp[i]=p[i]+l*q[i];
}
void CRosenBrock::transform1(float *temp,float *q,float l,int count)//公式:y(2)=y(2)+l*d(1)
{
	for(int i=0;i<count;i++)
		temp[i]+=l*q[i];
}
void CRosenBrock::transform3(float *temp,float *q,float l,int count)//公式:y(2)=l*d(1)
{
	for(int i=0;i<count;i++)
		temp[i]=l*q[i];
}
float CRosenBrock::transform2(float *p,float *q,int count)//公式:resual=y(2)*y(1)
{
	float temp=0;
	for(int i=0;i<count;i++)
		temp+=p[i]*q[i];
	return temp;
}
float CRosenBrock::fun(float *x1)
{
	float aa;
	//double *kk=new double[n];
	//aa=x[0]*x[0]+x[1]*x[1]-3*x[0]-x[0]*x[1]+3;
	aa=(x1[0]-3)*(x1[0]-3)+(x1[1]+2)*(x1[1]+2);
	//for(int k=0;k<n;k++)
	//	kk[k]=(double)x1[k];

	//aa=jisuan.computer(kk,jisuan.GetVariantTableSize());
    return aa;
}
void CRosenBrock::trans(float **p,float **q,float **d,float *r,int count)//构造新的正交方向
{
	float *temp_p;
	int j,i,k;
	for(j=1;j<=count;j++)       ///计算p(j)的值
	{
      temp_p=p[j];

	  if(r[j]==0)
	  {
		  transform1(temp_p,d[j],count);
	  }
	  else
	  {
		  for(k=0;k<count;k++)
			  temp_p[k]=0;

	      for(i=j;i<=count;i++)
	      {
			  transform1(temp_p,d[i],r[i],count);			  
		  }		  
	  }
	}
    for( j=1;j<=count;j++)       ///计算q(j)的值
	{
		if(j==1)
			transform1(q[j],p[j],count);
		else
		{
			float *temp_q=new float[count];
			float *temp_q1=new float[count];
			for(k=0;k<count;k++)
			  temp_q[k]=0;
			for(i=1;i<=j-1;i++)
			{
			    float temp;
			    temp=transform2(q[i],p[j],count)/transform2(q[i],q[i],count);
				transform1(temp_q,q[i],temp,count);
			}

			transform1(temp_q1,p[j],temp_q,-1,count);
			transform1(q[j],temp_q1,count);
		}
	}
     for( j=1;j<=count;j++)       ///计算d(j)的值
	 {
		 float *temp_q=q[j];
		 float aa=((float)1)/fanshu(q[j],count);
		 transform3(d[j],q[j],aa,count);
		 
	 }

}
float CRosenBrock::fanshu(float *x1,int count)//范数求值函数
{
   float temp=0;
	for(int i=0;i<count;i++)
		temp+=x1[i]*x1[i];
	return sqrtf(temp);
}


void CRosenBrock::string2(CString str1, float *zhi, int n0)
{
	int temp=0;
	CString str=str1;
	str.TrimLeft();
	str.TrimRight();
	str+=",";
		for(int e=0;e<n0;e++)      
		{
          zhi[e]=atof(str.Mid(temp,str.Find(",",temp)-temp));
          temp=str.Find(",",temp)+1;
		}

}