最优法与黄金分割

      优选法应用在我国从70年代初，华罗庚先生最初研究推广并大量应用。对编程领域某些解决问题的方法有所帮助。其实数学上的很多方法可以解决一个程序时间和空间最优复杂度的问题。希望对大家有用……
　　优选法,是以数学原理（即中国的古代黄金分割）为指导，用最可能少的试验次数，尽快找到生产和科学实验中最优方案的一种科学试验的方法。例如：在现代体育实践的科学实验中，怎样选取最合适的配方、配比；寻找最好的操作和工艺条件；找出产品的最合理的设计参数，使产品的质量最好，产量最多，或在一定条件下使成本最低，消耗原料最少，生产周期最短等。把这种最合适、最好、最合理的方案，一般总称为最优；把选取最合适的配方、配比，寻找最好的操作和工艺条件，给出产品最合理的设计参数，叫做优选。也就是根据问题的性质在一定条件下选取最优方案。最简单的最优化问题是极值问题，这样问题用微分学的知识即可解决。
　　实际工作中的优选问题 ，即最优化问题，大体上有两类：

           一类是求函数的极值；如果目标函数有明显的表达式，一般可用微分法、变分法、极大值原理或动态规划等分析方法求解（间接选优）。

           另一类是求泛函的极值。如果目标函数的表达式过于复杂或根本没有明显的表达式，则可用数值方法或试验最优化等直接方法求解（直接选优）。
　
        优选法的优点：用较少的试验次数和外在条件，缩短时间和空间的消费、提高质量，达到解决问题的最经济的方法。

优选法基本步骤
　　1）选定优化判据（试验指标），确定影响因素，优选数据是用来判断优选程度的依据。
　　2）优化判据与影响因素直接的关系称为目标函数。
　　3）优化计算。优化(选）试验方法一般分为两类： 

                         分析法：同步试验法 ；黑箱法：循序试验法

优选法的分类
　　优选法分为单因素方法和多因素方法两类。单因素方法有平分法、0.618法(黄金分割法)、分数法、分批试验法等；多因素方法很多．但在理论上都不完备．主要有降维法、爬山法、单纯形调优胜。随机试验法、试验设计法等。优选法已在体育领域得到广泛应用。
　　1.单因素优选法
　　 如果在试验时，只考虑一个对目标影响最大的因素，其它因素尽量保持不变，则称为单因素问题。这个概念太复杂  了，通俗一点，把解决问题的结果看做方案的函数，当可供选择的方案分布在一条直线上时，其结果可看作一个单变量函数，此时的最优化问题为单因素优化问题。一般步骤：
　　（1）首先应估计包含最优点的试验范围,如果用a表示下限，b表示上限，试验范围为[a，b];
　　（2）然后将试验结果和因素取值的关系写成数学表达式,不能写出表达式时，就要确定评定结果好坏的方法。

一般可采用对分法、黄金分割法及牛顿一元函数求根的切线法与割线法等方法解决．也可用最小二乘法解决。

 

当多种方案是在一个平面上更高维数的空间中连续分布时，叫做二因素或多因索优化问题，此时可采用牛顿方法或梯度方法求解。若选择方案的分布是离散的，则可采用单纯形等方法解决。

      2.多因素优选法
　　多因素问题：首先对各个因素进行分析，找出主要因素，略去次要因素，划“多”为“少”，以利于解决问题。

 

 

 

黄金分割法（又称0.618法）是用来求单峰函数的max（或min）的算法。这是一种搜索法，不需要利用函数的导数值。

所用到的 0.618 是黄金分割比的近似值。
黄金分割比  = (sqrt(5)-1)/2 ~=  0.61803398874989484820... 又等于 斐波那契数列的 a(n)/a(n+1), n->∞

下边是我检索到的一些范例程序，没有经过调试，有兴趣自己研究一下。

Const k2 = 0.618033988749895

 Const k1 = 1 - k2 

Function M618(ByVal a#, ByVal b#, r#, L%, Nr%) As Double ' [a,b]-区间，r-精度，L=1求峰，L=0求谷

 Dim u#, v#, tu#, tv# 

    If a > b Then 

      Text2.Text = a: 

      Text1.Text = b 

    a = Val(Text1.Text):

    b = Val(Text2.Text) 

End If '开始计算 

    u = a + k1 * (b - a) 

    v = a + b - u

    tu = ff(u)

    tv = ff(v)

    Nr = 0 ' 

Do Until Abs(u - v) < r Or Nr > 70

 If tu < tv Xor L = 0 Then '判断函数值的大小 

   a = u u = v: 

   tu = tv 'u 点不再计算，要求K2要严格等于黄金分割比 

    v = a + k2 * (b - a) '求出新的v

    tv = ff(v) Else b = v v = u: 

    tv = tu u = a + k1 * (b - a) 

    tu = ff(u) End If 

    Nr = Nr + 1 Loop 

    M618 = (u + v) / 2

 End Function

 

 

黄金分割法求解 一元二次函数f(x)=x2+2x+1

#include "math.h" 

#include "stdio.h" 

#define f(x) x*x+2*x+1 //函数功能是用黄金分割法实现求一元函数 的最优解 

double hj(double *a,double *b,double e,int *n) 

{ 

double x1,x2,s;

 if(fabs(*b-*a)<=e) s=f((*b+*a)/2);

 else { x1=*a+0.382*(*b-*a); x2=*a+0.618*(*b-*a);

 if(f(x1)>f(x2)) *a=x1; else *b=x2; 

*n=*n+1; s=hj(a,b,e,n); 

} 

return s; 

} 

main() {

 double s,a,b,e; int n=0; 

 scanf("%lf %lf %lf",&a,&b,&e); // 输入区间[a,b]和精度e的值 

 s=hj(&a,&b,e,&n); //调用hj函数，其中n代表迭代次数 

 printf("a=%lf,b=%lf,s=%lf,n=%d\n",a,b,s,n); 

} 

运行时：输入：0.6 0.5 0.1 输出结果为： 0.6 0.5 0.1 a=0.600000,b=0.500000,s=2.402500,n=0

 

 

 

 

0.618算法优化子程序，程序包含4个C文件mhjfgf.c、funct.c、jtf.c、hjfgf.c

//mhjfgf.c代码如下：
#include "hjfgf.c"
main()
{double xx[1],a[1],b[1],ff;
double s[]={1},x0[]={0};
jtf(x0,0.1,s,1,a,b);
ff=gold(a,b,0.00001,1,xx);
printf("\nx[0]=%f,,ff=%f",xx[0],ff);
getch();
}


//funct.c代码如下：
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
double objf(double x[])
{double ff;
ff=8*pow(x[0],3)-2*x[0]*x[0]-7*x[0]+3;
return(ff);
}


//jtf.c代码如下：
#include "funct.c"
void jtf(double x0[],double h0,double s[],int n,double a[],double b[])
{int i;
double *x[3],h,f1,f2,f3;
for(i=0;i<3;i++)
  x[i]=(double *)malloc(n*sizeof(double));
  h=h0;
for(i=0;i<n;i++)
  *(x[0]+i)=x0[i];
f1=objf(x[0]);
for(i=0;i<n;i++)
  *(x[1]+i)=*(x[0]+i)+h*s[i];
  f2=objf(x[1]);
if(f2>=f1)
  {h=-h0;
    for(i=0;i<n;i++)
    *(x[2]+i)=*(x[0]+i);
   f3=f1;
    for(i=0;i<n;i++)
    {*(x[0]+i)=*(x[1]+i);
     *(x[1]+i)=*(x[2]+i);
    }
   f1=f2;
   f2=f3;
   }
   for(;;)
   {h=2*h;
     for(i=0;i<n;i++)
     *(x[2]+i)=*(x[1]+i)+h*s[i];
   f3=objf(x[2]);
   if(f2<f3) break;
   else
    { for(i=0;i<n;i++)
       {*(x[0]+i)=*(x[1]+i);
        *(x[1]+i)=*(x[2]+i);
       }
      f1=f2;
      f2=f3;
    }
   }
   if(h<0)
    for(i=0;i<n;i++)
    {a[i]=*(x[2]+i);
     b[i]=*(x[0]+i);
    }
   else
    for(i=0;i<n;i++)
    {a[i]=*(x[0]+i);
     b[i]=*(x[2]+i);
     }
   for(i=0;i<3;i++)
   free(x[i]);
}


//hjfgf.c代码如下：
#include "jtf.c"
double gold(double a[],double b[],double eps,int n,double xx[])
{int i;
double f1,f2,*x[2],ff,q,w;
for(i=0;i<2;i++)
  x[i]=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
  {*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
   *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
  }
  f1=objf(x[0]);
  f2=objf(x[1]);
  do
   {if(f1>f2)
     {for(i=0;i<n;i++)
      {b[i]=*(x[0]+i);
       *(x[0]+i)=*(x[1]+i);
       }
     f1=f2;
     for(i=0;i<n;i++)
      *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
     f2=objf(x[1]);
     }
    else
     { for(i=0;i<n;i++)
       {a[i]=*(x[1]+i);
       *(x[1]+i)=*(x[0]+i);}
     f2=f1;
    for(i=0;i<n;i++)
     *(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
    f1=objf(x[0]);
     }
  q=0;
  for(i=0;i<n;i++)
   q=q+(b[i]-a[i])*(b[i]-a[i]);
  w=sqrt(q);
  }while(w>eps);
  for(i=0;i<n;i++)
   xx[i]=0.5*(a[i]+b[i]);
  ff=objf(xx);
  for(i=0;i<2;i++)
  free(x[i]);
  return(ff);
}

 

 

前面提到过黄金分割与斐波那契（Fibonacci）数列的相同之处，现在对斐波那契作一些了解，以及世界十大难题中其一“兔子问题”。

★ 扩展(摘录):
    斐波那契是意大利的数学家。他是一个商人的儿子。儿童时代跟随父亲到了阿尔及利亚，在那里学到了许多阿拉伯的算术和代数知识，从而对数学产生了浓厚的兴趣。
　　长大以后，因为商业贸易关系，他走遍了许多国家，到过埃及、叙利亚、希腊、西西里和法兰西。每到一处他都留心搜集数学知识。回国后，他把搜集到的算术和代数材料，进行研究、整理，编写成一本书，取名为《算盘之书》，于1202年正式出版。
　　这本书是欧洲人从亚洲学来的算术和代数知识的整理和总结，它推动了欧洲数学的发展。其中有一道“兔子数目”的问题是这样的：一个人到集市上买了一对小兔子，一个月后，这对小兔子长成一对大兔子。然后这对大兔子每过一个月就可以生一对小兔子，而每对小兔子也都是经过一个月可以长成大兔子，长成大兔后也是每经过一个月就可以生一对小兔子。那么，从此人在市场上买回那对小兔子算起，每个月后，他拥有多少对小兔子和多少对大兔子？
　　这是一个有趣的问题。当你将小兔子和大兔子的对数算出以后，你将发现这是一个很有规律的数列，而且这个数列与一些自然现象有关。人们为了纪念这位兔子问题的创始人，就把这个数列称为“斐波那契数列”。 又找到了这么一段话:

规律表:

月数 小兔 中兔 老兔 总数
 1    1    0    0    1
 2    0    1    0    1
 3    1    0    1    2
 4    1    1    1    3
 5    2    1    2    5
 6    3    2    3    8
 7    5    3    5   13

    在计算每一行时，大兔数为上月的大兔数加上月的中兔数，中兔数为上月的小兔数，小兔数为本月的大兔数，算总数为本月的小兔数加本月的中兔数加本月的大兔数。在观察总数的过程中找出了规律：总数的第一、二月都是1，以后的每一月是前两月的和。数列为1，1，2，3，5，8，13，21，34，55，……

    当n=50时，后项与前项的比是1.61803398874989，而前项与后项的比是0.61803398874989，即b/a的值与a/b的值相差1，假设后项与前项的比是φ，则有（φ-1）/φ=1，解这个方程得：φ= (√5+1) /2，这就是黄金分割。
    当n充分大时，斐波纳契数列后前项的比值，与前后项的比值，相差1，它们的比值是黄金分割！黄金分割是一个十分有用的无理数。据此，把黄金分割可用一个有理数近似表示，如斐波纳契数列的第七项与斐波纳契数列的第六项的比13/8，斐波纳契数列的第九项与斐波纳契数列的第八项的比34/21等都可以近似地表示为黄金分割，当然项数越后越精确。

 

#include <stdio.h>
 int main()
 { 
long double f1, f2; 
int i;
 f1 = 1;
 f2 = 1; 
for (i = 1; i <= 25; i++) 
{
 f1 = f1 + f2; 
f2 = f1 + f2; 
} 
printf("%lf\n", (long double)f2/f1); return 0; } 

运行结果:1.618034

 

 ★ 有人验证, 如果使用该程序循环7次,即: i<=7,则结果为1.618033,循环8次以上(有效范围内),均为1.618034, 由于"对于双精度数，使用%lf格式符输出时，前16位是有效数字，小数6位",如何显示更多的小数位不妨自己测试一下，一下午都在了解黄金分割与最优化法的内容……呵呵

 

解决这样一个问题：

      数据规律F1=1    (n=1)；F2=1   (n=2)；Fn=Fn-1+Fn-2    (n>=3)

      显示Fibonacci数列的前n个数。

 

int main()
 {
 long int f1, f2; /* 作为循环体的"前两项",会超过32767,所以用long */ 
int i; /* 计数 */
 f1 = 1;
 f2 = 1; 
for (i = 1; i <= n/2; i++) /* 循环一次输出2个,共进行n/2次 */ 
{ 
printf("%12ld %12ld", f1, f2); /* 第22个数之后值超过32767，必须用ld,而非d */ 
if ( i % 2 == 0) /* 每次输出两个, 两次输出4个, 便换行 */
 printf("\n"); /* 控制换行 */ 
f1 = f1 + f2; /* 计算下两个数 */
 f2 = f1 + f2;
 } 
return 0; 
} 若不考虑其他一切条件的话（没人傻得去测试一个循环n的数列）输出结果: ================================================   1         1             2          3             5               8           13
 21        34          55          89          144          233        377
 610      987       1597      2584       4181        6765      10946
 17711  28657   46368    75025    121393    196418   317811
 ………… ================================================ 另一种数组实现的方法，如果你有还好的方法何不开源一下，呵呵。
 #include <stdio.h>
 int main() {
 int i; 
int f[20] = {1,1}; 
for (i = 2; i < 20; i++) 
f[i] = f[i-2] + f[i-1];
 for (i = 1; i <= 20; i++) 
{ 
printf("%8d", f[i-1]); 
if (i % 5 == 0)
 printf("\n"); 
} 
return 0; 
}

 

 

黄金分割在搜索区域中缩短探索区间的思想能够优化响应时间的问题，你可以拿一个具体的例子然后自己分析，如比较黄金分割查找与二分查找（折半）的时间复杂度和空间复杂度……

引用一个例子：

 

#include<iostream.h> 
#include<math.h> /*黄金分割法求最小值的C++程序，部分变量及函数书写并不规范*/ 
int n = (lnδ/ln0.618 + 1) + 1;//δ为题给精度
int i;
float f(float ai, float bi) 
{ 
a(i + 1) = ai + 0.618(bi - ai); 
return ai + 1; 
}
 float g(float ai, float bi)
 {
 b(i + 1) = ai + 0.382(bi - ai);
 return b(i + 1); 
} 
float F(float ai, float bi) 
{ 
//题给的f(x)函数式; 
result; 
} 
float A(float ai, float bi) 
{ 
    int i = 1; 
    float result;
 L: do { 
   a(i + 1) = f(float ai, float bi);
   b(i + 1) = g(float ai, float bi); 
   float F1 = F(float ai, float bi); 
   float F2 = F(float a(i + 1), 
   float b(i + 1)); 
   ai = ai, bi = b(i + 1);
    i ++; 
}
while(i <= n && F1 >= F2)
 if(i < n){ 
  B(float ai, float bi); 
  } 
else result = F2;
   return result;
   } 
float B(float ai, float bi) 
{ 
   do{
         a(i + 1) = f(float ai, float bi); 
         b(i + 1) = g(float ai, float bi); 
         float F1 = F(float ai, float bi); 
         float F2 = F(float a(i + 1), float b(i + 1));
         ai = a(i + 1), bi = bi;
          i ++; 
        }
   while(i <= n && F1 <= F2)
    if(i < n) 
       { goto L; } 
    else 
    result = F1; 
return result; 
}
 void main() { 
  int i = 1;
  float A(float ai, float bi);
  cout<<"最小值为："<<result<<endl; 
} 

 

 

先了解到这里 ……