dfp法.cpp

用C++编的一些最优化作业中的程序,有Newton法,DFP法,共轭梯度法,单纯形法,内点法,外点法,内外点法,都能使用,我已经全部调试过了

#include<iostream.h>
#include<math.h>
#include<iomanip.h>
int const n=2;//正定二次函数的自变量个数
double fun(double x[n],double f_xs[n+n+1+(n-1)*n/2]);//输入变量为函数自变量初值
void Q_fun(double f_xs[n+n+1+(n-1)*n/2],double Q[n][n+1]);//Q[n][n]是二次函数的正定矩阵，但Q的第n+1列存储一次项的系数
void D_fun(double x[n],double Q[n][n+1],double g0[n]);//自变量初值，正定二次函数的各项系数，返回梯度的初值
int H(double g0[n],double c);//判别准则：返回1结束，返回0继续迭代
void abc(double x[n],double p[n],double f_xs[n+n+1+(n-1)*n/2],double t[3]);//t[3]中返回的是a,b,c的系数值.开始计算minf(Xk+tPk)时的步长t的值，由于这是n元二次函数所以minf(t)是关于t>0的二次函数，先求二次方程a,b,c系数,用一阶导为零。
double t_c(double t[3]);//二次函数一阶导为零计算t的值，t>0
void main()
{
	double f_xs[n+n+1+(n-1)*n/2]={4,1,-40,-12,136,0};//正定二次函数的各项系数，第一个为：X1^2系数，第二个为：X2^2系数，第三个为：X1系数，第四个为:X2系数，第五个为：常数项，第六个为：X1X2系数；
	//n元的系数存放类推
	double x[n]={8,9};//函数自变量初值
	double f0;//函数值
	double g0[n];//梯度的值
	double Q[n][n+1];//求梯度处值设置的中间变量，包含两部分：正定二次函数对应的实对称矩阵，还有一次项系数
	double c=0.01;//H准则限值
	double t[3];//返回求minf()时t的二次函数的a,b,c的系数值
	double t_bc;//步长
	double p[n];//保存下降方向
	double H0[n][n];//保存模拟Hesse矩阵的逆
	double y[n];//y(k)=g0(k+1)-g0(k)
	double s[n];//s(k)=X(k+1)-X(k)
	double s_temp[n][n]={0};//计算保存矩阵
	double s_temp2[n][n]={0};
	double s_temp3[n][n]={0};
	double s_tl[n]={0};
	double temp;//临时值
	int i,j,k,flag=0,tap=0;//迭代次数

	Q_fun(f_xs,Q);//计算正定二次函数对应的实对称矩阵
	f0=fun(x,f_xs);//求函数初值
	D_fun(x,Q,g0);//返回梯度的初值

	do
	{
		for(i=0;i<n;i++)//给H0[n][n]的处值赋单位矩阵
		{
			for(j=0;j<n;j++)
			{
				if(i==j)
					H0[i][j]=1;
				else
					H0[i][j]=0;
			}
		}
		for(i=0;i<n;i++)
			p[i]=(-1)*g0[i];
		k=0;//step 2;

		do
		{
			abc(x,p,f_xs,t);//开始计算minf(Xk+tPk)时的步长t的值，
			t_bc=t_c(t);//求一阶导来计算t

			for(i=0;i<n;i++)
			{
				x[i]=x[i]+t_bc*p[i];
				s[i]=t_bc*p[i];//保存计算之值X(k+1)-X(k)
			}
			for(i=0;i<n;i++)
				y[i]=g0[i];//保存之类的

			f0=fun(x,f_xs);
			D_fun(x,Q,g0);//step 3;

			for(i=0;i<n;i++)
				y[i]=g0[i]-y[i];//保存计算g0(k+1)-g0(k)

			if(H(g0,c)==0)//即不满足小于c
			{
				if(k!=n)
				{
					//y
					//s   have done!
					temp=0;//初值
					for(i=0;i<n;i++)
						temp+=s[i]*y[i];
					for(i=0;i<n;i++)
					{
						for(j=0;j<n;j++)
						{
							s_temp[i][j]=s[i]*s[j]/temp;
						}
					}//求出S(k)S(k)t/S(k)tyk

					for(i=0;i<n;i++)//初值
						s_tl[i]=0;

					for(i=0;i<n;i++)
					{
						for(j=0;j<n;j++)
						{
							s_tl[i]+=y[j]*H0[j][i];
						}
					}
					temp=0;//初值
					for(i=0;i<n;i++)
					{
						temp+=s_tl[i]*y[i];
					}//这时s_tl[n]和s_temp2[n][n]可以利用

					for(i=0;i<n;i++)//初值
						s_tl[i]=0;

					for(i=0;i<n;i++)
					{
						for(j=0;j<n;j++)
						{
							s_tl[i]+=H0[i][j]*y[j];
						}
					}

					for(i=0;i<n;i++)//初值
						for(j=0;j<n;j++)
							s_temp2[i][j]=0;

					for(i=0;i<n;i++)
					{
						for(j=0;j<n;j++)
						{
							s_temp2[i][j]=s_tl[i]*y[j];
						}
					}

					for(i=0;i<n;i++)//初值
						for(j=0;j<n;j++)
							s_temp3[i][j]=0;

					for(i=0;i<n;i++)
					{
						for(j=0;j<n;j++)
						{
							for(int ii=0;ii<n;ii++)
							{
								s_temp3[i][j]+=s_temp2[i][ii]*H0[ii][j];
							}
						}
					}

					for(i=0;i<n;i++)
					{
						for(j=0;j<n;j++)
						{
							H0[i][j]=H0[i][j]+s_temp[i][j]-s_temp3[i][j]/temp;
						}
					}
					for(i=0;i<n;i++)
					{
						temp=0;
						for(j=0;j<n;j++)
						{
							temp+=H0[i][j]*g0[j];
						}
						p[i]=(-1)*temp;
					}
					k++;
				}
				else
				{
					f0=fun(x,f_xs);
					//这里x[n],g0[n]的值已经修改
					break;
				}
			}
			else
			{
				flag=1;
				break;
			}
			tap++;
		}while(H(g0,c)==0);
		if(flag==1)//符合梯度准则跳出
		{
			flag=0;
			break;
		}
	}while(k==n);

	cout<<"DFP法"<<endl;
	cout<<"函数f(x1,x2)=4(X1-5)^2+(X2-6)^2.的极小点为："<<"f="<<f0<<endl;
	cout<<"自变量取值为："<<endl;
	for(i=0;i<n;i++)
	{
		cout<<"x"<<i+1<<"="<<x[i]<<endl;
	}
	cout<<"迭代次数为："<<tap<<endl;
}
double fun(double x[n],double f_xs[n+n+1+(n-1)*n/2])//输入变量为函数自变量初值
{
	int i,j;
	double f=0;
	for(i=0;i<n;i++)//计算X^2部分
	{
		f+=pow(x[i],2)*f_xs[i];
	}
	for(;i<2*n;i++)//计算X部分
	{
		f+=x[i%n]*f_xs[i];
	}
	f+=f_xs[i];//计算常数项部分
	for(i=0;i<n;i++)//计算XiXj部分
	{
		for(j=i+1;j<n;j++)
		{
			f+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*x[i]*x[j];
		}
	}
	return f;
}
void Q_fun(double f_xs[n+n+1+(n-1)*n/2],double Q[n][n+1])//Q[n][n]是二次函数的正定矩阵，但Q的第n+1列存储一次项的系数
{
	int i,j;
	for(i=0;i<n;i++)
	{
		Q[i][i]=2*f_xs[i];
	}
	for(i=0;i<n;i++)
	{
		Q[i][n]=f_xs[n+i];
	}
	for(i=0;i<n;i++)
	{
		for(j=i+1;j<n;j++)
		{
			Q[j][i]=Q[i][j]=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1];
		}
	}
}
void D_fun(double x[n],double Q[n][n+1],double g0[n])//自变量初值，正定二次函数的各项系数，返回梯度的初值
{
	int i,j;
	for(i=0;i<n;i++)
	{
		g0[i]=0;//清零
		for(j=0;j<n;j++)
		{
			if(i==j)
				g0[i]+=Q[i][j]*x[j];
			else 
				g0[i]+=Q[i][j]*x[j];
		}
	}
	for(i=0;i<n;i++)
	{
		g0[i]+=Q[i][n];
	}
	cout<<endl<<"梯度g0的值="<<endl;
	for(i=0;i<n;i++)
		cout<<setw(10)<<g0[i]<<endl;
	cout<<endl;
}
int H(double g0[n],double c)//判别准则：返回1结束，返回0继续迭代
{
	double s=0;
	for(int i=0;i<n;i++)
	{
		s+=pow(g0[i],2);
	}
	if(sqrt(s)<c)
		return 1;
	else 
		return 0;
}
void abc(double x[n],double p[n],double f_xs[n+n+1+(n-1)*n/2],double t[3])//t[3]中返回的是a,b,c的系数值
{
	int i,j;
	
	t[0]=t[1]=t[2]=0;
	t[0]=fun(p,f_xs)-f_xs[2*n];
	for(i=n;i<2*n;i++)
	{
		t[0]-=f_xs[i]*p[i%n];
	}	
	for(i=0;i<n;i++)
	{
		t[1]+=2*f_xs[i]*x[i]*p[i];
	}	
	for(;i<2*n;i++)
	{
		t[1]+=f_xs[i]*p[i%n];
	}
	for(i=0;i<n;i++)
	{
		for(j=i+1;j<n;j++)
		{
			t[1]+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*(x[i]*p[j]+x[j]*p[i]);
		}
	}
	t[2]=fun(x,f_xs);
	cout<<endl<<"二次函数minf(Xk)的二次项系数："<<t[0]<<" 一次项系数："<<t[1]<<" 常数项系数："<<t[2]<<endl<<endl;	
}
double t_c(double t[3])
{
	cout<<"用一阶导求得步长为：t="<<-t[1]/t[0]/2<<endl;
	return -t[1]/t[0]/2;
}