数学建模之数学模型—2：非线性规划

非线性规划

基本概念与结论

有约束非线性规划的一般数学模型为如下形式：
$$\begin{gathered} \min f\left(x\right) \\ s.t.\{\begin{array}{c}{g}_i\left(x\right)⩾0,i=1,2,...,m\\ h_j\left(x\right)=0,j=1,2,...,n\end{array} \end{gathered}$$
下面给出非线性规划的基本概念
定义1：设$f(x)$为目标函数，$S$为可行域，$x^{\ast }\in S$。若：$\forall x\in S,f\left(x\right)⩾f\left(x^{\ast }\right)$则称$x^{\ast }$为
$$\begin{gathered} \min f\left(x\right) \\ s.t.x\in S \end{gathered}$$
极小化问题的最优解(全局最优解)。
定义2：设$f(x)$为目标函数，$S$为可行域，$x^{\ast }\in S$。若： ∃ x ∗ ∈ { x ∣ ∥ x − x ∗ ∥ < ε , ε > 0 } , f ( x ) ⩾ f ( x ∗ ) \exist x^*\in \left\{ x\left| \left\| x-x^* \right\| \right. <\varepsilon ,\varepsilon >0 \right\} ,f\left( x \right) \geqslant f\left( x^* \right) ∃x∗∈{x∣∥x−x∗∥<ε,ε>0},f(x)⩾f(x∗)则称$x^{\ast }$为极小化问题：
$$\begin{gathered} \min f\left(x\right) \\ s.t.x\in S \end{gathered}$$
的局部最优解。

凸集与凸函数

定义：设$f：S\subset R^n\to R$,$S$是凸集，如果对所有的$x^{(1)},x^{(2)}\in S$,$\lambda \in (0,1)$,都有：$$f(\lambda x^{(1)}+(1-\lambda )x^{(2)})⩽\lambda f(x^{(1)})+(1-\lambda )f(x^{(2)})$$
称$f(x)$为$S$上的凸函数。
Hessian矩阵：$f$存在二阶偏导数，$x\in R^n$:$${\nabla }^{2}f\left(x\right)=\left(\begin{array}{cccc}\frac{{\partial }^{2}f\left(x\right)}{\partial x_1^2}& \frac{{\partial }^{2}f\left(x\right)}{\partial x_1\partial x_2}& ...& \frac{{\partial }^{2}f\left(x\right)}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f\left(x\right)}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f\left(x\right)}{\partial x_2^2}& ...& \frac{{\partial }^{2}f\left(x\right)}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f\left(x\right)}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f\left(x\right)}{\partial x_n\partial x_2}& ...& \frac{{\partial }^{2}f\left(x\right)}{\partial x_n^2}\end{array}\right)$$
下面给出一系列定理：

