?? test.cpp
字號:
//www.pudn.com > 最優化作業C++源代碼.rar > 單純形法.cpp
/*#include<iostream.h>
#include<math.h>
int const n=2;
double fun(double x[n]);
void compare(double x[n+1][n],double LGH[3][n],double LGH_tap[3]);//在調用之前LGH先給x[n+1]賦值
double H(double x[n+1][n],double LGH[3][n],double c);//判別準則
void zx(double x[n+1][n],double LGH[3][n],double LGH_tap[3],double xn_1[n],double xn_2[n]);
void main()
{
double x[n+1][n]={0};
double c=0.1;
double a=1.1;
double b=0.5;
double LGH[3][n]={0};
double LGH_tap[3]={0};
double xn_1[n]={0};
double xn_2[n]={0};
double xn_3[n]={0};
double xn_4[n]={0};
int i,j,tap=0;
//給出一個初始點X0,X!,X2
x[0][0]=0;x[0][1]=0;
x[1][0]=0.965;x[1][1]=0.259;
x[2][0]=0.259;x[2][1]=0.965;
cout<<"初始點："<<endl;
for(i=0;i<n+1;i++)
{
cout<<"x["<<i<<"]=(";
for(j=0;j<n;j++)
{
cout<<x[i][j];
if(j!=n-1)
cout<<",";
}
cout<<")"<<endl;
}
compare(x,LGH,LGH_tap);//返回Xl,Xh,Xg及它們在x[n+1][n]中的下標
zx(x,LGH,LGH_tap,xn_1,xn_2);//返回中心X(n+1)和反射點X(n+2)
while(H(x,LGH,c)==0)
{
if(fun(xn_2)<fun(LGH[0]))//if(fun(xn_2)<fun(LGH[0][n]))//////////////////////////////////////////////////?關于傳值疑問
{
for(i=0;i<n;i++)
{
xn_3[i]=xn_1[i]+a*(xn_2[i]-xn_1[i]);//擴張
}
if(fun(xn_3)<fun(LGH[0]))
{
for(i=0;i<n;i++)
{
x[int(LGH_tap[2])][i]=LGH[2][i]=xn_3[i];
}
//goto step 2;
}
else
{
for(i=0;i<n;i++)
{
x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];
}
//goto step 2;
}
}
else
{
if(fun(xn_2)<fun(LGH[1]))
{
for(i=0;i<n;i++)
{
xn_4[i]=xn_1[i]+b*(xn_2[i]-xn_1[i]);//壓縮
}
//goto step 2;
}
else
{
for(i=0;i<n;i++)
{
xn_4[i]=xn_1[i]-b*(xn_2[i]-xn_1[i]);//壓縮
}
}
if(fun(xn_4)<fun(LGH[2]))
{
for(i=0;i<n;i++)
{
x[int(LGH_tap[2])][i]=LGH[2][i]=xn_4[i];
}
}
else//縮邊
{
cout<<"進行縮邊。"<<endl;
for(i=0;i<n+1;i++)
{
if(i==LGH_tap[0])
continue;
for(j=0;j<n;j++)
{
x[i][j]=(LGH[0][j]+x[i][j])/2;
}
}
}
}
compare(x,LGH,LGH_tap);//返回Xl,Xh,Xg及它們在x[n+1][n]中的下標
zx(x,LGH,LGH_tap,xn_1,xn_2);//返回中心X(n+1)和反射點X(n+2)
tap++;
}
cout<<"單純形法求解minf(x)=X1^2+2X2^2-4X1-8X2+5;"<<endl;
cout<<"XL="<<endl;
for(i=0;i<n;i++)
{
cout<<LGH[0][i]<<endl;
}
cout<<"fun="<<fun(LGH[0])<<endl;
cout<<"迭代次數:"<<tap<<endl;
}
double fun(double x[n])//求函數值
{
double y;
y=pow(x[0],2)+2*pow(x[1],2)-4*x[0]-8*x[1]+5;
//cout<<"y="<<y<<endl;
return y;
}
void compare(double x[n+1][n],double LGH[3][n],double LGH_tap[3])//在調用之前LGH先給x[n+1]賦值
{
int i;
double y[n+1];
for(i=0;i<n+1;i++)
y[i]=fun(x[i]);//y[i]=fun(x[i][n]);//////////////////////////////////////////////////?關于傳值疑問
double temp[2];
temp[0]=y[0];
for(i=0;i<n+1;i++)//temp[0]存儲最小的值XL,temp[1]存儲對應x的下標
{
if(y[i]<=temp[0])
{
temp[0]=y[i];
temp[1]=i;
}
}
LGH_tap[0]=temp[1];//記錄標記
for(i=0;i<n;i++)
LGH[0][i]=x[int(temp[1])][i];
for(i=0;i<n+1;i++)//temp[0]存儲最大的值XH,temp[1]存儲對應x的下標
{
if(y[i]>=temp[0])
{
temp[0]=y[i];
temp[1]=i;
}
}
LGH_tap[2]=temp[1];//
for(i=0;i<n;i++)
LGH[2][i]=x[int(temp[1])][i];
temp[0]=y[0];
if(LGH_tap[2]==0)
{
temp[0]=y[1];
}
for(i=0;i<n+1;i++)//尋找次大值
{
if(i==LGH_tap[2])//if(i==temp[1])
continue;
if(y[i]>=temp[0])
{
temp[0]=y[i];
temp[1]=i;
}
}
LGH_tap[1]=temp[1];
for(i=0;i<n;i++)
LGH[1][i]=x[int(temp[1])][i];
}
void zx(double x[n+1][n],double LGH[3][n],double LGH_tap[3],double xn_1[n],double xn_2[n])
{
int i,j;
for(i=0;i<n;i++)
xn_1[i]=xn_2[i]=0;
for(i=0;i<n+1;i++)//返回重心
{
if(i==LGH_tap[2])//輪流加n次（去掉一次）
continue;
for(j=0;j<n;j++)
{
xn_1[j]+=x[i][j];
}
}
for(i=0;i<n;i++)
xn_1[i]/=n;
for(i=0;i<n;i++)//返回反射點
xn_2[i]=2*xn_1[i]-LGH[2][i];
}
double H(double x[n+1][n],double LGH[3][n],double c)//判別準則
{
int i;
double s=0;
for(i=0;i<n+1;i++)//小心i從0到n即(i=0;i<n+1;i++)
{
s+=pow(fun(x[i])-fun(LGH[0]),2);//s+=pow(fun(x[i][n])-fun(LGH[0][n]),2);//////////////////////////////////////////////////?關于傳值疑問
}
if(s<=c)
return 1;
else
return 0;
}
*/
//Part II
//www.pudn.com > 最優化作業C++源代碼.rar > DFP法.CPP
//
// #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++)
// {
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -