优化方法在油藏生产优化^[1]、驱油试验的参数优化^[2]以及天然气水合物藏的注热水开采试验的参数研究^[3]中发挥着重要作用。一维优化方法作为多维优化方法的核心算法, 具有较高的研究价值, 其基本假设是目标函数在寻优方向上是单峰函数。一维优化方法可以分为试探类方法、区间削去类方法、函数拟合类方法和导数辅助定点类方法。一维二阶近似式定点法(牛顿法)属于第四类方法, 其数学推导严密, 理论价值大。将其推广应用到多维优化问题的求解, 则得到多维二阶近似式定点法(多维牛顿法^[4])。多维二阶近似式定点法在通用Gibbs反应器的机理建模和求解^[5]、图像处理的压缩感知重建算法^[6]中也起到重要作用。盲人探路算法^[7-8]具有很好的寻优效果, 将其用于求解多维有约束优化问题的随机方向法改进, 取得了较好的效果^[9]。笔者将该寻优思想用于多维二阶近似式定点法的改进, 研究其算法步骤及算例分析。

1 一维二阶近似式定点法

一维二阶近似式定点法也称为牛顿法。通常, 根据无穷阶近似级数建立目标函数f(x)的近似表达式, 然后将其极值点作为新的探测点。取连续可微目标函数在当前点x⁽⁰⁾附近的二阶近似多项式作为该函数的近似表达式:

$ \varphi \left( x \right) = f\left( {{x^{\left( 0 \right)}}} \right) + f'\left( {{x^{\left( 0 \right)}}} \right)\left( {x - {x^{\left( 0 \right)}}} \right) + \frac{1}{{2!}}f''\left( {{x^{\left( 0 \right)}}} \right){\left( {x - {x^{\left( 0 \right)}}} \right)^2}. $

(1)

根据极值条件, 极值点满足

$ f'\left( {{x^{\left( 0 \right)}}} \right) + f''\left( {{x^{\left( 0 \right)}}} \right)\left( {x - {x^{\left( 0 \right)}}} \right) = 0, $

(2)

即为直线的零点。该直线为目标函数的导函数f'(x)在当前点x⁽⁰⁾的切线(图 1)。该极值点的计算式为