• 定理1：有限凸函数的线性组合为凸函数
• 定理2：凸函数在有限函数值处的沿任意方向导数都存在
• 定理3：凸函数在定义区间上必连续
• 定理4：设$S$是$R^n$上的非空凸集，$f(x)$是定义在$S$上的凸函数，则$f(x)$在$S$上的局限极小点是全局极小点，且极小点的集合为凸集
• 凸函数的一阶、二阶判断定理：
  • $$\exists x^{\left(1\right)}\ne x^{\left(2\right)}\in S\ne \oslash ,S\text{是一个凸集，}f\text{可微}:f\text{为凸函数}\iff f\left(x^{\left(2\right)}\right)⩾f\left(x^{\left(1\right)}\right)+\nabla f{\left(x^{\left(1\right)}\right)}^T\left(x^{\left(2\right)}-x^{\left(1\right)}\right)$$
  • 在凸集定义内的点的Hessian矩阵半正定$\iff$凸函数
  • 在凸集定义上所有的点在其处的Hessian矩阵正定$\Rightarrow$严格凸函数
    凸规划
    $$\begin{gathered} \min f\left(x\right) \\ s.t.\{\begin{array}{c}{g}_i\left(x\right)⩾0,i=1,2,...,m\\ h_j\left(x\right)=0,j=1,2,...,n\end{array} \end{gathered}$$
    $f\left(x\right)$为凸函数，$g_i\left(x\right)$为凹函数，$h_j\left(x\right)$为线性函数，称此类问题为凸规划

极值条件

无约束条件的极值判断条件

无约束非线性规划的一般数学模型为如下形式：$$\begin{gathered} \min f\left(x\right) \\ x\in {\mathbb{R}}^n \end{gathered}$$
定理1：$$\text{设函数}f\left(x\right)\text{在}\overline{x}\text{处可微},\text{若存在向量}d,\text{使得}\nabla f{\left(\overline{x}\right)}^{T}d<0,\text{则存在数}\delta >0,\text{使得对}\forall \lambda \in \left(0,\delta \right),\text{有}f\left(\overline{x}+\lambda d\right)<f\left(\overline{x}\right)$$
定理2：$$\text{设函数}f\left(x\right)\text{在点}\overline{x}\text{处可微，若}\overline{x}\text{是局部最小点},\text{则}\nabla f\left(\overline{x}\right)=0$$
定理3：$$\text{设函数}f\left(x\right)\text{在点}\overline{x}\text{处二次可微},\text{若}\overline{x}\text{是局部极小点，则}\nabla f\left(\overline{x}\right)=0,\text{且}{\nabla }^{2}f\left(\overline{x}\right)\text{半正定}$$
定理4——局部极小点：$$\text{设函数}f\left(x\right)\text{在点}\overline{x}\text{处二次可微},\text{若}\nabla f\left(\overline{x}\right)=0,\text{且}{\nabla }^{2}f\left(\overline{x}\right)\text{正定，则}\overline{x}\text{是局部极小点}$$

定理5——全局最小点：$$\text{设函数}f\left(x\right)\text{在}{R}^n\text{上的可微函数},\overline{x}\in R^n,\text{则}\overline{x}\text{是全局极小点}\iff \nabla f\left(\overline{x}\right)=0$$

有约束条件的极值判断条件

有约束非线性规划的一般数学模型为如下表现形式：
$$\begin{gathered} \min f\left(x\right) \\ s.t.\{\begin{array}{c}{g}_i\left(x\right)⩾0,i=1,2,...,m\\ h_j\left(x\right)=0,j=1,2,...,n\end{array} \end{gathered}$$
其中$x\in R^n$，$f(x)$，$g_i(x)$，$h_j(x)$都是定义在$R^n$上的实值函数。
FRitz John条件——一般约束的一阶必要条件：
x ‾ 为可行点 , I = { i ∣ g i ( x ‾ ) = 0 } , f 和 g i ( x ‾ ) ( i ∈ I ) 在点 x ‾ 处可微 , g i ( x ‾ ) ( i ∉ I ) 在点 x ‾ 处连续， h j ( x ) ( j = 1 , 2 , . . . , l ) 在点 x ‾ 处连续可微。如果 x ‾ 是局部最优解 , 则存在不全为 0 的数 u 0 , u i ( i ∈ I ) 和 v j ( j = 1 , 2 , . . . , l ) , 使得： u 0 ∇ f ( x ‾ ) + ∑ i ∈ I u i ∇ g i ( x ∗ ) − ∑ j = 1 l μ j ∇ h j ( x ∗ ) = 0 , u 0 , u i ⩾ 0 \overline{x}\text{为可行点},I=\left\{ i|g_i\left( \overline{x} \right) =0 \right\} ,f\text{和}g_i\left( \overline{x} \right) \left( i\in I \right) \text{在点}\overline{x}\text{处可微},g_i\left( \overline{x} \right) \left( i\notin I \right) \text{在点}\overline{x}\text{处连续，}h_j\left( x \right) \left( j=1,2,...,l \right) \text{在点}\overline{x}\text{处连续可微。如果}\overline{x}\text{是局部最优解},\text{则存在不全为}0\text{的数}u_0,u_i\left( i\in I \right) \text{和}v_j\left( j=1,2,...,l \right) ,\text{使得：} \\ u_0\nabla f(\overline{x})+\sum_{i\in I}^{}{u_i\nabla g_i(x^*)}-\sum_{j=1}^l{\mu _j\nabla h_j(x^*)}=0, u_0,u_i\geqslant 0 x为可行点,I={i∣gi​(x)=0},f和gi​(x)(i∈I)在点x处可微,gi​(x)(i∈/I)在点x处连续，hj​(x)(j=1,2,...,l)在点x处连续可微。如果x是局部最优解,则存在不全为0的数u0​,ui​(i∈I)和vj​(j=1,2,...,l),使得：u0​∇f(x)+i∈I∑​ui​∇gi​(x∗)−j=1∑l​μj​∇hj​(x∗)=0,u0​,ui​⩾0

无约束非线性规划

一维搜索

一维搜索又称直线搜索，是指单变量函数最优化，即求一维问题的最优解方法，这里介绍黄金分割法

• 先在区间$\left[a,b\right]$上插入两个试探点：$t_l=a+\left(1-\beta \right)\left(b-a\right),t_r=a+\beta \left(b-a\right)$
• 由单峰函数性质丢掉$\frac{1}{3}$区间，逐步搜索
  黄金分割法是一维搜索中用于求解单变量函数最优解的一种高效方法，以下是对它更详细的介绍：

算法步骤

1. 初始化
   • 给定初始搜索区间$[a,b]$，以及精度要求$ϵ$。
2. 计算试探点
   • 计算两个试探点：$t_l=a+(1-\beta )(b-a)$，$t_r=a+\beta (b-a)$。
3. 比较函数值
   • 计算$f(t_l)$和$f(t_r)$，比较它们的大小。
   • 若$f(t_l)<f(t_r)$，则新的搜索区间为$[a,t_r]$。
   • 若$f(t_l)\ge f(t_r)$，则新的搜索区间为$[t_l,b]$。
4. 判断终止条件
   • 检查新的搜索区间长度是否小于精度要求$ϵ$，即$|b-a|<ϵ$。
   • 若是，则取当前区间的中点或两个试探点中函数值较好的点作为近似最优解，算法结束。
   • 若否，则返回步骤2，继续迭代。

示例

假设要求函数$f(x)=-x^2+4x-3$在区间$[0,4]$上的最大值。

1. 首先，初始化$a=0$，$b=4$，$\beta =0.618$。
2. 计算试探点：
   • $t_l=0+(1-0.618)\times (4-0)=1.528$。
   • $t_r=0+0.618\times (4-0)=2.472$。
3. 计算函数值：
   • $f(t_l)=f(1.528)=-1.528^2+4\times 1.528-3=1.315$。
   • $f(t_r)=f(2.472)=-2.472^2+4\times 2.472-3=1.985$。
4. 因为$f(t_l)<f(t_r)$，所以新的搜索区间为$[0,2.472]$。
5. 继续迭代，直到满足精度要求。

特点

• 优点
  • 不需要计算函数的导数，适用于导数不存在或难以计算的函数。
  • 算法简单，易于实现。
  • 收敛速度较快，每次迭代都能将搜索区间缩小一定比例。
• 缺点
  • 对于非单峰函数可能无法找到全局最优解。
  • 收敛速度相对较慢，与一些基于导数的方法相比，需要更多的迭代次数才能达到相同的精度。

代码模板

import math

# 定义黄金分割比
GOLDEN_RATIO = (math.sqrt(5) - 1) / 2

def golden_section_search(func, a, b, epsilon=1e-6):
    """
    黄金分割法进行一维搜索
    :param func: 目标函数
    :param a: 搜索区间左端点
    :param b: 搜索区间右端点
    :param epsilon: 精度要求
    :return: 近似最优解
    """
    # 计算初始试探点
    t_l = a + (1 - GOLDEN_RATIO) * (b - a)
    t_r = a + GOLDEN_RATIO * (b - a)
    # 计算初始试探点的函数值
    f_l = func(t_l)
    f_r = func(t_r)

    while (b - a) > epsilon:
        if f_l < f_r:
            # 新的搜索区间为 [a, t_r]
            b = t_r
            t_r = t_l
            f_r = f_l
            t_l = a + (1 - GOLDEN_RATIO) * (b - a)
            f_l = func(t_l)
        else:
            # 新的搜索区间为 [t_l, b]
            a = t_l
            t_l = t_r
            f_l = f_r
            t_r = a + GOLDEN_RATIO * (b - a)
            f_r = func(t_r)

    # 返回近似最优解
    return (a + b) / 2

# 示例目标函数
def example_function(x):
    return -x**2 + 4*x - 3

# 测试代码
if __name__ == "__main__":
    a = 0
    b = 4
    result = golden_section_search(example_function, a, b)
    print(f"近似最优解: {result}")
    print(f"最优解处的函数值: {example_function(result)}")

最速下降法

最速下降法属于迭代算法。求解无约束非线性规划的迭代算法一般具有如下形式

• 选定初始点$x^{(0)},k=0$
• $k=k+1$
• 寻找一个合适的方向$P^{(k)}$
• 求出沿着这个方向的前进步长${\lambda }^{(k)}$
• 得到新的点$x^{(k+1)}=x^{(k)}+{\lambda }^{(k)}{P}^{(k)}$
  最速下降算法的基本思想是沿着负梯度方向搜索极值点，即搜索方向：
  $$P^{(k)}=-\nabla f(x^{(k)})$$
  搜索步长满足以下条件:
  $$f(x^{(k)}+{\lambda }^{(k)}{P}^{(k)})=\underset{\lambda ⩾0}{\min }f(x^{(k)}-\lambda \nabla f(x^{(k)}))$$

算法详细步骤

1. 初始化
   选定初始点 $x^{(0)}$，并令迭代计数器 $k=0$。初始点的选择通常会影响算法的收敛速度和最终结果，在实际应用中，可以根据问题的特点或经验来选择合适的初始点。
2. 迭代过程

• 步骤一：更新迭代计数器
  • $k=k+1$
• 步骤二：确定搜索方向
  • 计算当前点 $x^{(k)}$ 处目标函数 $f(x)$ 的梯度 $\nabla f(x^{(k)})$，搜索方向 $P^{(k)}$ 取为负梯度方向，即 $P^{(k)}=-\nabla f(x^{(k)})$。梯度 $\nabla f(x)$ 是一个向量，它的每个分量是目标函数关于各个变量的偏导数。
• 步骤三：确定搜索步长
  • 搜索步长 ${\lambda }^{(k)}$ 要满足 $f(x^{(k)}+{\lambda }^{(k)}{P}^{(k)})=\underset{\lambda ⩾0}{\min }f(x^{(k)}-\lambda \nabla f(x^{(k)}))$。这意味着要在负梯度方向上找到一个最优的步长，使得目标函数值下降最多。在实际计算中，精确求解这个一维搜索问题可能比较困难，因此常常采用一些近似的方法，如黄金分割法、斐波那契法等进行一维搜索来确定步长。
• 步骤四：更新迭代点
  • 得到新的点 $x^{(k+1)}=x^{(k)}+{\lambda }^{(k)}{P}^{(k)}$。

3. 终止条件
   可以根据不同的准则来判断迭代是否终止，常见的终止条件有：

• 梯度的模小于某个给定的小正数 $ϵ$，即 $\Vert \nabla f(x^{(k)})\Vert <ϵ$。这表示目标函数在当前点的变化率已经很小，可能已经接近极值点。
• 相邻两次迭代点的距离小于某个给定的小正数 $ϵ$，即 $\Vert x^{(k+1)}-x^{(k)}\Vert <ϵ$。
• 目标函数值的变化小于某个给定的小正数 $ϵ$，即 $\left|f(x^{(k+1)})-f(x^{(k)})\right|<ϵ$。

代码实现示例

import numpy as np
import math

# 定义黄金分割比
GOLDEN_RATIO = (math.sqrt(5) - 1) / 2

# 黄金分割法进行一维搜索
def golden_section_search(func, a=0, b=10, epsilon=1e-6):
    t_l = a + (1 - GOLDEN_RATIO) * (b - a)
    t_r = a + GOLDEN_RATIO * (b - a)
    f_l = func(t_l)
    f_r = func(t_r)

    while (b - a) > epsilon:
        if f_l < f_r:
            b = t_r
            t_r = t_l
            f_r = f_l
            t_l = a + (1 - GOLDEN_RATIO) * (b - a)
            f_l = func(t_l)
        else:
            a = t_l
            t_l = t_r
            f_l = f_r
            t_r = a + GOLDEN_RATIO * (b - a)
            f_r = func(t_r)

    return (a + b) / 2

# 最速下降法
def steepest_descent(f, grad_f, x0, epsilon=1e-6, max_iter=1000):
    """
    最速下降法求解无约束非线性规划问题
    :param f: 目标函数
    :param grad_f: 目标函数的梯度函数
    :param x0: 初始点
    :param epsilon: 收敛精度
    :param max_iter: 最大迭代次数
    :return: 最优解
    """
    x = np.array(x0, dtype=np.float64)
    k = 0
    while k < max_iter:
        # 计算当前点的梯度
        grad = grad_f(x)
        # 判断是否满足收敛条件
        if np.linalg.norm(grad) < epsilon:
            break
        # 确定搜索方向，为负梯度方向
        p = -grad

        # 定义一维搜索的目标函数
        def one_dim_func(lambda_):
            return f(x + lambda_ * p)

        # 使用黄金分割法确定步长
        lambda_ = golden_section_search(one_dim_func)

        # 更新迭代点
        x = x + lambda_ * p
        k += 1
    return x

# 示例目标函数
def example_f(x):
    """
    示例目标函数：f(x) = x1^2 + x2^2
    """
    return np.sum(np.square(x))

# 示例目标函数的梯度
def example_grad_f(x):
    """
    示例目标函数的梯度：grad f(x) = [2*x1, 2*x2]
    """
    return 2 * x

# 主程序
if __name__ == "__main__":
    # 初始点
    x0 = [1.0, 1.0]
    # 调用最速下降法求解
    result = steepest_descent(example_f, example_grad_f, x0)
    print("最优解:", result)
    print("最优值:", example_f(result))

最优步长的求解

以下是使用黄金分割法和斐波那契法求步长${\lambda }^{(k)}$的具体方法：

黄金分割法

1. 原理
   • 黄金分割法是基于黄金分割比例来逐步缩小搜索区间，以逼近函数的最优解。对于一个单峰函数，在搜索区间$[a,b]$内，通过在两个特定点$x_1$和$x_2$处计算函数值，根据函数值的大小关系，不断舍弃不包含最优解的区间，从而使搜索区间不断缩小。
   • 这两个特定点$x_1$和$x_2$与区间端点的距离满足黄金分割比例，即$\frac{x_1-a}{b-a}=\frac{b-x_2}{b-a}=\frac{\sqrt{5}-1}{2}\approx 0.618$。
2. 计算步骤
   • 首先确定初始搜索区间$[a,b]$，即确定步长$\lambda$的取值范围，保证目标函数$f(x^{(k)}+\lambda P^{(k)})$在该区间内是单峰的。
   • 计算区间内的两个试探点$x_1=a+0.382(b-a)$，$x_2=a+0.618(b-a)$，并计算对应的函数值$f_1=f(x^{(k)}+x_1{P}^{(k)})$，$f_2=f(x^{(k)}+x_2{P}^{(k)})$。
   • 比较$f_1$和$f_2$的大小，如果$f_1<f_2$，则新的搜索区间为$[a,x_2]$；否则，新的搜索区间为$[x_1,b]$。
   • 重复上述步骤，不断缩小搜索区间，直到满足预设的精度要求，此时的区间中点或其中一个端点即可作为近似的最优步长${\lambda }^{(k)}$。

斐波那契法

1. 原理
   • 斐波那契法也是一种用于一维搜索的区间收缩方法，它与黄金分割法类似，但在确定试探点的位置时，是根据斐波那契数列来确定的。斐波那契数列满足$F_n=F_{n-1}+F_{n-2}$，其中$F_0=0$，$F_1=1$。
   • 在搜索过程中，利用斐波那契数列的性质来确定试探点的位置，使得在给定的搜索次数下，能够最大程度地缩小搜索区间。
2. 计算步骤
   • 同样先确定初始搜索区间$[a,b]$和搜索次数$n$。
   • 根据斐波那契数列计算出第$n$项$F_n$以及前一项$F_{n-1}$，第一个试探点为$x_1=a+\frac{F_{n-1}}{F_n}(b-a)$，第二个试探点为$x_2=a+\frac{F_{n-2}}{F_n}(b-a)$，计算对应的函数值$f_1=f(x^{(k)}+x_1{P}^{(k)})$，$f_2=f(x^{(k)}+x_2{P}^{(k)})$。
   • 比较$f_1$和$f_2$的大小，根据大小关系舍弃相应的区间，更新搜索区间和试探点。
   • 随着搜索次数的增加，不断重复上述步骤，当达到预定的搜索次数$n$时，得到的搜索区间的中点或其中一个端点就是近似的最优步长${\lambda }^{(k)}$。

牛顿法

牛顿法基本思想以目标函数的一个二次函数作为其近似，然后求出这个二次函数的极小值点，从而该极小值点近似为原目标函数的一个局部极小值点。
设$f(x)$是二次可微函数。假设上一步迭代到了$x^{(k)}$。我们在点$x^{(k)}$处$Taylor$展开，取二阶近似:
$$f\left(x\right)\approx Q\left(x\right)=f\left(x^{\left(k\right)}\right)+\nabla f{\left(x^{\left(k\right)}\right)}^T\left(x-x^{\left(k\right)}\right)+\frac{1}{2}{\left(x-x^{\left(k\right)}\right)}^T{\nabla }^{2}f\left(x^{\left(k\right)}\right)\left(x-x^{\left(k\right)}\right)$$
$$g\left(x^{\left(k\right)}\right)\stackrel{\Delta }{=}\nabla f\left(x^{\left(k\right)}\right),H\left(x^{\left(k\right)}\right)\stackrel{\Delta }{=}{\nabla }^{2}f\left(x^{\left(k\right)}\right)$$
反解可得到：
$$x=x^{\left(k\right)}-{\nabla }^{2}f{\left(x^{\left(k\right)}\right)}^{-1}\nabla f\left(x^{\left(k\right)}\right)$$
基于上面给出的信息，我们得到牛顿算法：
（1）给定初始点$x^{\left(1\right)}\in R^n$，给定误差$ϵ>0$,令$k=1$
（2）由${\nabla }^{2}f\left(x^{\left(k\right)}\right)d^{\left(k\right)}=-\nabla f\left(x^{\left(k\right)}\right)$，解线性方程组得到搜索方向$d^{(k)}$
（3）置$x^{\left(k+1\right)}=x^{\left(k\right)}+d^{\left(k\right)}$，即$x^{(k+1)}=x^{\left(k\right)}-{\nabla }^{2}f{\left(x^{\left(k\right)}\right)}^{-1}\nabla f\left(x^{\left(k\right)}\right)$
（4）当$\Vert d^{\left(k\right)}\Vert <ϵ$，迭代终止，否则，$k=k+1$,跳转到$2$

import numpy as np

def newton_method(f, grad_f, hess_f, x1, epsilon=1e-6, max_iter=100):
    """
    牛顿算法实现

    参数:
    f: 目标函数
    grad_f: 目标函数的梯度
    hess_f: 目标函数的海森矩阵
    x1: 初始点，numpy 数组
    epsilon: 误差阈值，默认为 1e-6
    max_iter: 最大迭代次数，默认为 100

    返回:
    x: 最优解
    iter_count: 迭代次数
    """
    x = x1.copy()
    k = 0
    while k < max_iter:
        # 步骤 2: 求解线性方程组得到搜索方向
        grad = grad_f(x)
        hess = hess_f(x)
        try:
            # 解线性方程组 H * d = -grad
            d = np.linalg.solve(hess, -grad)
        except np.linalg.LinAlgError:
            print("海森矩阵不可逆，无法求解线性方程组。")
            break

        # 步骤 3: 更新 x
        x_new = x + d

        # 步骤 4: 判断是否满足终止条件
        if np.linalg.norm(d) < epsilon:
            break

        x = x_new
        k += 1

    return x, k

# 示例：定义一个简单的目标函数及其梯度和海森矩阵
def f(x):
    # 目标函数 f(x) = x1^2 + x2^2
    return x[0]**2 + x[1]**2

def grad_f(x):
    # 目标函数的梯度
    return np.array([2 * x[0], 2 * x[1]])

def hess_f(x):
    # 目标函数的海森矩阵
    return np.array([[2, 0], [0, 2]])

# 初始点
x1 = np.array([1.0, 1.0])

# 调用牛顿算法
optimal_x, iter_count = newton_method(f, grad_f, hess_f, x1)

print("最优解:", optimal_x)
print("迭代次数:", iter_count)

阻尼牛顿法

阻尼牛顿法与原始牛顿法的主要区别在于增加了沿着牛顿方向进行一维搜索，其迭代公式:
$$x^{\left(k+1\right)}=x^{\left(k\right)}+{\lambda }^{\left(k\right)}{d}^{\left(k\right)}$$

• （1）给定初始点$x^{\left(1\right)}\in R^n$，给定误差$ϵ>0$,令$k=1$

• （2）由${\nabla }^{2}f\left(x^{\left(k\right)}\right)d^{\left(k\right)}=-\nabla f\left(x^{\left(k\right)}\right)$，解线性方程组得到搜索方向$d^{(k)}$

• （3） 从$x^{(k)}$出发，沿着$d^{(k)}$方向一维搜索：$$\underset{\lambda }{\min }f\left(x^{\left(k\right)}+\lambda d^{\left(k\right)}\right)$$

• （4）置$x^{\left(k+1\right)}=x^{\left(k\right)}+d^{\left(k\right)}$，即$x^{(k+1)}=x^{\left(k\right)}-{\lambda }^{\left(k\right)}{\nabla }^{2}f{\left(x^{\left(k\right)}\right)}^{-1}\nabla f\left(x^{\left(k\right)}\right)$

• （5）当$\Vert d^{\left(k\right)}\Vert <ϵ$，迭代终止，否则，$k=k+1$,跳转到$2$

模式搜索法

Powell方法

Powell方法是研究对称正定的二次函数：
$$f(x)=\frac{1}{2}{x}^{T}Qx+b^{T}x+c$$
极小化问题时候形成的。

• （1）选定初始点$x^{\left(0\right)}$，$n$个线性无关向量$d^{(i)}$,，给定误差$ϵ>0$,令$k=1$,$i=1,2,...,n$
• （2）依次对目标函数作直线搜索$$\{\begin{array}{c}f\left(x^{\left(i-1\right)}+{\lambda }_i{d}^{\left(i\right)}\right)=\underset{\lambda \in R}{\min }f\left(x^{\left(i-1\right)}+\lambda d^{\left(i\right)}\right)\\ x^{\left(i\right)}=x^{\left(i-1\right)}+{\lambda }_i{d}^{\left(i\right)}\end{array}$$
• （3） 对$i=1,2,...,n-1$，$d^{(i)}=d^{(i+1)}$,$d^{(n+1)}=x^{(n)}-x^{(0)}$,置$d^{(n)}=d^{(n+1)}$
• （4）直线搜索$x^{\left(n+1\right)}=x^{\left(n\right)}+{\lambda }_{n+1}{d}^{\left(n+1\right)}$，$$f\left(x^{\left(n\right)}+{\lambda }_{n+1}{d}^{\left(n+1\right)}\right)=\min f\left(x^{\left(n\right)}+\lambda d^{\left(n+1\right)}\right)$$
• （5）$\Vert x^{\left(n+1\right)}-x^{\left(0\right)}\Vert <ϵ,x^{\left(n+1\right)}$就是所求的极小值点，则停止计算，否则转第$6$步
• （6）置$x^{(0)}=x^{(n+1)}$转$(2)$

有约束非线性规划

有约束非线性规划的一般数学模型为如下形式：
$$\begin{gathered} \min f\left(x\right) \\ s.t.\{\begin{array}{c}{g}_i\left(x\right)⩾0,i=1,2,...,m\\ h_j\left(x\right)=0,j=1,2,...,n\end{array} \end{gathered}$$
其中$f(x)，g_i(x),h_i(x)$至少有一个非线性函数。

外点罚函数法

外点罚函数法是一种将有约束优化问题转化为无约束优化问题的方法。其核心思想是通过在目标函数中添加一个罚项，当约束条件不满足时，罚项的值会变得很大，从而迫使无约束优化问题的解逐渐逼近原约束优化问题的可行域。

1. 对于只含等式约束的优化问题：
   $$\begin{gathered} \min f\left(x\right) \\ s.t.h_j\left(x\right)=0,j=1,2,...,l \end{gathered}$$
   采用如下方法建立罚函数：
   $$F\left(x,m_k\right)=f\left(x\right)+m_k\sum _{j=1}^l{h}_j^2\left(x\right)$$
   这样，我们就把原约束问题化为如下的无约束问题：
   $$minF(x,m_k)$$
2. 对于不等式条件：
   构造罚函数：
   $$F\left(x,m_k\right)=f\left(x\right)+m_k\sum _{i=1}^m{g}_i^2\left(x\right)u_i\left(g_i\right),u\text{为单位跃阶函数}$$
3. 综合以上思想，我们可以处理一般的罚函数
   对于包含等式约束 $h_j(x)=0$（$j=1,2,\cdots ,l$）和不等式约束 $g_i(x)\le 0$（$i=1,2,\cdots ,m$）的一般约束优化问题：
   $$\begin{array}{rl}\min & f(x)\\ \text{s.t.}& h_j(x)=0,j=1,2,\cdots ,l\\ & g_i(x)\le 0,i=1,2,\cdots ,m\end{array}$$
   可以构造外点罚函数：
   F ( x , m k ) = f ( x ) + m k ( ∑ j = 1 l h j 2 ( x ) + ∑ i = 1 m [ max ⁡ { 0 , g i ( x ) } ] 2 ) F(x, m_k) = f(x) + m_k \left( \sum_{j = 1}^l h_j^2(x) + \sum_{i = 1}^m [\max\{0, g_i(x)\}]^2 \right) F(x,mk​)=f(x)+mk​(j=1∑l​hj2​(x)+i=1∑m​[max{0,gi​(x)}]2)
   其中，$m_k$ 是罚因子，且满足 $0<m_1<m_2<\cdots <m_k<\cdots$，${\lim }_{k\to \infty }{m}_k=\infty$。

算法步骤

1. 初始化：选择初始点 $x^{(0)}$，初始罚因子 $m_1>0$，罚因子增长系数 $c>1$，收敛精度 $ϵ>0$，令 $k=1$。
2. 求解无约束优化问题：以 $x^{(k-1)}$ 为初始点，求解无约束优化问题 $\min F(x,m_k)$，得到最优解 $x^{(k)}$。
3. 判断收敛条件：如果满足 $\left|F(x^{(k)},m_k)-F(x^{(k-1)},m_{k-1})\right|<ϵ$ 或者约束条件的违反程度足够小（例如， ∑ j = 1 l h j 2 ( x ( k ) ) + ∑ i = 1 m [ max ⁡ { 0 , g i ( x ( k ) ) } ] 2 < ϵ \sum_{j = 1}^l h_j^2(x^{(k)}) + \sum_{i = 1}^m [\max\{0, g_i(x^{(k)})\}]^2 < \epsilon ∑j=1l​hj2​(x(k))+∑i=1m​[max{0,gi​(x(k))}]2<ϵ），则停止迭代，$x^{(k)}$ 即为原约束优化问题的近似最优解；否则，进入下一步。
4. 更新罚因子：令 $m_{k+1}=c{m}_k$，$k=k+1$，返回步骤 2。

Python 代码实现

import numpy as np
from scipy.optimize import minimize

def penalty_function(x, f, h, g, mk):
    """
    计算外点罚函数的值
    :param x: 当前点
    :param f: 目标函数
    :param h: 等式约束函数列表
    :param g: 不等式约束函数列表
    :param mk: 罚因子
    :return: 罚函数的值
    """
    penalty_term = 0
    # 等式约束的罚项
    for hj in h:
        penalty_term += hj(x) ** 2
    # 不等式约束的罚项
    for gi in g:
        penalty_term += max(0, gi(x)) ** 2
    return f(x) + mk * penalty_term

def outer_penalty_method(f, h, g, x0, m1=1.0, c=10.0, epsilon=1e-6, max_iter=100):
    """
    外点罚函数法求解约束优化问题
    :param f: 目标函数
    :param h: 等式约束函数列表
    :param g: 不等式约束函数列表
    :param x0: 初始点
    :param m1: 初始罚因子
    :param c: 罚因子增长系数
    :param epsilon: 收敛精度
    :param max_iter: 最大迭代次数
    :return: 最优解
    """
    mk = m1
    x = np.array(x0)
    for k in range(max_iter):
        # 定义当前罚函数
        def current_penalty(x):
            return penalty_function(x, f, h, g, mk)
        # 求解无约束优化问题
        result = minimize(current_penalty, x)
        x_new = result.x
        # 判断收敛条件
        if k > 0 and np.abs(current_penalty(x_new) - current_penalty(x)) < epsilon:
            break
        x = x_new
        # 更新罚因子
        mk *= c
    return x

# 示例：求解一个简单的约束优化问题
# 目标函数
def f(x):
    return (x[0] - 1) ** 2 + (x[1] - 2.5) ** 2

# 等式约束
def h1(x):
    return x[0] + x[1] - 3

# 不等式约束
def g1(x):
    return -x[0] + 2

def g2(x):
    return -x[1] + 2

# 初始点
x0 = [0, 0]

# 调用外点罚函数法
optimal_x = outer_penalty_method(f, [h1], [g1, g2], x0)

print("最优解:", optimal_x)

注意事项

• 罚因子的增长系数 $c$ 和初始罚因子 $m_1$ 的选择会影响算法的收敛速度和稳定性，需要根据具体问题进行调整。
• 无约束优化问题的求解可以使用不同的优化算法，这里使用了 scipy.optimize.minimize 函数，你可以根据需要选择其他合适的算法。

内点罚函数法

内点罚函数法（Interior Point Penalty Method），也称为障碍函数法，是一种用于求解有约束优化问题的重要方法。它主要适用于只包含不等式约束的优化问题，基本思想是在可行域内部设置障碍，使得迭代点始终保持在可行域内，并且随着迭代的进行，障碍的作用逐渐减弱，从而使迭代点逼近最优解。以下将详细介绍内点罚函数法的原理、算法步骤，并给出 Python 代码示例。

原理

考虑如下只含不等式约束的优化问题：
$$\begin{array}{rl}\min & f(x)\\ \text{s.t.}& g_i(x)\le 0,i=1,2,\cdots ,m\end{array}$$
内点罚函数法通过构造障碍函数来实现。常见的障碍函数有倒数障碍函数和对数障碍函数。

倒数障碍函数

$$B(x,r_k)=f(x)+r_k\sum _{i=1}^m\frac{1}{-g_i(x)}$$

对数障碍函数

$$B(x,r_k)=f(x)-r_k\sum _{i=1}^m\ln (-g_i(x))$$

其中，$r_k$ 是障碍因子，且满足 $r_1>r_2>\cdots >r_k>\cdots$，${\lim }_{k\to \infty }{r}_k=0$。当迭代点接近可行域边界时，障碍函数的值会变得很大，从而阻止迭代点越过边界。

算法步骤

1. 初始化：选择一个初始可行点 $x^{(0)}$（即满足所有不等式约束 $g_i(x^{(0)})\le 0$，$i=1,2,\cdots ,m$），初始障碍因子 $r_1>0$，障碍因子缩减系数 $c\in (0,1)$，收敛精度 $ϵ>0$，令 $k=1$。
2. 求解无约束优化问题：以 $x^{(k-1)}$ 为初始点，求解无约束优化问题 $\min B(x,r_k)$，得到最优解 $x^{(k)}$。
3. 判断收敛条件：如果满足 $\left|r_k{\sum }_{i=1}^m\frac{1}{-g_i(x^{(k)})}\right|<ϵ$（对于倒数障碍函数）或者 $\left|-r_k{\sum }_{i=1}^m\ln (-g_i(x^{(k)})\right|<ϵ$（对于对数障碍函数），则停止迭代，$x^{(k)}$ 即为原约束优化问题的近似最优解；否则，进入下一步。
4. 更新障碍因子：令 $r_{k+1}=c{r}_k$，$k=k+1$，返回步骤 2。

Python 代码示例（使用对数障碍函数）

import numpy as np
from scipy.optimize import minimize

def logarithmic_barrier_function(x, f, g, rk):
    """
    计算对数障碍函数的值
    :param x: 当前点
    :param f: 目标函数
    :param g: 不等式约束函数列表
    :param rk: 障碍因子
    :return: 对数障碍函数的值
    """
    barrier_term = 0
    for gi in g:
        if gi(x) >= 0:
            # 如果不满足不等式约束，返回一个很大的值
            return np.inf
        barrier_term += np.log(-gi(x))
    return f(x) - rk * barrier_term

def interior_point_method(f, g, x0, r1=1.0, c=0.1, epsilon=1e-6, max_iter=100):
    """
    内点罚函数法求解约束优化问题
    :param f: 目标函数
    :param g: 不等式约束函数列表
    :param x0: 初始可行点
    :param r1: 初始障碍因子
    :param c: 障碍因子缩减系数
    :param epsilon: 收敛精度
    :param max_iter: 最大迭代次数
    :return: 最优解
    """
    rk = r1
    x = np.array(x0)
    for k in range(max_iter):
        # 定义当前对数障碍函数
        def current_barrier(x):
            return logarithmic_barrier_function(x, f, g, rk)
        # 求解无约束优化问题
        result = minimize(current_barrier, x)
        x_new = result.x
        # 判断收敛条件
        barrier_term = 0
        for gi in g:
            barrier_term += np.log(-gi(x_new))
        if np.abs(-rk * barrier_term) < epsilon:
            break
        x = x_new
        # 更新障碍因子
        rk *= c
    return x

# 示例：求解一个简单的约束优化问题
# 目标函数
def f(x):
    return (x[0] - 1) ** 2 + (x[1] - 2.5) ** 2

# 不等式约束
def g1(x):
    return -x[0] + 2

def g2(x):
    return -x[1] + 2

# 初始可行点
x0 = [1, 1]

# 调用内点罚函数法
optimal_x = interior_point_method(f, [g1, g2], x0)

print("最优解:", optimal_x)

代码说明

1. logarithmic_barrier_function 函数：用于计算对数障碍函数的值。如果迭代点不满足不等式约束，返回一个很大的值，以避免该点被选为最优解。
2. interior_point_method 函数：实现了内点罚函数法的主要步骤，包括初始化、求解无约束优化问题、判断收敛条件和更新障碍因子。
3. 示例：定义了一个简单的约束优化问题，包含两个不等式约束，调用 interior_point_method 函数求解最优解。

注意事项

• 初始可行点的选择非常重要，如果选择不当，可能会导致算法无法收敛。
• 障碍因子的缩减系数 $c$ 和初始障碍因子 $r_1$ 的选择会影响算法的收敛速度和稳定性，需要根据具体问题进行调整。