#include<iostream.h>
#include<math.h>
double lamta[10]={0, 1.0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,1};//鮑威爾方法初始化方向，線性無關
double lamta1[3]={0, 0 , 0};//暫存新的搜索方向
double x1[4]={0, 0 ,0, 0 };//x1到x3用于存儲各共軛方向的點
double x2[4]={0, 0 ,0, 0 };
double x3[4]={0, 0 ,0, 0 };
double x4[4]={0, 0 ,0, 0 };//x4用于中間判斷
double x5[4]={0, 0 ,0, 0 };//x5用存放于更換方向后產生的新點
int m=0;//標志
double x_[4]={0, 0, 0, 0};//暫存鮑威爾最優解
double x0[4]={0, 2, 2 , 2};//初值 
double c=10;//遞減系數
double e=0.00001;//精度控制
double r0=1;//初始懲罰因子
double r=1;
//函數聲明部分
void Powell(double r);   //鮑威爾方法函數
double fxy(double x1,double x2,double x3,double r);  //待求函數
double ysearch(double x);  //一維搜索的目標函數
void search(double &a,double &b,double h);   //區間搜索
double  yellowcut(double &a,double &b);  //黃金分割
void sort(double *p,int size);//選擇法排序
void main()                     //約束優化方法主函數入口
{
 cout<<"請輸入精度"<<endl;
  cin>>e;
changyan:Powell(r);
 double cmpare[4];
 int flag1=0;
 for (int i=1;i<=3;i++)
 {
  cmpare[i]=x_[i]-x0[i];
  if (fabs(cmpare[i])<e)
  {flag1++;}
 }
  if (flag1==3)
  {
   cout<<"最優解為："<<"x1="<<x_[1]<<"      "<<"x2="<<x_[2]<<"     "<<"x3="<<x_[3]<<endl;
   cout<<"最小值為"<<fxy(x_[1],x_[2],x_[3],r)<<endl;
  }
  else
  {
   for (int j=1;j<=3;j++)
   {
    x0[j]=x_[j];
   }
   r=c*r;
   goto changyan;
  }
}
//子函數定義部分
double fxy(double x1,double x2,double x3,double r)//待求函數
{
 double m,n,p;
 m=(-x1>0)?(-x1):0;
 n=(-x2>0)?(-x2):0;
 p=(-x3>0)?(-x3):0;
 return    //懲罰函數
1000-x1*x1-2*x2*x2-x3*x3-x1*x2-x1*x3+r*(m*m+n*n+p*p)+r*((x1*x1+x2*x2+x3*x3-25)*(x1*x1+x2*x2+x3*x3-25)+(8*x1+14*x2+7*x3-56)*(8*x1+14*x2+7*x3-56));
}
void Powell(double r)     //鮑威爾方法函數定義
{
  double det=0.0001;      //迭代精度
 int k;
my1: for (k=1;k<=3;k++)
  {
   m=3*k-2;
   double a=0,b=0,xo=0;
   search(a,b,1);  //完成區間搜索
   double temp;
   temp=yellowcut(a,b);//黃金分割法
   int n=3*k-2;
   for (int i=1;i<=3;i++)
   {
    switch (k)
    {
    case 1:x1[i]=x0[i]+temp*lamta[n++];break;
    case 2:x2[i]=x1[i]+temp*lamta[n++];break;
    case 3:x3[i]=x2[i]+temp*lamta[n++];break;
    default :break;
    }
   }
  }
  double cmp[4];
  int flag=0;
  for (int i=1;i<=3;i++)
  {
   cmp[i]=x3[i]-x0[i];
   if (fabs(cmp[i])<det)
   {flag++;}}
  if (flag==3)     //找到最優解
  {
   x_[1]=x3[1];
   x_[2]=x3[2];
   x_[3]=x3[3];
   
  }
  else
  {
   double fy[4];
   fy[0]=fxy(x0[1],x0[2],x0[3],r);
   fy[1]=fxy(x1[1],x1[2],x1[3],r);
   fy[2]=fxy(x2[1],x2[2],x2[3],r);
   fy[3]=fxy(x3[1],x3[2],x3[3],r); double fyy[3];
   for (int ii=0;ii<3;ii++)
   {fyy[ii]=fy[ii]-fy[ii+1];}
   sort(fyy,3);
   for (ii=1;ii<=3;ii++)
   {x4[ii]=2*x3[ii]-x0[ii];}
   double f0,f3,f4;
   f0=fy[0];
   f3=fy[3];
   f4=fxy(x4[1],x4[2],x4[3],r);
   if ((f0+f4-2*f3)/2>=fyy[2])
   {
    if (f3<f4)
    {
     for (int t=1;t<=3;t++)
     {x0[t]=x3[t];}
    }
    else
    {
     for (int t=1;t<=3;t++)
     {x0[t]=x4[t];
    }}
    goto my1;
   }
   else{
    for (int t=0;t<3;t++)
    {lamta1[t]=x3[t+1]-x0[t+1];}
    m=0;                   //switch 標志！
    double aa=0,bb=0;
    search(aa,bb,1);
    double temp1;
    temp1=yellowcut(aa,bb);
    for (int i=1;i<=3;i++)
    {x5[i]=x3[i]+temp1*lamta1[i-1];}
    for (i=1;i<=3;i++)
    {x0[i]=x5[i];}
    for (i=1;i<=6;i++)
    {lamta[i]=lamta[i+3];}
    for (i=1;i<=3;i++)
    {
     lamta[6+i]=lamta1[i-1];}
    goto my1;
}}}

double ysearch(double x)  //一維搜索的目標函數
{
switch (m)
{
case 1: return fxy(x0[1]+x*lamta[m],x0[2]+x*lamta[m+1],x0[3]+x*lamta[m+2],r);break;
case 4: return fxy(x1[1]+x*lamta[m],x1[2]+x*lamta[m+1],x1[3]+x*lamta[m+2],r);break;
case 7: return fxy(x2[1]+x*lamta[m],x2[2]+x*lamta[m+1],x2[3]+x*lamta[m+2],r);break;
case 0: return fxy(x3[1]+x*lamta1[0],x3[2]+x*lamta1[1],x3[3]+x*lamta1[2],r);break;//更改方向后的一維搜索
default:return 0; break;
}
}
void search(double &a,double &b,double h)   //區間搜索
{double a1,a2,a3,y1,y2,y3;
h=1;
a1=a,y1=ysearch(a1);
a2=a+h,y2=ysearch(a2);
if(y2>=y1){
 h=-h,a3=a1,y3=y1;
 a1=a2,y1=y2,a2=a3,y2=y3;}
a3=a2+h,y3=ysearch(a3);
while(y3<=y2){
 h=2*h;
 a1=a2,y1=y2,a2=a3,y2=y3;
 a3=a2+h,y3=ysearch(a3);
}
if(h<0)a=a3,b=a1;
else a=a1,b=a3;}
double  yellowcut(double &a,double &b){
 double e;                                        //黃金分割法求解
 e=0.001;
 double c,fc;
 c=a+0.382*(b-a);
 fc=ysearch(c);
 double d,fd;
 double xo;
 d=a+0.618*(b-a);
 fd=ysearch(d);
label2: if (fc<=fd)
  {b=d;
  d=c;
  fd=fc;
  c=a+0.382*(b-a);
  fc=ysearch(c);}
  else
  {a=c;
  c=d;
  fc=fd;
  d=a+0.618*(b-a);
  fd=ysearch(d);}
  if ((b-a)<=e)
  {xo=(a+b)/2;}
  else
   goto label2;
  return xo;
}
void sort(double *p,int size){//選擇法排序
 int i,j;
 double k;
 for(i=0;i<size-1;i++)
  for(j=i+1;j<size;j++)
   if(*(p+i)>*(p+j)){k=*(p+i);*(p+i)=*(p+j);*(p+j)=k;}
}