$ x = {x^{\left( 0 \right)}} - \frac{{f'\left( {{x^{\left( 0 \right)}}} \right)}}{{f''\left( {{x^{\left( 0 \right)}}} \right)}}. $

(3)

在当前点x^(k)处, 由式(3)获得新的探测点为

$ {x^{\left( {k + 1} \right)}} = {x^{\left( k \right)}} - \frac{{f'\left( {{x^{\left( k \right)}}} \right)}}{{f''\left( {{x^{\left( k \right)}}} \right)}}. $

(4)

相当于在目标函数导函数曲线上做切线, 以该切线的零点作为新的探测点。其寻优过程如图 1所示。该方法也称为导函数切线零点法。

对于多峰函数, 如果探测点到达极大点或拐点, 则由于目标函数一阶导数为零或二阶导数为零而导致寻优失败。因此, 须取迭代点为

$ {x^{\left( {k + 1} \right)}} = {x^{\left( k \right)}} + \Delta x. $

(5)

式中, Δx为给定小值。其收敛精度准则宜选用梯度准则和点距准则。

对于复杂的目标函数, 如果初始点离极值点太远, 采用一维二阶近似式定点法易发生迭代点在极值点两边交替出现的现象, 从而导致寻优失败。为增强收敛稳定性, 引入阻尼因子(0 < α < 1), 减小新点相对于当前点的变化值, 得到一维二阶近似式阻尼定点法, 其迭代公式为

$ {x^{\left( {k + 1} \right)}} = \alpha {x^{\left( k \right)}} + \left( {1 - \alpha } \right)\left( {{x^{\left( k \right)}} - f'\left( {{x^{\left( k \right)}}} \right)/f''\left( {{x^{\left( k \right)}}} \right)} \right). $

(6)

上述阻尼法的步长具有一定的模式, 其适用性必然受到限制。如果基于盲人探路思想在探测点劣于当前点时将步长逐步减半试探, 则可增强算法的稳定性。

2 多维二阶近似式定点法的迭代公式

与一维二阶近似式定点法类似, 在当前点x^(k)处取目标函数的二阶近似多项式作为目标函数的近似表达式, 表示为

$ \begin{array}{l} f\left( \mathit{\boldsymbol{x}} \right) \approx J\left( \mathit{\boldsymbol{x}} \right) = f\left( {{\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right) + {\nabla ^{\rm{T}}}f\left( {{\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right)\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right) + \\ \frac{1}{2}{\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right)^{\rm{T}}}{\nabla ^{\rm{2}}}f\left( {{\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right)\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right). \end{array} $

(7)

其中, 二阶偏导数矩阵(海森矩阵)∇²f(x^(k))为

$ \mathit{\boldsymbol{H}}\left( f \right) = {\left[{\begin{array}{*{20}{c}} {\frac{{{\partial ^2}f}}{{\partial x_1^2}}}&{\frac{{{\partial ^2}f}}{{\partial {x_1}\partial {x_2}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {x_1}\partial {x_N}}}}\\ {\frac{{{\partial ^2}f}}{{\partial {x_2}\partial {x_1}}}}&{\frac{{{\partial ^2}f}}{{\partial x_2^2}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {x_2}\partial {x_N}}}}\\ \vdots&\vdots &{}& \vdots \\ {\frac{{{\partial ^2}f}}{{\partial {x_N}\partial {x_1}}}}&{\frac{{{\partial ^2}f}}{{\partial {x_N}\partial {x_2}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial x_N^2}}} \end{array}} \right]_{x = {x^{\left( k \right)}}}}. $

式(7)的极值点满足:

$ \nabla f\left( {{\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right) + \mathit{\boldsymbol{H}}\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right) = 0. $

(8)

令该极值点为下一个迭代点x^(k+1), 得

$ {\mathit{\boldsymbol{x}}^{\left( {k + 1} \right)}} = {\mathit{\boldsymbol{x}}^{\left( k \right)}} + \left[{-{\mathit{\boldsymbol{H}}^{-1}}\nabla f\left( {{\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right)} \right]. $

(9)

式(9)即多维二阶近似式定点法的迭代公式。

对于非二次函数, 有时迭代点目标函数值会上升, 从而丧失收敛稳定性, 也可引入阻尼因子(0 < α < 1)得到多维二阶近似式阻尼定点法, 其迭代公式为

$ {\mathit{\boldsymbol{x}}^{\left( {k + 1} \right)}} = {\mathit{\boldsymbol{x}}^{\left( k \right)}} - \alpha {\mathit{\boldsymbol{H}}^{ - 1}}\nabla f\left( {{\mathit{\boldsymbol{x}}^{\left( k \right)}}} \right). $

(10)

如果将迭代点修正量作为寻优方向进行一维寻优, 该方向上的最优点作为新的迭代点x^(k+1), 则称为多维二阶近似式方向法。其寻优效果较好^[10]。

多维二阶近似式定点法及其阻尼法的步长具有固定的模式, 与一维方法相同, 基于盲人探路思想在探测点劣于当前点时将步长逐步减半试探, 也可增强算法的稳定性。

3 基于盲人探路寻优思想的改进 3.1 寻优思想流程

只有二阶近似式极值点优于当前点时才更新当前点, 否则退半步探测, 如果探测点优于当前点, 则将其置为当前点。基于盲人探路寻优思想的二阶近似式定点法的寻优思想流程图见图 2。

3.2 算法步骤

(1) 给定初始点, 并设置为当前点; 给定终止条件。

(2) 根据式(9)计算当前点处二阶近似式的极值点, 设该极值点与当前点之差为Δ; 令寻优步长α为1。

(3) 计算探测点x⁽¹⁾=x⁽⁰⁾+αΔ。

(4) 如果x⁽¹⁾劣于x⁽⁰⁾, 根据一维盲人探路寻优思想交替令α=α/2或Δ=-Δ, 转(3);若α足够小, 则结束。

(5) 令探测点x⁽¹⁾为当前点x⁽⁰⁾。

(6) 如果步长曾减半, 转(8)。

(7) 按步长Δ一步步探测, 若探测点优于当前点则将其置为当前点, 直到探测点劣于当前点为止。

(8) 如果当前点满足终止条件(当前点梯度的模足够小或偏导数不存在), 则转(10)。

(9) 如果二阶偏导数矩阵奇异, 则转(11), 否则转(2)。

(10)若当前点与极小值足够接近, 则结束寻优。

(11)参照式(5), 当前点平移一步, 转(2)。

针对具体的优化问题, 有的步骤是多余的, 徒增计算量, 可舍去。不基于盲人探路寻优思想的算法为步骤(1)、(2)、(3)、(8)~(11)。

3.3 计算子程序部分代码

int Second_appro(x0, y0, e, gr);

  /*二阶近似式定点法子程序*/

double x0[2], y0[1], e, gr[2];/*初始点, 初始解(返回最优点), 终止条件, 梯度*/

    {double x1[2], y1, H_ni[2][2], dx[2];

/*试探点, 二阶偏导数矩阵的逆阵, 单步修正量*/

  double h; /*探测步长*/

  int i, j, k=0;/*当前点更新的次数*/

  int n=2;/*以二维为例*/

  y0[0]=func(x0[0], x0[1]);

    /*当前点的目标函数值*/

  for(; ;)/*计算循环*/

  {

  gradu(x0[0], x0[1], gr); /*目标函数的梯度*/

  Hessain(x0[0], x0[1], H_ni);

  /*目标函数的二阶偏导数矩阵及其逆阵*/

dx[0]=-(H_ni[0][0]*gr[0]+H_ni[0][1]*gr[1]);

    /*单步修正量*/

dx[1]=-(H_ni[1][0]*gr[0]+H_ni[1][1]*gr[1]);

  h=1.0; i=1;/*步长, 方向标志*/

  for(; ;)/*减半步长循环:获得优于当前点的探测点*/

{for(j=0;j < n; j++)

  x1[j]=x0[j]+h*dx[j];

  y1=func(x1[0], x1[1]);

    /*探测点及其目标函数值*/

  if(y0[0]>y1)

   break;

    /*如果探测点比当前点好*/

  if(i==-1)

    {h=0.5*h;

    i=-i;

    }

   else

    {i=-i;

for(j=0;j < n; j++)

    dx[j]=-dx[j];

    }

  }

  for(i=0;i < n; i++)

    x0[i]=x1[i];

  y0[0]=y1;/*步骤(5)*/

  if(h>0.9)/*如果没有步长减半成功的操作, 则沿原来的方向一步步走下去*/

   {for(; ;)/*步骤(7)*/

  {for(i=0;i < n; i++)

     x1[i]=x0[i]+h*dx[i];

     y1=func(x1[0], x1[1]);

       /*探测点及其目标函数值*/

     if(y0[0]>y1)

  {for(i=0;i < n; i++)

     x0[i]=x1[i]; y0[0]=y1;

     }

     else break; /*如果探测点比当前点好*/

    }

   }

   k=k+1;

   gradu(x0[0], x0[1], gr); /*目标函数的梯度*/

   if((gr[0]*gr[0]+gr[1]*gr[1]) < e*e)

  {if(zhengding(x0[0], x0[1])==1)

  return; }

  /*二阶偏导数矩阵正定, 则取当前点为最优点, 结束寻优*/

   else

  {if(zhengding(x0[0], x0[1])!=0)

  continue; }

  /*如果二阶偏导数矩阵不奇异, 则下一次迭代*/

  /*当前点平移一小步, 百分之一当前步长*/

  for(i=0;i < n; i++)

    x1[i]=x0[i]+0.01*h*dx[i];

  if(zhengding(x1[0], x1[1])==0)

      /*如果仍然奇异, 则沿垂直方向平移。*/

  {x1[0]=x0[0]+0.01*h*dx[1];

  x1[1]=x0[1]-0.01*h*dx[0];

      }

    for(i=0;i < n; i++)

    x0[i]=x1[i];

    y0[0]=func(x1[0], x1[1]);

  }/*大的计算循环结束*/

  }

4 算例分析 4.1 算例一

$ \min f\left( \mathit{\boldsymbol{x}} \right) = \left\{ \begin{array}{l} \sin \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right], \varphi \left( {{x_1}, {x_2}} \right) < 1.5;\\ \sin 1.5, \varphi \left( {{x_1}, {x_2}} \right) \ge 1.5. \end{array} \right. $

其中

$ \varphi \left( {{x_1}, {x_2}} \right) = 0.01{\left( {{x_1} + {x_2} - 5} \right)^2} + 0.001{\left( {{x_1} - {x_2}} \right)^2}. $

在φ(x₁, x₂)∈[0, 1.5)内, 点x处的梯度∇f$\left( \mathit{\boldsymbol{x}} \right) = {\left[{\frac{{\partial f\left( \mathit{\boldsymbol{x}} \right)}}{{\partial {x_1}}}\;\;\frac{{\partial f\left( \mathit{\boldsymbol{x}} \right)}}{{\partial {x_2}}}} \right]^{\rm{T}}} $为

$ \left[\begin{array}{l} \cos \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right]\left( {0.022{x_1} + 0.018{x_2} - 0.1} \right)\\ \cos \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right]\left( {0.018{x_1} + 0.022{x_2} -0.1} \right) \end{array} \right]. $

二阶偏导数矩阵(海森矩阵)∇²f(x)为

$ \mathit{\boldsymbol{H = }}\left[{\begin{array}{*{20}{c}} A&B\\ B&C \end{array}} \right], {\mathit{\boldsymbol{H}}^{ - 1}} = \frac{1}{{AC - {B^2}}}\left[{\begin{array}{*{20}{c}} C&{-B}\\ {-B}&A \end{array}} \right]. $

(11)

其中

$ \begin{array}{l} A = - \sin \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right]{\left( {0.022{x_1} + 0.018{x_2} - 0.1} \right)^2} + \\ 0.022\cos \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right], \end{array} $

$ \begin{array}{l} B = - \sin \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right]\left( {0.018{x_1} + 0.022{x_2} - 0.1} \right) \times \\ \left( {0.022{x_1} + 0.018{x_2} - 0.1} \right) + 0.018\cos \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right], \end{array} $

$ \begin{array}{l} C = - \sin \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right]{\left( {0.018{x_1} + 0.022{x_2} - 0.1} \right)^2} + \\ 0.022\cos \left[{\varphi \left( {{x_1}, {x_2}} \right)} \right]. \end{array} $

不采用步骤(7)的算法, 如果目标函数不采用分段形式, 则寻优结果如表 1所示。

如果目标函数采用分段形式, 则寻优结果如表 2所示。

4.2 算例二

$ \min f\left( \mathit{\boldsymbol{x}} \right) = \sqrt[8]{{{\varphi _2}\left( {{x_1}, {x_2}} \right)}}. $

其中

$ {\varphi _2}\left( {{x_1}, {x_2}} \right) = 10{\left( {{x_1} + {x_2} - 5} \right)^2} + {\left( {{x_1} - {x_2}} \right)^2}. $

在x处的梯度∇f(x)为

$ \left[{\begin{array}{*{20}{c}} {\frac{1}{8}{{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]}^{ - \frac{7}{8}}}\left( {22{x_1} + 18{x_2} - 100} \right)}\\ {\frac{1}{8}{{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]}^{ -\frac{7}{8}}}\left( {18{x_1} + 22{x_2} -100} \right)} \end{array}} \right]. $

二阶偏导数矩阵(海森矩阵)∇²f(x)为

$ \mathit{\boldsymbol{H = }}\left[{\begin{array}{*{20}{c}} A&B\\ B&C \end{array}} \right]. $

其中

$ \begin{array}{l} A = - \frac{7}{{64}}{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^{ - \frac{{15}}{8}}}{\left( {22{x_1} + 18{x_2} - 100} \right)^2} + \\ \frac{{22}}{8}{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^{ - \frac{7}{8}}}, \end{array} $

$ \begin{array}{l} B = - \frac{7}{{64}}{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^{ - \frac{{15}}{8}}}\left( {22{x_1} + 18{x_2} - 100} \right)\left( {18{x_1} + } \right.\\ \left. {22{x_2} - 100} \right) + \frac{{18}}{8}{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^{ - \frac{7}{8}}}, \end{array} $

$ \begin{array}{l} C = - \frac{7}{{64}}{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^{ - \frac{{15}}{8}}}{\left( {18{x_1} + 22{x_2} - 100} \right)^2} + \\ \frac{{22}}{8}{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^{ - \frac{7}{8}}}. \end{array} $

不采用步骤(7)的算法, 寻优过程如图 3所示。

迭代点距离极值点较远, 但是均指向极值点。在寻优过程中须反向才能靠近极值点。

4.3 算例三

$ \min f\left( \mathit{\boldsymbol{x}} \right) = {\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^4}. $

在x处的梯度∇f(x)为

$ \left[\begin{array}{l} 4{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^3}\left( {22{x_1} + 18{x_2} - 100} \right)\\ 4{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^3}\left( {18{x_1} + 22{x_2} -100} \right) \end{array} \right]. $

二阶偏导数矩阵(海森矩阵)∇²f(x)为

$ \mathit{\boldsymbol{H = }}\left[{\begin{array}{*{20}{c}} A&B\\ B&C \end{array}} \right]. $

其中

$ \begin{array}{l} A = 12{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^2}{\left( {22{x_1} + 18{x_2} - 100} \right)^2} + \\ 88{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^3}, \end{array} $

$ \begin{array}{l} B = 12{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^2}\left( {22{x_1} + 18{x_2} - 100} \right)\left( {18{x_1} + 22{x_2} - } \right.\\ \left. {100} \right) + 72{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^3}, \end{array} $

$ \begin{array}{l} C = 12{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^2}{\left( {18{x_1} + 22{x_2} - 100} \right)^2} + \\ 88{\left[{{\varphi _2}\left( {{x_1}, {x_2}} \right)} \right]^3}. \end{array} $

采用不基于盲人探路寻优思想的多维二阶近似式定点法, 经过30次寻优方向(或步长)的确定和31次目标函数的调用, 获得最优点为[2.377 42.553 9]^T的寻优结果, 寻优过程如图 4所示。

采用基于盲人探路寻优思想的多维二阶近似式定点法, 经过1次寻优方向的确定, 8次目标函数调用, 获得最优点为[2.5 2.5]^T的寻优结果, 寻优过程如图 5所示。当前点的步长为与极值点距离的七分之一。

这3个算例的目标函数等值线均为椭圆族, 因此寻优方向无异议。如果目标函数为

$ f\left( \mathit{\boldsymbol{x}} \right) = x_1^4 + x_1^2 + 3x_2^2. $

(12)

初始点不变, 基于盲人探路思想的寻优过程如图 6所示。每条寻优方向上只计算一次目标函数的梯度、海森矩阵及其逆阵, 第一个点为不采用本文算法的迭代点, 可见本文算法的计算量小。

5 结束语

根据算法的特点以“多维二阶近似式定点法”更名“牛顿型多维优化方法”, 给出的一维和多维迭代公式可用于相关优化问题的求解。提出的基于盲人探路寻优思想的优化算法具有实用性强、计算量小的特点。对于算例一和算例二的部分情况, 按照式(9)所确定的步长反向才能获得靠近极值点的迭代点, 这是因为目标函数的最高幂次小于1, 导致牛顿方向不是下降方向。算例二和算例三的极值点在当前点指向迭代点的方向上或其相反方向上。