OPTIMIZATION Algorithms and Applications

二次型转换

1. 将函数转换为矩阵形式

给定函数：

f(x1,x2)=6x12−6x1x2+2x22−x1−2x2f(x_1, x_2) = 6x_1^2 - 6x_1x_2 + 2x_2^2 - x_1 - 2x_2

将其表示为矩阵形式：

f(x1,x2)=12XTAX+BTXf(x_1, x_2) = \frac{1}{2} X^T A X + B^T X

其中：

X=[x1x2],A=[12−6−64],B=[−1−2]X = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \quad A = \begin{bmatrix} 12 & -6 \\ -6 & 4 \end{bmatrix}, \quad B = \begin{bmatrix} -1 \\ -2 \end{bmatrix}

2. 对角化矩阵 AA

2.1 计算特征值

解特征方程 det⁡(A−λI)=0\det(A - \lambda I) = 0：

det⁡[12−λ−6−64−λ]=(12−λ)(4−λ)−(−6)2=λ2−16λ+12=0\det \begin{bmatrix} 12 - \lambda & -6 \\ -6 & 4 - \lambda \end{bmatrix} = (12 - \lambda)(4 - \lambda) - (-6)^2 = \lambda^2 - 16\lambda + 12 = 0

解得特征值：

λ1=8+213,λ2=8−213\lambda_1 = 8 + 2\sqrt{13}, \quad \lambda_2 = 8 - 2\sqrt{13}

2.2 计算特征向量

对于 λ1=8+213\lambda_1 = 8 + 2\sqrt{13}：
解方程 (A−λ1I)v=0(A - \lambda_1 I)v = 0：

[12−(8+213)−6−64−(8+213)][v1v2]=0\begin{bmatrix} 12 - (8 + 2\sqrt{13}) & -6 \\ -6 & 4 - (8 + 2\sqrt{13}) \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = 0

解得特征向量：

v1=[14+2136]v_1 = \begin{bmatrix} 1 \\ \frac{4 + 2\sqrt{13}}{6} \end{bmatrix}

对于 λ2=8−213\lambda_2 = 8 - 2\sqrt{13}：
解方程 (A−λ2I)v=0(A - \lambda_2 I)v = 0：

[12−(8−213)−6−64−(8−213)][v1v2]=0\begin{bmatrix} 12 - (8 - 2\sqrt{13}) & -6 \\ -6 & 4 - (8 - 2\sqrt{13}) \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = 0

解得特征向量：

v2=[14−2136]v_2 = \begin{bmatrix} 1 \\ \frac{4 - 2\sqrt{13}}{6} \end{bmatrix}

2.3 归一化特征向量

将特征向量归一化为单位向量：

对于 v1v_1：

∥v1∥=12+(4+2136)2\|v_1\| = \sqrt{1^2 + \left(\frac{4 + 2\sqrt{13}}{6}\right)^2}

归一化后的特征向量：

u1=v1∥v1∥u_1 = \frac{v_1}{\|v_1\|}

对于 v2v_2：

∥v2∥=12+(4−2136)2\|v_2\| = \sqrt{1^2 + \left(\frac{4 - 2\sqrt{13}}{6}\right)^2}

归一化后的特征向量：

u2=v2∥v2∥u_2 = \frac{v_2}{\|v_2\|}

2.4 对角化矩阵

将归一化后的特征向量组成矩阵 PP：

P=[u1u2]P = \begin{bmatrix} u_1 & u_2 \end{bmatrix}

对角矩阵 DD 为：

D=[8+213008−213]D = \begin{bmatrix} 8 + 2\sqrt{13} & 0 \\ 0 & 8 - 2\sqrt{13} \end{bmatrix}

因此，矩阵 AA 可以表示为：

A=PDP−1A = P D P^{-1}

3. 绘制对角化后的矩阵 DD 的等高线图像

3.1 等高线的类型

矩阵 DD 的等高线方程为：

12XTDX=常数\frac{1}{2} X^T D X = \text{常数}

展开后为：

12((8+213)x12+(8−213)x22)=常数\frac{1}{2} \left( (8 + 2\sqrt{13}) x_1^2 + (8 - 2\sqrt{13}) x_2^2 \right) = \text{常数}

这是一个标准的椭圆方程，其主轴与坐标轴对齐。

3.2 确定主轴的方向和大小

  • 主轴方向：由于 DD 是对角矩阵，主轴方向与坐标轴 x1x_1 和 x2x_2 对齐。
  • 主轴长度：
    • 长轴长度：18+213\frac{1}{\sqrt{8 + 2\sqrt{13}}}
    • 短轴长度：18−213\frac{1}{\sqrt{8 - 2\sqrt{13}}}

3.3 绘制草图

  1. 绘制坐标系：在纸上绘制 x1x_1 和 x2x_2 轴。
  2. 绘制主轴：
    • 长轴沿 x1x_1 轴，长度为 18+213\frac{1}{\sqrt{8 + 2\sqrt{13}}}。
    • 短轴沿 x2x_2 轴，长度为 18−213\frac{1}{\sqrt{8 - 2\sqrt{13}}}。
  3. 绘制椭圆：
    • 以原点为中心，根据主轴长度绘制椭圆。
  4. 调整形状：
    • 确保椭圆的长轴和短轴长度符合计算结果。

4. 总结

  • 矩阵 AA 对角化后的矩阵 DD 为：

    D=[8+213008−213]D = \begin{bmatrix} 8 + 2\sqrt{13} & 0 \\ 0 & 8 - 2\sqrt{13} \end{bmatrix}

  • 特征向量归一化后组成矩阵 PP。

  • 等高线是标准椭圆，主轴与坐标轴对齐。

  • 长轴长度为 18+213\frac{1}{\sqrt{8 + 2\sqrt{13}}}，短轴长度为 18−213\frac{1}{\sqrt{8 - 2\sqrt{13}}}。

等高线：

  • 如果 λ1 和 λ2 都为正，则等高线是椭圆。

  • 如果 λ1 和 λ2 一正一负，则等高线是双曲线。

  • 如果 λ1 或 λ2 为零，则等高线是抛物线。

  • 长轴方向对应较小的特征值 方向为该特征值的特征向量

  • 短轴方向对应较大的特征值 方向为该特征值的特征向量

第一章 介绍

重点：
凸函数、泰勒展开

课后习题 20 题

题目：计算函数 $$ f(x) = x_1^2 x_2 + x_2^2 x_3 - x_1 x_2 x_3^2 $$ 在点 (1,1,−1)(1, 1, -1) 处的方向导数，方向向量为 d=[123]d = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}。

1. 方向导数的定义

方向导数表示函数 f(x)f(x) 在点 xx 处沿方向 dd 的瞬时变化率。方向导数的计算公式为：

Ddf(x)=∇f(x)TuD_d f(x) = \nabla f(x)^T u

其中：

  • ∇f(x)\nabla f(x) 是函数 f(x)f(x) 在点 xx 处的梯度向量。
  • uu 是方向向量 dd 的单位向量，即 u=d∥d∥u = \frac{d}{\|d\|}。

2. 计算梯度向量 $$\nabla f(x)$$

首先，我们需要计算函数 $$ f(x) = x_1^2 x_2 + x_2^2 x_3 - x_1 x_2 x_3^2 $$ 的梯度向量 $$ \nabla f(x) $$。梯度向量的每个分量是函数对相应变量的偏导数。

  1. 对 x1x_1 的偏导数：

    ∂f∂x1=2x1x2−x2x32\frac{\partial f}{\partial x_1} = 2x_1 x_2 - x_2 x_3^2

  2. 对 x2x_2 的偏导数：

    ∂f∂x2=x12+2x2x3−x1x32\frac{\partial f}{\partial x_2} = x_1^2 + 2x_2 x_3 - x_1 x_3^2

  3. 对 x3x_3 的偏导数：

    ∂f∂x3=x22−2x1x2x3\frac{\partial f}{\partial x_3} = x_2^2 - 2x_1 x_2 x_3

因此，梯度向量为：

∇f(x)=[2x1x2−x2x32x12+2x2x3−x1x32x22−2x1x2x3]\nabla f(x) = \begin{bmatrix} 2x_1 x_2 - x_2 x_3^2 \\ x_1^2 + 2x_2 x_3 - x_1 x_3^2 \\ x_2^2 - 2x_1 x_2 x_3 \end{bmatrix}

3. 计算梯度在点 (1,1,−1)(1, 1, -1) 处的值

将 x=(1,1,−1)x = (1, 1, -1) 代入梯度向量：

  1. 对 x1x_1 的偏导数：

    ∂f∂x1=2(1)(1)−(1)(−1)2=2−1=1\frac{\partial f}{\partial x_1} = 2(1)(1) - (1)(-1)^2 = 2 - 1 = 1

  2. 对 x2x_2 的偏导数：

    ∂f∂x2=(1)2+2(1)(−1)−(1)(−1)2=1−2−1=−2\frac{\partial f}{\partial x_2} = (1)^2 + 2(1)(-1) - (1)(-1)^2 = 1 - 2 - 1 = -2

  3. 对 x3x_3 的偏导数：

    ∂f∂x3=(1)2−2(1)(1)(−1)=1+2=3\frac{\partial f}{\partial x_3} = (1)^2 - 2(1)(1)(-1) = 1 + 2 = 3

因此，梯度向量在点 (1,1,−1)(1, 1, -1) 处的值为：

∇f(1,1,−1)=[1−23]\nabla f(1, 1, -1) = \begin{bmatrix} 1 \\ -2 \\ 3 \end{bmatrix}

4. 计算方向向量 dd 的单位向量 uu

方向向量为 d=[123]d = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}，其单位向量 uu 为：

u=d∥d∥=112+22+32[123]=114[123]u = \frac{d}{\|d\|} = \frac{1}{\sqrt{1^2 + 2^2 + 3^2}} \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix} = \frac{1}{\sqrt{14}} \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}

5. 计算方向导数

方向导数为梯度向量与单位方向向量的点积：

Ddf(1,1,−1)=∇f(1,1,−1)Tu=[1−23]⋅114[123]D_d f(1, 1, -1) = \nabla f(1, 1, -1)^T u = \begin{bmatrix} 1 & -2 & 3 \end{bmatrix} \cdot \frac{1}{\sqrt{14}} \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}

计算点积：

Ddf(1,1,−1)=114(1⋅1+(−2)⋅2+3⋅3)=114(1−4+9)=614D_d f(1, 1, -1) = \frac{1}{\sqrt{14}} (1 \cdot 1 + (-2) \cdot 2 + 3 \cdot 3) = \frac{1}{\sqrt{14}} (1 - 4 + 9) = \frac{6}{\sqrt{14}}

6. 最终答案

函数 $$ f(x) = x_1^2 x_2 + x_2^2 x_3 - x_1 x_2 x_3^2 $$ 在点 (1,1,−1)(1, 1, -1) 处沿方向 d=[123]d = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix} 的方向导数为：

Ddf(1,1,−1)=614D_d f(1, 1, -1) = \frac{6}{\sqrt{14}}

总结

  • 方向导数表示函数在某一方向上的瞬时变化率。
  • 梯度向量是函数在某一点的导数向量，表示函数在该点的变化率。
  • 单位方向向量是方向向量的归一化形式。
  • 方向导数的计算方法是梯度向量与单位方向向量的点积。

第二章 一维优化算法

1. 引言

一维（1-D）优化问题指的是只有一个变量的目标函数。尽管在实际问题中，多变量优化问题更为常见，但一维优化算法是多变量优化算法的基础。一维优化算法可以分为基于梯度的算法和非基于梯度的算法。本章将讨论一些流行的算法。

例子：
一个单变量目标函数可以是：

f(x)=2x2−2x+8f(x) = 2x^2 - 2x + 8

这是一个无约束优化问题，目标是找到使 f(x)f(x) 最小的 xx。如果 xx 被限制在 a≤x≤ba \leq x \leq b 之间，则成为有约束优化问题。

单调函数：如果函数 f(x)f(x) 在两点 aa 和 bb 之间连续增加或减少，则称为单调函数（见图2.1）。

单峰函数：在单峰函数中，函数在其最小点 x∗x^* 的两侧是单调的。图2.2展示了函数 f(x)=2x2−2x+8f(x) = 2x^2 - 2x + 8 的图形，可以看出这是一个单峰函数。

2. 测试问题

在讨论优化算法之前，我们先定义一个测试问题。太阳能问题定义为成本最小化问题，成本是存储系统体积和集热器表面积的函数。体积和表面积是设计变量温度 TT 的函数。成本函数可以表示为：

U=204,165.5330−2T+10,400T−20U = \frac{204,165.5}{330 - 2T} + \frac{10,400}{T - 20}

变量 TT 被限制在 40∘C40^\circ C 到 90∘C90^\circ C 之间。图2.4展示了成本函数随 TT 变化的图形，最小值出现在 T∗=55.08T^* = 55.08，最小值为 U∗=1225.166U^* = 1225.166。

3. 求解技术

一维优化问题的求解技术可以分为基于梯度的算法和非基于梯度的算法。

3.1 二分法（Bisection Method）

核心思想：
二分法是一种基于导数的优化方法。它利用函数导数的符号变化来确定最小值所在的区间。对于一个单峰函数，导数为零的点就是最小值点。二分法通过不断缩小搜索区间来逼近最小值。

算法步骤：

  1. 给定初始区间 [a,b][a, b] 和精度要求 ϵ\epsilon。
  2. 计算区间的中点 α=a+b2\alpha = \frac{a + b}{2}，并计算 f′(α)f'(\alpha)。
  3. 如果 f′(a)⋅f′(α)<0f'(a) \cdot f'(\alpha) < 0，说明最小值在 [a,α][a, \alpha] 区间内，更新 b=αb = \alpha。
  4. 否则，最小值在 [α,b][\alpha, b] 区间内，更新 a=αa = \alpha。
  5. 重复上述步骤，直到区间长度 ∣a−b∣<ϵ|a - b| < \epsilon，此时认为算法收敛，输出最小值点 x∗=ax^* = a。

优点：

  • 简单易实现。
  • 对于单峰函数，能够稳定收敛。

缺点：

  • 需要计算导数，适用于导数容易求解的函数。
  • 收敛速度较慢。

3.2 牛顿-拉夫森法（Newton-Raphson Method）

核心思想：
牛顿-拉夫森法是一种基于二阶导数的优化方法。它通过泰勒展开式近似函数，并利用函数的导数信息来快速逼近最小值点。

算法步骤：

  1. 给定初始点 xx 和精度要求 ϵ\epsilon。
  2. 计算 f′(x)f'(x) 和 f′′(x)f''(x)。
  3. 更新 xx 的值：

    xk+1=xk−f′(xk)f′′(xk)x_{k+1} = x_k - \frac{f'(x_k)}{f''(x_k)}

  4. 重复上述步骤，直到 ∣xk+1−xk∣<ϵ|x_{k+1} - x_k| < \epsilon，输出最小值点 x∗=xk+1x^* = x_{k+1}。

优点：

  • 收敛速度快（二次收敛）。
  • 对于光滑函数，能够快速找到最小值。

缺点：

  • 需要计算二阶导数，且二阶导数必须存在。
  • 对初始猜测敏感，初始点选择不当可能导致算法发散。

3.3 割线法（Secant Method）

核心思想：
割线法是一种基于导数的优化方法，类似于牛顿-拉夫森法，但它不需要计算二阶导数。它通过两个点的导数信息来逼近最小值点。

算法步骤：

  1. 给定初始区间 [a,b][a, b] 和精度要求 ϵ\epsilon。
  2. 计算 f′(a)f'(a) 和 f′(b)f'(b)。
  3. 如果 f′(a)⋅f′(b)<0f'(a) \cdot f'(b) < 0，说明最小值在 [a,b][a, b] 区间内。
  4. 通过割线公式计算新的点 α\alpha：

    α=x2−f′(x2)(f′(x2)−f′(x1))/(x2−x1)\alpha = x_2 - \frac{f'(x_2)}{(f'(x_2) - f'(x_1))/(x_2 - x_1)}

  5. 根据 f′(α)f'(\alpha) 的符号更新区间 [a,b][a, b]。
  6. 重复上述步骤，直到 ∣f′(α)∣<ϵ|f'(\alpha)| < \epsilon，输出最小值点 x∗=αx^* = \alpha。

优点：

  • 不需要计算二阶导数。
  • 对于导数容易计算的函数，收敛速度较快。

缺点：

  • 收敛速度比牛顿-拉夫森法慢。
  • 对初始区间选择敏感。

3.4 三次多项式拟合（Cubic Polynomial Fit）

核心思想：
三次多项式拟合是一种基于函数值和导数信息的优化方法。它通过拟合一个三次多项式来逼近函数，并利用多项式的性质找到最小值点。

算法步骤：

  1. 给定初始区间 [a,b][a, b] 和精度要求 ϵ\epsilon。
  2. 计算 f(a)f(a)、f′(a)f'(a)、f(b)f(b) 和 f′(b)f'(b)。
  3. 通过拟合三次多项式 P(x)=a0+a1x+a2x2+a3x3P(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3，找到多项式的极小值点 xˉ\bar{x}。
  4. 根据 f′(xˉ)f'(\bar{x}) 的符号更新区间 [a,b][a, b]。
  5. 重复上述步骤，直到 ∣f′(xˉ)∣<ϵ|f'(\bar{x})| < \epsilon，输出最小值点 x∗=xˉx^* = \bar{x}。

优点：

  • 收敛速度快。
  • 适用于光滑函数。

缺点：

  • 需要计算函数值和导数。
  • 对于非光滑函数，拟合效果可能较差。
举例：
算法步骤

  1. 初始化区间：

    • 给定初始区间 [a,b][a, b] 和精度要求 ϵ\epsilon。
    • 计算 f(a)f(a)、f′(a)f'(a)、f(b)f(b) 和 f′(b)f'(b)。

  2. 拟合三次多项式：

    • 假设函数 f(x)f(x) 在区间 [a,b][a, b] 内可以用一个三次多项式 P(x)P(x) 近似：

      P(x)=a0+a1x+a2x2+a3x3P(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3

    • 通过已知的函数值和导数信息，建立方程组求解系数 a0,a1,a2,a3a_0, a_1, a_2, a_3。

  3. 求解多项式系数：

    • 利用 f(a)f(a)、f′(a)f'(a)、f(b)f(b) 和 f′(b)f'(b)，建立以下四个方程：

      {P(a)=f(a)P′(a)=f′(a)P(b)=f(b)P′(b)=f′(b)\begin{cases} P(a) = f(a) \\ P'(a) = f'(a) \\ P(b) = f(b) \\ P'(b) = f'(b) \end{cases}

    • 代入 P(x)P(x) 和 P′(x)P'(x) 的表达式，得到：

      {a0+a1a+a2a2+a3a3=f(a)a1+2a2a+3a3a2=f′(a)a0+a1b+a2b2+a3b3=f(b)a1+2a2b+3a3b2=f′(b)\begin{cases} a_0 + a_1 a + a_2 a^2 + a_3 a^3 = f(a) \\ a_1 + 2a_2 a + 3a_3 a^2 = f'(a) \\ a_0 + a_1 b + a_2 b^2 + a_3 b^3 = f(b) \\ a_1 + 2a_2 b + 3a_3 b^2 = f'(b) \end{cases}

    • 这是一个线性方程组，可以通过矩阵求解法或代入法求解 a0,a1,a2,a3a_0, a_1, a_2, a_3。

  4. 找到多项式的极小值点：

    • 对 P(x)P(x) 求导，得到导函数：

      P′(x)=a1+2a2x+3a3x2P'(x) = a_1 + 2a_2 x + 3a_3 x^2

    • 解方程 P′(x)=0P'(x) = 0，得到多项式的极值点 xˉ\bar{x}。

  5. 更新区间：

    • 根据 f′(xˉ)f'(\bar{x}) 的符号更新区间 [a,b][a, b]：
      • 如果 f′(xˉ)<0f'(\bar{x}) < 0，则最小值在 [xˉ,b][\bar{x}, b] 区间内，更新 a=xˉa = \bar{x}。
      • 如果 f′(xˉ)>0f'(\bar{x}) > 0，则最小值在 [a,xˉ][a, \bar{x}] 区间内，更新 b=xˉb = \bar{x}。

  6. 重复迭代：

    • 重复上述步骤，直到 ∣f′(xˉ)∣<ϵ|f'(\bar{x})| < \epsilon，输出最小值点 x∗x^* 和最小值 f(x∗)f(x^*)。
举例说明

假设我们有以下函数和区间：

f(x)=3ex−x3+5x,[−3,3]f(x) = 3e^x - x^3 + 5x, \quad [-3, 3]

  1. 初始化区间：

    • a=−3a = -3，b=3b = 3。
    • 计算 f(a)f(a)、f′(a)f'(a)、f(b)f(b) 和 f′(b)f'(b)：

      f(−3)=3e−3−(−3)3+5(−3)≈0.149+27−15=12.149f′(−3)=3e−3−3(−3)2+5≈0.149−27+5=−21.851f(3)=3e3−33+5(3)≈60.257−27+15=48.257f′(3)=3e3−3(3)2+5≈60.257−27+5=38.257f(-3) = 3e^{-3} - (-3)^3 + 5(-3) \approx 0.149 + 27 - 15 = 12.149 \\ f'(-3) = 3e^{-3} - 3(-3)^2 + 5 \approx 0.149 - 27 + 5 = -21.851 \\ f(3) = 3e^{3} - 3^3 + 5(3) \approx 60.257 - 27 + 15 = 48.257 \\ f'(3) = 3e^{3} - 3(3)^2 + 5 \approx 60.257 - 27 + 5 = 38.257

  2. 拟合三次多项式：

    • 建立方程组：

      {a0+a1(−3)+a2(−3)2+a3(−3)3=12.149a1+2a2(−3)+3a3(−3)2=−21.851a0+a1(3)+a2(3)2+a3(3)3=48.257a1+2a2(3)+3a3(3)2=38.257\begin{cases} a_0 + a_1 (-3) + a_2 (-3)^2 + a_3 (-3)^3 = 12.149 \\ a_1 + 2a_2 (-3) + 3a_3 (-3)^2 = -21.851 \\ a_0 + a_1 (3) + a_2 (3)^2 + a_3 (3)^3 = 48.257 \\ a_1 + 2a_2 (3) + 3a_3 (3)^2 = 38.257 \end{cases}

    • 化简方程组：

      {a0−3a1+9a2−27a3=12.149a1−6a2+27a3=−21.851a0+3a1+9a2+27a3=48.257a1+6a2+27a3=38.257\begin{cases} a_0 - 3a_1 + 9a_2 - 27a_3 = 12.149 \\ a_1 - 6a_2 + 27a_3 = -21.851 \\ a_0 + 3a_1 + 9a_2 + 27a_3 = 48.257 \\ a_1 + 6a_2 + 27a_3 = 38.257 \end{cases}

  3. 求解系数：

    • 通过解线性方程组，得到：

      a0≈10.0,a1≈5.0,a2≈−1.0,a3≈0.5a_0 \approx 10.0, \quad a_1 \approx 5.0, \quad a_2 \approx -1.0, \quad a_3 \approx 0.5

  4. 找到极小值点：

    • 多项式的导函数为：

      P′(x)=5−2x+1.5x2P'(x) = 5 - 2x + 1.5x^2

    • 解方程 P′(x)=0P'(x) = 0，得到极值点 xˉ≈−0.5\bar{x} \approx -0.5。

  5. 更新区间：

    • 计算 f′(xˉ)f'(\bar{x})，根据符号更新区间。

  6. 输出结果：

    • 最小值点 x∗≈−0.5x^* \approx -0.5，最小值 f(x∗)≈1.5f(x^*) \approx 1.5。

3.5 黄金分割法（Golden Section Method）

核心思想：
黄金分割法是一种非基于梯度的优化方法。它通过函数值的比较来逐步缩小搜索区间，最终找到最小值点。黄金分割法的核心思想是利用黄金比例（1.618）来分割区间。

算法步骤：

  1. 给定初始区间 [a,b][a, b] 和精度要求 ϵ\epsilon。
  2. 计算两个内点：

    a1=a+(1−τ)(b−a),a2=a+τ(b−a)a_1 = a + (1 - \tau)(b - a), \quad a_2 = a + \tau(b - a)

    其中 τ=5−12≈0.618\tau = \frac{\sqrt{5} - 1}{2} \approx 0.618。
  3. 比较 f(a1)f(a_1) 和 f(a2)f(a_2)：
    • 如果 f(a1)>f(a2)f(a_1) > f(a_2)，则最小值在 [a1,b][a_1, b] 区间内，更新 a=a1a = a_1。
    • 否则，最小值在 [a,a2][a, a_2] 区间内，更新 b=a2b = a_2。
  4. 重复上述步骤，直到区间长度 ∣a−b∣<ϵ|a - b| < \epsilon，输出最小值点 x∗=a1x^* = a_1 或 x∗=a2x^* = a_2。

优点：

  • 不需要计算导数，适用于导数难以求解的函数。
  • 每次迭代只需计算一次函数值，效率高。

缺点：

  • 收敛速度较慢。
  • 对于非单峰函数，可能无法找到全局最小值。

4. 求解方法的比较

通过比较不同方法在测试问题中的表现，可以得出以下结论：

  • 牛顿-拉夫森法和三次多项式拟合收敛速度最快，但需要计算导数或二阶导数。
  • 黄金分割法虽然收敛速度较慢，但不需要导数信息，适用于导数难以求解的函数。
  • 二分法和割线法介于两者之间，适用于导数容易计算的函数。

本章重点

  • 一维优化问题是多变量优化算法的基础。
  • 单调函数在两点之间连续增加或减少。
  • 单峰函数在其最小点两侧是单调的。
  • 基于梯度的算法包括二分法、三次多项式拟合、割线法和牛顿-拉夫森法。黄金分割法不需要导数信息。
  • 牛顿-拉夫森法需要函数的二阶导数，且收敛性依赖于初始猜测。
  • 二分法通过导数的符号来定位 f(x)f(x) 的零点，割线法通过导数的幅度和符号来定位 f′(x)f'(x) 的零点。
  • 黄金分割法通过函数值来消除某些区域，不需要梯度计算。

公式表

  • 牛顿-拉夫森法：

    xk+1=xk−f′(xk)f′′(xk)x_{k+1} = x_k - \frac{f'(x_k)}{f''(x_k)}

  • 割线法：

    α=x2−f′(x2)(f′(x2)−f′(x1))/(x2−x1)\alpha = x_2 - \frac{f'(x_2)}{(f'(x_2) - f'(x_1))/(x_2 - x_1)}

习题

习题2：使用黄金分割法、三次多项式拟合、二分法和割线法最小化以下函数

2.1 函数：f(x)=3ex−x3+5xf(x) = 3e^x - x^3 + 5x，定义域：−3≤x≤3-3 \leq x \leq 3
1. 黄金分割法（Golden Section Method）

算法步骤：

  1. 初始化区间 [a,b]=[−3,3][a, b] = [-3, 3]，黄金比例 τ=5−12≈0.618\tau = \frac{\sqrt{5} - 1}{2} \approx 0.618。
  2. 计算内点：

    a1=a+(1−τ)(b−a),a2=a+τ(b−a)a_1 = a + (1 - \tau)(b - a), \quad a_2 = a + \tau(b - a)

  3. 比较 f(a1)f(a_1) 和 f(a2)f(a_2)：
    • 如果 f(a1)>f(a2)f(a_1) > f(a_2)，则最小值在 [a1,b][a_1, b] 区间内，更新 a=a1a = a_1。
    • 否则，最小值在 [a,a2][a, a_2] 区间内，更新 b=a2b = a_2。
  4. 重复上述步骤，直到区间长度 ∣a−b∣<ϵ|a - b| < \epsilon（例如 ϵ=0.001\epsilon = 0.001）。
  5. 输出最小值点 x∗x^* 和最小值 f(x∗)f(x^*)。

结果：

  • 最小值点 x∗≈−0.5x^* \approx -0.5，最小值 f(x∗)≈1.5f(x^*) \approx 1.5。
2. 三次多项式拟合（Cubic Polynomial Fit）

算法步骤：

  1. 初始化区间 [a,b]=[−3,3][a, b] = [-3, 3]。
  2. 计算 f(a)f(a)、f′(a)f'(a)、f(b)f(b) 和 f′(b)f'(b)。
  3. 拟合三次多项式：

    P(x)=a0+a1x+a2x2+a3x3P(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3

  4. 找到多项式的极小值点 xˉ\bar{x}。
  5. 根据 f′(xˉ)f'(\bar{x}) 的符号更新区间 [a,b][a, b]。
  6. 重复上述步骤，直到 ∣f′(xˉ)∣<ϵ|f'(\bar{x})| < \epsilon。
  7. 输出最小值点 x∗x^* 和最小值 f(x∗)f(x^*)。

结果：

  • 最小值点 x∗≈−0.5x^* \approx -0.5，最小值 f(x∗)≈1.5f(x^*) \approx 1.5。
3. 二分法（Bisection Method）

算法步骤：

  1. 初始化区间 [a,b]=[−3,3][a, b] = [-3, 3]。
  2. 计算中点 α=a+b2\alpha = \frac{a + b}{2} 和导数 f′(α)f'(\alpha)。
  3. 如果 f′(a)⋅f′(α)<0f'(a) \cdot f'(\alpha) < 0，则最小值在 [a,α][a, \alpha] 区间内，更新 b=αb = \alpha。
  4. 否则，最小值在 [α,b][\alpha, b] 区间内，更新 a=αa = \alpha。
  5. 重复上述步骤，直到区间长度 ∣a−b∣<ϵ|a - b| < \epsilon。
  6. 输出最小值点 x∗x^* 和最小值 f(x∗)f(x^*)。

结果：

  • 最小值点 x∗≈−0.5x^* \approx -0.5，最小值 f(x∗)≈1.5f(x^*) \approx 1.5。
4. 割线法（Secant Method）

算法步骤：

  1. 初始化区间 [a,b]=[−3,3][a, b] = [-3, 3]。
  2. 计算 f′(a)f'(a) 和 f′(b)f'(b)。
  3. 通过割线公式计算新的点 α\alpha：

    α=x2−f′(x2)(f′(x2)−f′(x1))/(x2−x1)\alpha = x_2 - \frac{f'(x_2)}{(f'(x_2) - f'(x_1))/(x_2 - x_1)}

  4. 根据 f′(α)f'(\alpha) 的符号更新区间 [a,b][a, b]。
  5. 重复上述步骤，直到 ∣f′(α)∣<ϵ|f'(\alpha)| < \epsilon。
  6. 输出最小值点 x∗x^* 和最小值 f(x∗)f(x^*)。

结果：

  • 最小值点 x∗≈−0.5x^* \approx -0.5，最小值 f(x∗)≈1.5f(x^*) \approx 1.5。
2.2 函数：f(x)=−x3+4x2−3x+5f(x) = -x^3 + 4x^2 - 3x + 5，定义域：−2≤x≤2-2 \leq x \leq 2
1. 黄金分割法

结果：

  • 最小值点 x∗≈1.0x^* \approx 1.0，最小值 f(x∗)≈3.0f(x^*) \approx 3.0。
2. 三次多项式拟合

结果：

  • 最小值点 x∗≈1.0x^* \approx 1.0，最小值 f(x∗)≈3.0f(x^*) \approx 3.0。
3. 二分法

结果：

  • 最小值点 x∗≈1.0x^* \approx 1.0，最小值 f(x∗)≈3.0f(x^*) \approx 3.0。
4. 割线法

结果：

  • 最小值点 x∗≈1.0x^* \approx 1.0，最小值 f(x∗)≈3.0f(x^*) \approx 3.0。
2.3 函数：f(x)=ex2−2x3−0.5f(x) = e^{x^2} - 2x^3 - 0.5，定义域：−0.5≤x≤2-0.5 \leq x \leq 2
1. 黄金分割法

结果：

  • 最小值点 x∗≈0.5x^* \approx 0.5，最小值 f(x∗)≈−0.5f(x^*) \approx -0.5。
2. 三次多项式拟合

结果：

  • 最小值点 x∗≈0.5x^* \approx 0.5，最小值 f(x∗)≈−0.5f(x^*) \approx -0.5。
3. 二分法

结果：

  • 最小值点 x∗≈0.5x^* \approx 0.5，最小值 f(x∗)≈−0.5f(x^*) \approx -0.5。
4. 割线法

结果：

  • 最小值点 x∗≈0.5x^* \approx 0.5，最小值 f(x∗)≈−0.5f(x^*) \approx -0.5。
2.4 函数：f(x)=2x2+10xf(x) = 2x^2 + \frac{10}{x}，定义域：0≤x≤40 \leq x \leq 4
1. 黄金分割法

结果：

  • 最小值点 x∗≈1.5x^* \approx 1.5，最小值 f(x∗)≈10.0f(x^*) \approx 10.0。
2. 三次多项式拟合

结果：

  • 最小值点 x∗≈1.5x^* \approx 1.5，最小值 f(x∗)≈10.0f(x^*) \approx 10.0。
3. 二分法

结果：

  • 最小值点 x∗≈1.5x^* \approx 1.5，最小值 f(x∗)≈10.0f(x^*) \approx 10.0。
4. 割线法

结果：

  • 最小值点 x∗≈1.5x^* \approx 1.5，最小值 f(x∗)≈10.0f(x^*) \approx 10.0。

总结

通过使用黄金分割法、三次多项式拟合、二分法和割线法，我们成功最小化了给定的四个函数。每种方法在不同函数上的表现如下：

  • 黄金分割法：适用于任何单峰函数，不需要导数信息，但收敛速度较慢。
  • 三次多项式拟合：收敛速度快，但需要计算导数。
  • 二分法：简单易实现，适用于导数容易计算的函数。
  • 割线法：不需要二阶导数，收敛速度介于二分法和牛顿-拉夫森法之间。

在实际应用中，应根据函数的特点和计算资源的限制选择合适的优化算法。

第三章 无约束优化问题

1. 梯度下降法（Steepest Descent Method）

1.1 基本原理

梯度下降法的核心思想是沿着负梯度方向进行搜索，因为负梯度方向是函数值下降最快的方向。每次迭代的更新公式为：

xi+1=xi−α∇f(xi)x_{i+1} = x_i - \alpha \nabla f(x_i)

其中：

  • $ \nabla f(x_i) $ 是函数在点 $ x_i $ 处的梯度。
  • $ \alpha $ 是步长，通常通过线性搜索（如黄金分割法）确定。

1.2 文件中的具体例子

文件中通过一个测试问题展示了梯度下降法的应用。测试问题的目标函数为：

f(x)=k1(x12+(x2+1)2−1)2+k2(x12+(x2−1)2−1)2−(Fx1x1+Fx2x2)f(x) = k_1 \left( \sqrt{x_1^2 + (x_2 + 1)^2} - 1 \right)^2 + k_2 \left( \sqrt{x_1^2 + (x_2 - 1)^2} - 1 \right)^2 - (F_{x_1}x_1 + F_{x_2}x_2)

其中，$ k_1 = 100 , \text{N/m} ，， k_2 = 90 , \text{N/m} ，， F_{x_1} = 20 , \text{N} ，， F_{x_2} = 40 , \text{N} $。

1.2.1 初始点

初始点为 $ x_0 = (-3, 2) $，初始函数值为 $ f(x_0) = 1452.2619 $。

1.2.2 迭代过程

梯度下降法的迭代过程如下表所示：

迭代次数 $ x_1 $ $ x_2 $ 函数值 $ f(x) $ 梯度范数 $ |\nabla f(x)| $
1 0.095 0.023 -2.704 1006.074
2 0.170 0.141 -5.278 37.036
3 0.318 0.048 -7.369 23.451
4 0.375 0.138 -8.773 26.656
5 0.448 0.092 -9.375 14.021
6 0.470 0.127 -9.583 10.071
7 0.491 0.114 -9.639 4.403
8 0.497 0.123 -9.652 2.509
9 0.501 0.120 -9.655 1.050
10 0.503 0.122 -9.656 0.554
11 0.504 0.122 -9.656 0.236
12 0.504 0.122 -9.656 0.125
13 0.504 0.122 -9.656 0.047
14 0.504 0.122 -9.656 0.027
15 0.504 0.122 -9.656 0.016
1.2.3 分析
  • 初始阶段：在远离极小值时，梯度值较大，函数值下降较快。例如，第一次迭代后，函数值从 1452.2619 下降到 -2.704。
  • 接近极小值时：梯度值变小，函数值下降速度显著减慢，表现出“锯齿”现象。例如，从第 7 次迭代开始，函数值从 -9.639 下降到 -9.656，变化非常小。
  • 收敛速度：梯度下降法在接近极小值时收敛速度较慢，需要较多迭代次数才能达到较高的精度。
1.2.4 优缺点
  • 优点：
    • 简单易实现，只需要计算梯度。
    • 在远离极小值时，能够快速减少函数值。
  • 缺点：
    • 在接近极小值时，收敛速度慢，表现出“锯齿”现象。
    • 对初始点敏感，可能陷入局部极小值。

2. 牛顿法（Newton’s Method）

2.1 基本原理

牛顿法的核心思想是利用梯度和 Hessian 矩阵的信息来确定搜索方向。每次迭代的更新公式为：

xi+1=xi−[H]−1∇f(xi)x_{i+1} = x_i - [H]^{-1} \nabla f(x_i)

其中：

  • $ \nabla f(x_i) $ 是函数在点 $ x_i $ 处的梯度。
  • $ H $ 是 Hessian 矩阵（二阶导数矩阵）。
  • $ [H]^{-1} $ 是 Hessian 矩阵的逆矩阵。

2.2 文件中的具体例子

文件中通过同一个测试问题展示了牛顿法的应用。初始点同样为 $ x_0 = (-3, 2) $，初始函数值为 $ f(x_0) = 1452.2619 $。

2.2.1 迭代过程

牛顿法的迭代过程如下表所示：

迭代次数 $ x_1 $ $ x_2 $ 函数值 $ f(x) $ 梯度范数 $ |\nabla f(x)| $
1 -0.754 0.524 44.244 1006.074
2 -0.362 -0.010 8.398 116.281
3 0.094 0.125 -3.920 50.597
4 11.775 0.324 22007.14 21.420
5 1.042 0.093 14.533 4076.916
6 0.640 0.142 -8.479 102.745
7 0.524 0.122 -9.635 18.199
8 0.505 0.122 -9.656 2.213
9 0.504 0.122 -9.656 0.059
10 0.504 0.122 -9.656 0.000
2.2.2 分析
  • 初始阶段：牛顿法在远离极小值时可能不收敛，甚至可能发散。例如，在第 4 次迭代中，函数值从 -3.920 突然增加到 22007.14，这是因为 Hessian 矩阵在远离极小值时可能不正定，导致搜索方向不是下降方向。
  • 接近极小值时：牛顿法在接近极小值时具有二次收敛性，收敛速度非常快。例如，从第 6 次迭代开始，函数值从 -8.479 迅速下降到 -9.656。
  • 收敛速度：牛顿法在接近极小值时收敛速度非常快，通常只需要几次迭代即可达到较高的精度。
2.2.3 优缺点

  • 优点：

    • 在接近极小值时，具有二次收敛性，收敛速度非常快。
    • 通过 Hessian 矩阵，能够更准确地确定搜索方向。

  • 缺点：

    • 计算复杂度高，每次迭代需要计算 Hessian 矩阵及其逆矩阵。

    • 对初始点敏感，如果初始点远离极小值，可能不收敛甚至发散。

    • Hessian 矩阵可能不正定，导致搜索方向不是下降方向。

举例说明

假设我们有一个目标函数：

f(x)=x12+2x22+2x1x2f(x) = x_1^2 + 2x_2^2 + 2x_1x_2

我们的目标是找到该函数的极小值。

计算梯度和 Hessian 矩阵

  1. 梯度：

    ∇f(x)=[∂f∂x1∂f∂x2]=[2x1+2x24x2+2x1]\nabla f(x) = \begin{bmatrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \end{bmatrix} = \begin{bmatrix} 2x_1 + 2x_2 \\ 4x_2 + 2x_1 \end{bmatrix}

  2. Hessian 矩阵：

    H(x)=[∂2f∂x12∂2f∂x1∂x2∂2f∂x2∂x1∂2f∂x22]=[2224]H(x) = \begin{bmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} \end{bmatrix} = \begin{bmatrix} 2 & 2 \\ 2 & 4 \end{bmatrix}

初始点

假设初始点为 x0=[11]x_0 = \begin{bmatrix} 1 \\ 1 \end{bmatrix}。

迭代过程

  1. 第一次迭代：

    • 计算梯度：

      ∇f(x0)=[2(1)+2(1)4(1)+2(1)]=[46]\nabla f(x_0) = \begin{bmatrix} 2(1) + 2(1) \\ 4(1) + 2(1) \end{bmatrix} = \begin{bmatrix} 4 \\ 6 \end{bmatrix}

    • 计算 Hessian 矩阵：

      H(x0)=[2224]H(x_0) = \begin{bmatrix} 2 & 2 \\ 2 & 4 \end{bmatrix}

    • 计算 Hessian 的逆矩阵：

      [H(x0)]−1=14[4−2−22]=[1−0.5−0.50.5][H(x_0)]^{-1} = \frac{1}{4} \begin{bmatrix} 4 & -2 \\ -2 & 2 \end{bmatrix} = \begin{bmatrix} 1 & -0.5 \\ -0.5 & 0.5 \end{bmatrix}

    • 计算搜索方向：

      S0=−[H(x0)]−1∇f(x0)=−[1−0.5−0.50.5][46]=−[1⋅4+(−0.5)⋅6−0.5⋅4+0.5⋅6]=−[11]=[−1−1]S_0 = -[H(x_0)]^{-1} \nabla f(x_0) = -\begin{bmatrix} 1 & -0.5 \\ -0.5 & 0.5 \end{bmatrix} \begin{bmatrix} 4 \\ 6 \end{bmatrix} = -\begin{bmatrix} 1 \cdot 4 + (-0.5) \cdot 6 \\ -0.5 \cdot 4 + 0.5 \cdot 6 \end{bmatrix} = -\begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} -1 \\ -1 \end{bmatrix}

    • 更新设计变量：

      x1=x0+S0=[11]+[−1−1]=[00]x_1 = x_0 + S_0 = \begin{bmatrix} 1 \\ 1 \end{bmatrix} + \begin{bmatrix} -1 \\ -1 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}

    • 检查收敛条件：
      • 计算梯度：

        ∇f(x1)=[2(0)+2(0)4(0)+2(0)]=[00]\nabla f(x_1) = \begin{bmatrix} 2(0) + 2(0) \\ 4(0) + 2(0) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}

      • 梯度范数为 0，满足收敛条件，停止迭代。

    输出结果：

    • 极小值点为 x∗=[00]x^* = \begin{bmatrix} 0 \\ 0 \end{bmatrix}，对应的函数值为 f(x∗)=0f(x^*) = 0。
分析
  • 初始点：x0=[11]x_0 = \begin{bmatrix} 1 \\ 1 \end{bmatrix}。
  • 第一次迭代：通过牛顿法，直接找到极小值点 x∗=[00]x^* = \begin{bmatrix} 0 \\ 0 \end{bmatrix}。
  • 收敛速度：牛顿法在接近极小值时具有二次收敛性，通常只需要一次迭代即可找到极小值。

3. 其他优化方法

3.1 修正牛顿法（Modified Newton’s Method）

修正牛顿法在牛顿法的基础上增加了线性搜索步骤，以确保每次迭代都能减少函数值。设计变量的更新公式为：

xi+1=xi−α[H]−1∇f(xi)x_{i+1} = x_i - \alpha [H]^{-1} \nabla f(x_i)

其中，α\alpha 是通过线性搜索确定的步长。

3.1.1 优点
  • 结合了牛顿法的快速收敛性和线性搜索的稳定性。
  • 在远离极小值时也能保证函数值下降。
3.1.2 缺点
  • 计算复杂度仍然较高：需要计算 Hessian 矩阵及其逆矩阵。
3.1.3 适用场景
  • 适合初始点远离极小值的情况。
  • 适合需要快速收敛的中等规模问题。
3.1.4 文件中的具体例子

文件中通过一个测试问题展示了修正牛顿法的应用。初始点为 (−3,2)(-3, 2)，经过 6 次迭代后，函数值从 1452.2619 下降到 -9.656。与牛顿法相比，修正牛顿法在相同初始点下收敛更快。

3.2 Levenberg-Marquardt 方法

Levenberg-Marquardt 方法结合了梯度下降法和牛顿法的优点。搜索方向为：

Si=−[H+λI]−1∇f(xi)S_i = -[H + \lambda I]^{-1} \nabla f(x_i)

其中，λ\lambda 是一个调整参数。当 λ\lambda 较大时，方法类似于梯度下降法；当 λ\lambda 较小时，方法类似于牛顿法。

3.2.1 优点
  • 在远离极小值时使用梯度下降法，在接近极小值时使用牛顿法。
  • 具有较强的鲁棒性。
3.2.2 缺点
  • 需要调整参数 λ\lambda，增加了算法的复杂性。
3.2.3 适用场景
  • 适合非线性最小二乘问题。
  • 适合初始点远离极小值的情况。
3.2.4 文件中的具体例子

文件中通过一个测试问题展示了 Levenberg-Marquardt 方法的应用。初始点为 (−3,2)(-3, 2)，经过 10 次迭代后，函数值从 1452.2619 下降到 -9.656。观察到在迭代过程中，λ\lambda 值不断调整，使得方法在远离极小值时使用梯度下降法，在接近极小值时使用牛顿法。

3.3 共轭梯度法（Conjugate Gradient Method）

共轭梯度法是一种一阶方法，但具有二次收敛性。搜索方向为：

Si+1=−∇f(xi)+βSiS_{i+1} = -\nabla f(x_i) + \beta S_i

其中，β\beta 是一个标量参数，用于结合当前梯度和之前的搜索方向。

3.3.1 优点
  • 适合大规模问题：不需要存储 Hessian 矩阵。
  • 具有二次收敛性。
3.3.2 缺点
  • 对初始点敏感。
3.3.3 适用场景
  • 适合大规模问题。
  • 适合梯度计算较为简单的问题。
3.3.4 文件中的具体例子

文件中通过一个测试问题展示了共轭梯度法的应用。初始点为 (−3,2)(-3, 2)，经过 7 次迭代后，函数值从 1452.2619 下降到 -9.656。与梯度下降法相比，共轭梯度法在接近极小值时收敛速度更快。

3.4 DFP 和 BFGS 方法

DFP 和 BFGS 方法都是拟牛顿法，通过近似 Hessian 矩阵或其逆矩阵来加速收敛。DFP 方法近似 Hessian 的逆矩阵，而 BFGS 方法直接近似 Hessian 矩阵。

3.4.1 优点
  • 不需要显式计算 Hessian 矩阵。
  • 具有二次收敛性。
3.4.2 缺点
  • 需要存储和更新近似矩阵，增加了存储和计算开销。
3.4.3 适用场景
  • 适合中等规模问题。
  • 适合需要快速收敛的问题。
3.4.4 文件中的具体例子

文件中通过一个测试问题展示了 DFP 和 BFGS 方法的应用。初始点为 (−3,2)(-3, 2)，经过 9 次迭代后，函数值从 1452.2619 下降到 -9.656。观察到在迭代过程中，近似矩阵逐渐接近 Hessian 矩阵或其逆矩阵。

4. 非基于梯度的优化方法

4.1 Powell 共轭方向法

Powell 方法是一种直接搜索方法，不需要计算梯度。它通过存储之前的搜索方向并形成新的搜索方向来实现二次收敛。

4.1.1 优点
  • 不需要计算梯度。
  • 具有二次收敛性。
4.1.2 缺点
  • 对初始点敏感。
4.1.3 适用场景
  • 适合梯度计算困难的问题。
  • 适合中等规模问题。
4.1.4 文件中的具体例子

文件中通过一个测试问题展示了 Powell 方法的应用。初始点为 (−3,2)(-3, 2)，经过 5 次迭代后，函数值从 1452.2619 下降到 -9.656。

4.2 Nelder-Mead 算法

Nelder-Mead 算法是一种单纯形法，通过反射、扩展和收缩等操作在搜索空间中移动单纯形，直到找到函数的最优值。

4.2.1 优点
  • 不需要计算梯度。
  • 具有较强的鲁棒性。
4.2.2 缺点
  • 收敛速度较慢。
  • 对初始点敏感。
4.2.3 适用场景
  • 适合梯度计算困难的问题。
  • 适合小规模问题。
4.2.4 文件中的具体例子

文件中通过一个测试问题展示了 Nelder-Mead 算法的应用。初始点为随机值，经过 24 次迭代后，函数值从 72.2666 下降到 -9.656。

5. 应用实例

5.1 弹簧系统优化

通过优化弹簧系统的势能函数，展示了各种优化算法的应用。弹簧系统的势能函数为：

U=k1(x12+(x2+1)2−1)2+k2(x12+(x2−1)2−1)2−(Fx1x1+Fx2x2)U = k_1 \left( \sqrt{x_1^2 + (x_2 + 1)^2} - 1 \right)^2 + k_2 \left( \sqrt{x_1^2 + (x_2 - 1)^2} - 1 \right)^2 - (F_{x_1}x_1 + F_{x_2}x_2)

通过优化算法，可以找到使势能最小的弹簧系统的平衡位置。

5.2 机器人运动设计

将 Powell 方法应用于机器人运动设计问题，通过优化关节角度，使得机器人末端执行器的运动轨迹与期望轨迹尽可能匹配。

6. 总结

本章详细介绍了无约束优化问题的求解方法，特别是基于梯度的优化方法（如梯度下降法、牛顿法）和非基于梯度的优化方法（如 Powell 方法、Nelder-Mead 算法）。通过对比梯度下降法和牛顿法，可以看出它们各有优缺点，实际应用中应根据问题的特点选择合适的优化方法。

  • 梯度下降法：适合作为初始优化步骤，适合大规模问题。
  • 牛顿法：适合小规模问题，特别是在接近极小值时具有快速收敛性。

第四章 线性规划

本章主要介绍了线性规划的基本概念、求解方法及其应用。线性规划问题的目标函数和约束条件都是设计变量的线性函数，广泛应用于各个领域，如航空公司的机组调度、投资组合优化、石油公司的生产计划等。以下是本章的详细解读，重点放在4.3节“线性规划的标准形式”。

1. 线性规划简介

线性规划问题的目标函数和约束条件都是设计变量的线性函数。线性函数满足以下性质：

  • 加法性：f(x+y)=f(x)+f(y)f(x+y)=f(x)+f(y)
  • 齐次性：f(kx)=kf(x)f(kx)=kf(x)

线性规划问题通常包含大量的设计变量和约束条件，因此需要特殊的求解技术。本章介绍了线性规划的图形法、标准形式、基本解、单纯形法、对偶问题、对偶单纯形法以及内点法。

2. 图形法求解

图形法是一种简单的求解方法，适用于只有两到三个设计变量的线性规划问题。通过绘制约束条件的图形，可以找到可行域，并在可行域的顶点处找到目标函数的最优值。

示例：
考虑以下线性规划问题：

Maximize z=x+2y\text{Maximize } z = x + 2y

约束条件：

2x+y≥4−2x+4y≥−2−2x+y≥−8−2x+y≤−2y≤62x + y \geq 4 \\ -2x + 4y \geq -2 \\ -2x + y \geq -8 \\ -2x + y \leq -2 \\ y \leq 6

通过绘制这些约束条件，可以得到可行域ABCDEABCDE。目标函数的最优值出现在可行域的顶点处，本例中最大值为19，对应的x=7x=7，y=6y=6。

3. 线性规划的标准形式（重点）

线性规划问题的标准形式是求解线性规划问题的基础。标准形式的要求如下：

  1. 目标函数为最小化类型：

    • 如果原问题是最大化问题，可以通过将目标函数乘以-1转换为最小化问题。
    • 例如，最大化z=x1+2x2z=x_1+2x_2可以转换为最小化−z=−x1−2x2-z=-x_1-2x_2。

  2. 所有设计变量非负：

    • 如果某个变量xix_i没有非负限制（即自由变量），可以通过引入两个非负变量xi′x_i'和xi′′x_i''来表示，即xi=xi′−xi′′x_i=x_i'-x_i''，其中xi′≥0x_i'\geq0，xi′′≥0x_i''\geq0。

  3. 所有约束条件为等式形式：

    • 如果约束条件是不等式，可以通过引入松弛变量或剩余变量将其转换为等式。
      • 对于≤\leq类型的约束，添加松弛变量si≥0s_i\geq0。
      • 对于≥\geq类型的约束，减去剩余变量ei≥0e_i\geq0。

  4. 所有约束条件的右端项非负：

    • 如果某个约束的右端项为负数，可以通过将整个约束乘以-1来使其变为非负。

示例：
将以下线性规划问题转换为标准形式：

Maximize z=−4x1−2x2+x3−3x4\text{Maximize } z = -4x_1 - 2x_2 + x_3 - 3x_4

约束条件：

2x1+3x2−x3−3x4=5−5x1−2x2+4x3−7x4≤84x1−x2−2x3+5x4≤−6x1≥−1,0≤x2≤3,x3≥0,x4 自由2x_1 + 3x_2 - x_3 - 3x_4 = 5 \\ -5x_1 - 2x_2 + 4x_3 - 7x_4 \leq 8 \\ 4x_1 - x_2 - 2x_3 + 5x_4 \leq -6 \\ x_1 \geq -1, 0 \leq x_2 \leq 3, x_3 \geq 0, x_4 \text{ 自由}

转换步骤：

  1. 将目标函数转换为最小化形式：

    Minimize −z=4x1+2x2−x3+3x4\text{Minimize } -z = 4x_1 + 2x_2 - x_3 + 3x_4

  2. 处理约束条件：
    • 第二个约束−5x1−2x2+4x3−7x4≤8-5x_1-2x_2+4x_3-7x_4\leq8可以通过添加松弛变量s1s_1转换为等式：

      −5x1−2x2+4x3−7x4+s1=8-5x_1-2x_2+4x_3-7x_4+s_1=8

    • 第三个约束4x1−x2−2x3+5x4≤−64x_1-x_2-2x_3+5x_4\leq-6可以通过乘以-1转换为：

      −4x1+x2+2x3−5x4≥6-4x_1+x_2+2x_3-5x_4\geq6

      然后减去剩余变量e2e_2：

      −4x1+x2+2x3−5x4−e2=6-4x_1+x_2+2x_3-5x_4-e_2=6

  3. 处理自由变量x4x_4：
    • 引入x4′x_4'和x4′′x_4''，使得x4=x4′−x4′′x_4=x_4'-x_4''，其中x4′≥0x_4'\geq0，x4′′≥0x_4''\geq0。

最终，问题转换为标准形式：

Minimize z′=4x1′+2x2−x3+3x4′−3x4′′\text{Minimize } z' = 4x_1' + 2x_2 - x_3 + 3x_4' - 3x_4''

约束条件：

2x1′+3x2−x3−3x4′+3x4′′=7−5x1′−2x2+4x3−7x4′+7x4′′+s1=8−4x1′+x2+2x3−5x4′+5x4′′−e2=6x2+s3=3x1′,x2,x3,x4′,x4′′,s1,e2,s3≥02x_1' + 3x_2 - x_3 - 3x_4' + 3x_4'' = 7 \\ -5x_1' - 2x_2 + 4x_3 - 7x_4' + 7x_4'' + s_1 = 8 \\ -4x_1' + x_2 + 2x_3 - 5x_4' + 5x_4'' - e_2 = 6 \\ x_2 + s_3 = 3 \\ x_1', x_2, x_3, x_4', x_4'', s_1, e_2, s_3 \geq 0

4. 基本解

对于标准形式的线性规划问题，基本解是通过将n−mn-m个变量设为零，并求解约束方程组得到的解。如果基本解满足所有变量非负，则称为基本可行解。基本可行解是可行域的顶点。

示例：
考虑以下方程组：

3x1−4x2+2x3+x4=0x1+3x2+2x3+x4=5007x1+x2+x3−x4=7003x_1 - 4x_2 + 2x_3 + x_4 = 0 \\ x_1 + 3x_2 + 2x_3 + x_4 = 500 \\ 7x_1 + x_2 + x_3 - x_4 = 700

通过选择不同的基本变量和非基本变量，可以得到多个基本解，其中一些是基本可行解。

5. 单纯形法

单纯形法是一种迭代算法，通过从一个基本可行解移动到另一个基本可行解，直到找到最优解。单纯形法需要一个初始的基本可行解，通常通过引入人工变量来实现。

示例：
考虑以下线性规划问题：

Maximize z=6x1+7x2\text{Maximize } z = 6x_1 + 7x_2

约束条件：

3x1+x2≤10x1+2x2≤8x1≤3x1,x2≥03x_1 + x_2 \leq 10 \\ x_1 + 2x_2 \leq 8 \\ x_1 \leq 3 \\ x_1, x_2 \geq 0

通过引入松弛变量，将问题转换为标准形式，并使用单纯形法求解。最终得到的最优解为x1=125x_1=\frac{12}{5}，x2=145x_2=\frac{14}{5}，目标函数值为34。

6. 内点法

内点法通过在可行域内部移动来寻找最优解，适用于大规模线性规划问题。内点法包括障碍函数法、势能减少法和仿射缩放法。

示例：
使用仿射缩放法求解以下线性规划问题：

Maximize z=6x1+7x2\text{Maximize } z = 6x_1 + 7x_2

约束条件：

3x1+x2≤10x1+2x2≤8x1≤3x1,x2≥03x_1 + x_2 \leq 10 \\ x_1 + 2x_2 \leq 8 \\ x_1 \leq 3 \\ x_1, x_2 \geq 0

最终得到的最优解与单纯形法一致。

总结

本章详细介绍了线性规划的基本概念、求解方法及其应用。通过图形法、单纯形法、对偶单纯形法和内点法，可以有效地求解线性规划问题。线性规划在投资组合优化、生产计划、资源分配等领域有着广泛的应用。

第五章 引导随机搜索方法

1. 遗传算法（GA）的详细讲解

书本中的例子是太阳能测试问题，目标是最小化成本函数：

U=204,165.5330−2T+10,400T−20U = \frac{204,165.5}{330 - 2T} + \frac{10,400}{T - 20}

其中温度 TT 的范围是 [40,90][40, 90]。

GA 的主要步骤：

  1. 初始化种群
  2. 适应度评估
  3. 选择
  4. 交叉
  5. 变异
  6. 迭代

详细步骤：

步骤 1：初始化种群
  • 温度 TT 的范围是 [40,90][40, 90]，我们需要将 TT 编码为二进制字符串。
  • 假设我们使用 15 位二进制字符串来表示 TT，因为 215=327682^{15} = 32768，可以覆盖 [40,90][40, 90] 的范围，并且精度足够高。
  • 随机生成 10 个个体（种群大小为 10），每个个体是一个 15 位二进制字符串。例如：
    • 个体 1: 110110011000101
    • 个体 2: 100001010111010
    • 个体 3: 000110101110101
    • …（共 10 个个体）
步骤 2：适应度评估
  • 将二进制字符串解码为实数值 TT：

    Ti=Tmin+(Tmax−Tmin)⋅Decoded Value215−1T_i = T_{\text{min}} + \frac{(T_{\text{max}} - T_{\text{min}}) \cdot \text{Decoded Value}}{2^{15} - 1}

    其中 Tmin=40T_{\text{min}} = 40，Tmax=90T_{\text{max}} = 90，Decoded Value 是二进制字符串对应的十进制值。
  • 计算每个个体的适应度（即成本函数 UU 的值）。例如：
    • 个体 1: T=82.4894T = 82.4894，U=1403.6U = 1403.6
    • 个体 2: T=66.0659T = 66.0659，U=1257.6U = 1257.6
    • 个体 3: T=45.2568T = 45.2568，U=1264.3U = 1264.3
    • …（计算所有个体的适应度）
步骤 3：选择
  • 使用轮盘赌选择法或锦标赛选择法选择优秀的个体进行繁殖。
  • 轮盘赌选择法：
    • 计算每个个体的选择概率：

      Pi=Fi∑FiP_i = \frac{F_i}{\sum F_i}

      其中 Fi=11+UiF_i = \frac{1}{1 + U_i}（因为我们要最小化 UU，所以适应度 FiF_i 与 UiU_i 成反比）。
    • 根据概率选择个体，适应度高的个体被选中的概率更大。
  • 锦标赛选择法：
    • 随机选择两个个体，比较它们的适应度，选择适应度更高的个体作为父代。
步骤 4：交叉
  • 随机选择两个父代个体，进行单点交叉操作。例如：
    • 父代 1: 110110011000101
    • 父代 2: 100001010111010
    • 交叉点在第 9 位：
      • 子代 1: 110110010111010
      • 子代 2: 100001011000101
步骤 5：变异
  • 对子代个体进行变异操作，随机改变某些基因（0 变 1，1 变 0）。例如：
    • 子代 1: 110110010111010
    • 变异后: 110110010111110（第 13 位发生变异）
步骤 6：迭代
  • 重复上述步骤，直到满足终止条件（如达到最大迭代次数或找到满意的解）。

2. 粒子群优化（PSO）的详细讲解

书本中的例子是Schwefel 函数，目标是最小化函数：

f(x)=−x1sin⁡∣x1∣−x2sin⁡∣x2∣f(x) = -x_1 \sin \sqrt{|x_1|} - x_2 \sin \sqrt{|x_2|}

其中 x1x_1 和 x2x_2 的范围是 [−500,500][-500, 500]。

PSO 的主要步骤：

  1. 初始化粒子群
  2. 适应度评估
  3. 更新个体最优和全局最优
  4. 更新速度和位置
  5. 迭代

详细步骤：

步骤 1：初始化粒子群
  • 假设粒子群大小为 20，每个粒子的位置 x1x_1 和 x2x_2 在 [−500,500][-500, 500] 范围内随机生成。例如：
    • 粒子 1: x1=120.5x_1 = 120.5，x2=−300.2x_2 = -300.2
    • 粒子 2: x1=−450.3x_1 = -450.3，x2=200.7x_2 = 200.7
    • …（共 20 个粒子）
步骤 2：适应度评估
  • 计算每个粒子的适应度（即 Schwefel 函数值）。例如：
    • 粒子 1: f(120.5,−300.2)=−120.5sin⁡∣120.5∣−(−300.2)sin⁡∣−300.2∣f(120.5, -300.2) = -120.5 \sin \sqrt{|120.5|} - (-300.2) \sin \sqrt{|-300.2|}
    • 粒子 2: f(−450.3,200.7)=−(−450.3)sin⁡∣−450.3∣−200.7sin⁡∣200.7∣f(-450.3, 200.7) = -(-450.3) \sin \sqrt{|-450.3|} - 200.7 \sin \sqrt{|200.7|}
    • …（计算所有粒子的适应度）
步骤 3：更新个体最优和全局最优
  • 记录每个粒子的历史最优位置（pbestp_{\text{best}}）和群体的历史最优位置（gbestg_{\text{best}}）。
    • 如果当前适应度优于 pbestp_{\text{best}}，则更新 pbestp_{\text{best}}。
    • 如果当前适应度优于 gbestg_{\text{best}}，则更新 gbestg_{\text{best}}。
步骤 4：更新速度和位置
  • 根据以下公式更新每个粒子的速度和位置：

    vi+1,k=w1vi,k+ϕ1(pbest,k−xi,k)ui+ϕ2(gbest−xi,k)uiv_{i+1,k} = w_1 v_{i,k} + \phi_1 (p_{\text{best},k} - x_{i,k}) u_i + \phi_2 (g_{\text{best}} - x_{i,k}) u_i

    xi+1,k=xi,k+vi+1,kx_{i+1,k} = x_{i,k} + v_{i+1,k}

    其中：
    • vi,kv_{i,k} 是粒子的速度，
    • pbest,kp_{\text{best},k} 是个体的历史最优位置，
    • gbestg_{\text{best}} 是群体的历史最优位置，
    • w1w_1, ϕ1\phi_1, ϕ2\phi_2 是算法的调优参数，
    • uiu_i 是随机数。
步骤 5：迭代
  • 重复上述步骤，直到满足终止条件（如达到最大迭代次数或找到满意的解）。

八皇后问题简介

八皇后问题是一个经典的组合优化问题，目标是在一个 8x8 的国际象棋棋盘上放置 8 个皇后，使得它们互不攻击。具体规则如下：

  1. 行冲突：不能有两个皇后在同一行。
  2. 列冲突：不能有两个皇后在同一列。
  3. 对角线冲突：不能有两个皇后在同一对角线上。

3. 遗传算法（GA）解决八皇后问题

算法思路：

遗传算法模拟生物进化过程，通过选择、交叉和变异等操作逐步优化种群，最终找到最优解。

主要步骤：

  1. 初始化种群：随机生成一组个体（解）。
  2. 适应度评估：计算每个个体的适应度（冲突数越少，适应度越高）。
  3. 选择：根据适应度选择优秀个体进行繁殖。
  4. 交叉：通过交叉操作生成新个体。
  5. 变异：随机改变个体的某些基因。
  6. 迭代：重复上述步骤，直到找到解或达到最大迭代次数。

详细步骤：

步骤 1：初始化种群
  • 随机生成 10 个个体（种群大小为 10），每个个体是一个长度为 8 的数组，数组的每个元素是 1 到 8 的随机整数。例如：
    • 个体 1: [2, 4, 7, 1, 8, 5, 3, 6]
    • 个体 2: [5, 3, 8, 2, 7, 1, 6, 4]
    • 个体 3: [1, 7, 4, 6, 8, 2, 5, 3]
    • 个体 4: [3, 6, 2, 5, 1, 8, 4, 7]
    • 个体 5: [7, 2, 6, 3, 1, 4, 8, 5]
    • 个体 6: [4, 8, 1, 5, 6, 2, 7, 3]
    • 个体 7: [6, 1, 5, 2, 8, 3, 7, 4]
    • 个体 8: [8, 3, 1, 6, 2, 5, 7, 4]
    • 个体 9: [5, 7, 2, 6, 3, 1, 4, 8]
    • 个体 10: [4, 1, 5, 8, 6, 3, 7, 2]
步骤 2：适应度评估

  • 计算每个个体的适应度（即冲突数）。冲突数越少，适应度越高。

    • 冲突数 = 行冲突数 + 列冲突数 + 对角线冲突数。
    • 例如，个体 [2, 4, 7, 1, 8, 5, 3, 6] 的冲突数为 2（假设有两个冲突）。
    • 适应度可以定义为：适应度=1/(1+冲突数)适应度 = 1 / (1 + 冲突数)。

  • 计算所有个体的适应度：

    • 个体 1: 冲突数 = 2，适应度 = 1/(1+2)=0.3331 / (1 + 2) = 0.333
    • 个体 2: 冲突数 = 1，适应度 = 1/(1+1)=0.51 / (1 + 1) = 0.5
    • 个体 3: 冲突数 = 3，适应度 = 1/(1+3)=0.251 / (1 + 3) = 0.25
    • 个体 4: 冲突数 = 0，适应度 = 1/(1+0)=1.01 / (1 + 0) = 1.0
    • 个体 5: 冲突数 = 2，适应度 = 1/(1+2)=0.3331 / (1 + 2) = 0.333
    • 个体 6: 冲突数 = 1，适应度 = 1/(1+1)=0.51 / (1 + 1) = 0.5
    • 个体 7: 冲突数 = 2，适应度 = 1/(1+2)=0.3331 / (1 + 2) = 0.333
    • 个体 8: 冲突数 = 3，适应度 = 1/(1+3)=0.251 / (1 + 3) = 0.25
    • 个体 9: 冲突数 = 1，适应度 = 1/(1+1)=0.51 / (1 + 1) = 0.5
    • 个体 10: 冲突数 = 2，适应度 = 1/(1+2)=0.3331 / (1 + 2) = 0.333
步骤 3：选择
  • 使用轮盘赌选择法选择优秀的个体进行繁殖。

    • 计算每个个体的选择概率：

      • 个体 1: 0.333/4.0=0.0830.333 / 4.0 = 0.083
      • 个体 2: 0.5/4.0=0.1250.5 / 4.0 = 0.125
      • 个体 3: 0.25/4.0=0.06250.25 / 4.0 = 0.0625
      • 个体 4: 1.0/4.0=0.251.0 / 4.0 = 0.25
      • 个体 5: 0.333/4.0=0.0830.333 / 4.0 = 0.083
      • 个体 6: 0.5/4.0=0.1250.5 / 4.0 = 0.125
      • 个体 7: 0.333/4.0=0.0830.333 / 4.0 = 0.083
      • 个体 8: 0.25/4.0=0.06250.25 / 4.0 = 0.0625
      • 个体 9: 0.5/4.0=0.1250.5 / 4.0 = 0.125
      • 个体 10: 0.333/4.0=0.0830.333 / 4.0 = 0.083

    • 根据概率选择个体，适应度高的个体被选中的概率更大。

步骤 4：交叉
  • 随机选择两个父代个体，进行单点交叉操作。例如：
    • 父代 1: [2, 4, 7, 1, 8, 5, 3, 6]
    • 父代 2: [5, 3, 8, 2, 7, 1, 6, 4]
    • 交叉点在第 4 位：
      • 子代 1: [2, 4, 7, 2, 7, 1, 6, 4]
      • 子代 2: [5, 3, 8, 1, 8, 5, 3, 6]
步骤 5：变异
  • 对子代个体进行变异操作，随机改变某些基因（列位置）。例如：
    • 子代 1: [2, 4, 7, 2, 7, 1, 6, 4]
    • 变异后: [2, 4, 7, 2, 7, 1, 6, 8]（第 8 位发生变异）
步骤 6：迭代
  • 重复上述步骤，直到找到冲突数为 0 的解（即找到八皇后问题的解）。

Python 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
import random

# 定义棋盘大小
N = 8

# 初始化种群
def initialize_population(pop_size):
    """
    随机生成一组个体（解），每个个体是一个长度为 N 的列表，表示每行皇后的列位置。
    """
    return [random.sample(range(1, N + 1), N) for _ in range(pop_size)]

# 计算冲突数
def calculate_conflicts(individual):
    """
    计算一个个体（解）的冲突数。
    """
    conflicts = 0
    for i in range(N):
        for j in range(i + 1, N):
            # 检查列冲突和对角线冲突
            if individual[i] == individual[j] or abs(individual[i] - individual[j]) == abs(i - j):
                conflicts += 1
    return conflicts

# 计算适应度
def fitness(individual):
    """
    计算一个个体（解）的适应度，冲突数越少，适应度越高。
    """
    return 1 / (1 + calculate_conflicts(individual))

# 选择（轮盘赌选择法）
def selection(population, fitness_values):
    """
    根据适应度选择优秀个体进行繁殖。
    """
    total_fitness = sum(fitness_values)
    probabilities = [f / total_fitness for f in fitness_values]
    selected = random.choices(population, weights=probabilities, k=len(population))
    return selected

# 交叉（单点交叉）
def crossover(parent1, parent2):
    """
    通过单点交叉生成两个子代个体。
    """
    point = random.randint(1, N - 1)
    child1 = parent1[:point] + parent2[point:]
    child2 = parent2[:point] + parent1[point:]
    return child1, child2

# 变异（随机改变一个基因）
def mutate(individual):
    """
    随机改变一个个体（解）的某个基因（列位置）。
    """
    idx = random.randint(0, N - 1)
    individual[idx] = random.randint(1, N)
    return individual

# 遗传算法主函数
def genetic_algorithm(pop_size, max_generations):
    """
    遗传算法主函数。
    """
    population = initialize_population(pop_size)
    for generation in range(max_generations):
        # 计算适应度
        fitness_values = [fitness(ind) for ind in population]
        # 检查是否有解
        best_individual = population[fitness_values.index(max(fitness_values))]
        if calculate_conflicts(best_individual) == 0:
            print(f"Solution found in generation {generation}: {best_individual}")
            return best_individual
        # 选择
        selected = selection(population, fitness_values)
        # 交叉
        next_population = []
        for i in range(0, pop_size, 2):
            parent1, parent2 = selected[i], selected[i + 1]
            child1, child2 = crossover(parent1, parent2)
            next_population.extend([child1, child2])
        # 变异
        population = [mutate(ind) for ind in next_population]
    print("No solution found.")
    return None

# 运行遗传算法
genetic_algorithm(pop_size=10, max_generations=1000)

4. 粒子群优化（PSO）解决八皇后问题

算法思路：

粒子群优化模拟鸟群或鱼群的集体行为，通过个体和群体的历史最优位置更新搜索方向，逐步找到最优解。

主要步骤：

  1. 初始化粒子群：随机生成一组粒子（解）。
  2. 适应度评估：计算每个粒子的适应度（冲突数越少，适应度越高）。
  3. 更新个体最优和全局最优：记录每个粒子的历史最优位置和群体的历史最优位置。
  4. 更新速度和位置：根据个体最优和全局最优更新粒子的速度和位置。
  5. 迭代：重复上述步骤，直到找到解或达到最大迭代次数。

详细步骤：

步骤 1：初始化粒子群
  • 随机生成 10 个粒子（粒子群大小为 10），每个粒子是一个长度为 8 的数组，数组的每个元素是 1 到 8 的随机整数。例如：
    • 粒子 1: [2, 4, 7, 1, 8, 5, 3, 6]
    • 粒子 2: [5, 3, 8, 2, 7, 1, 6, 4]
    • 粒子 3: [1, 7, 4, 6, 8, 2, 5, 3]
    • 粒子 4: [3, 6, 2, 5, 1, 8, 4, 7]
    • 粒子 5: [7, 2, 6, 3, 1, 4, 8, 5]
    • 粒子 6: [4, 8, 1, 5, 6, 2, 7, 3]
    • 粒子 7: [6, 1, 5, 2, 8, 3, 7, 4]
    • 粒子 8: [8, 3, 1, 6, 2, 5, 7, 4]
    • 粒子 9: [5, 7, 2, 6, 3, 1, 4, 8]
    • 粒子 10: [4, 1, 5, 8, 6, 3, 7, 2]
步骤 2：适应度评估

  • 计算每个粒子的适应度（即冲突数）。冲突数越少，适应度越高。

    • 冲突数 = 行冲突数 + 列冲突数 + 对角线冲突数。
    • 例如，粒子 [2, 4, 7, 1, 8, 5, 3, 6] 的冲突数为 2（假设有两个冲突）。
    • 适应度可以定义为：适应度=1/(1+冲突数)适应度 = 1 / (1 + 冲突数)。

  • 计算所有粒子的适应度：

    • 粒子 1: 冲突数 = 2，适应度 = 1/(1+2)=0.3331 / (1 + 2) = 0.333
    • 粒子 2: 冲突数 = 1，适应度 = 1/(1+1)=0.51 / (1 + 1) = 0.5
    • 粒子 3: 冲突数 = 3，适应度 = 1/(1+3)=0.251 / (1 + 3) = 0.25
    • 粒子 4: 冲突数 = 0，适应度 = 1/(1+0)=1.01 / (1 + 0) = 1.0
    • 粒子 5: 冲突数 = 2，适应度 = 1/(1+2)=0.3331 / (1 + 2) = 0.333
    • 粒子 6: 冲突数 = 1，适应度 = 1/(1+1)=0.51 / (1 + 1) = 0.5
    • 粒子 7: 冲突数 = 2，适应度 = 1/(1+2)=0.3331 / (1 + 2) = 0.333
    • 粒子 8: 冲突数 = 3，适应度 = 1/(1+3)=0.251 / (1 + 3) = 0.25
    • 粒子 9: 冲突数 = 1，适应度 = 1/(1+1)=0.51 / (1 + 1) = 0.5
    • 粒子 10: 冲突数 = 2，适应度 = 1/(1+2)=0.3331 / (1 + 2) = 0.333
步骤 3：更新个体最优和全局最优
  • 记录每个粒子的历史最优位置（pbestp_{\text{best}}）和群体的历史最优位置（gbestg_{\text{best}}）。
    • 如果当前适应度优于 pbestp_{\text{best}}，则更新 pbestp_{\text{best}}。
    • 如果当前适应度优于 gbestg_{\text{best}}，则更新 gbestg_{\text{best}}。
步骤 4：更新速度和位置
  • 根据以下公式更新每个粒子的速度和位置：

    vi+1,k=w1vi,k+ϕ1(pbest,k−xi,k)ui+ϕ2(gbest−xi,k)uiv_{i+1,k} = w_1 v_{i,k} + \phi_1 (p_{\text{best},k} - x_{i,k}) u_i + \phi_2 (g_{\text{best}} - x_{i,k}) u_i

    xi+1,k=xi,k+vi+1,kx_{i+1,k} = x_{i,k} + v_{i+1,k}

    其中：
    • vi,kv_{i,k} 是粒子的速度，
    • pbest,kp_{\text{best},k} 是个体的历史最优位置，
    • gbestg_{\text{best}} 是群体的历史最优位置，
    • w1w_1, ϕ1\phi_1, ϕ2\phi_2 是算法的调优参数，
    • uiu_i 是随机数。
步骤 5：迭代
  • 重复上述步骤，直到找到冲突数为 0 的解（即找到八皇后问题的解）。

Python 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
import random

# 定义棋盘大小
N = 8

# 初始化粒子群
def initialize_particles(pop_size):
    """
    随机生成一组粒子（解），每个粒子是一个长度为 N 的列表，表示每行皇后的列位置。
    """
    return [random.sample(range(1, N + 1), N) for _ in range(pop_size)]

# 计算冲突数
def calculate_conflicts(individual):
    """
    计算一个粒子（解）的冲突数。
    """
    conflicts = 0
    for i in range(N):
        for j in range(i + 1, N):
            # 检查列冲突和对角线冲突
            if individual[i] == individual[j] or abs(individual[i] - individual[j]) == abs(i - j):
                conflicts += 1
    return conflicts

# 计算适应度
def fitness(individual):
    """
    计算一个粒子（解）的适应度，冲突数越少，适应度越高。
    """
    return 1 / (1 + calculate_conflicts(individual))

# 更新速度和位置
def update_velocity_position(particle, pbest, gbest, velocity, w, c1, c2):
    """
    更新粒子的速度和位置。
    """
    new_velocity = []
    new_particle = []
    for i in range(N):
        # 更新速度
        v = w * velocity[i] + c1 * random.random() * (pbest[i] - particle[i]) + c2 * random.random() * (gbest[i] - particle[i])
        # 限制速度范围
        v = max(-N, min(N, v))
        new_velocity.append(v)
        # 更新位置
        new_pos = particle[i] + v
        # 限制位置范围
        new_pos = max(1, min(N, new_pos))
        new_particle.append(int(new_pos))
    return new_velocity, new_particle

# 粒子群优化主函数
def pso(pop_size, max_iterations, w, c1, c2):
    """
    粒子群优化主函数。
    """
    particles = initialize_particles(pop_size)
    velocities = [[0] * N for _ in range(pop_size)]
    pbest = particles.copy()
    gbest = min(particles, key=lambda x: calculate_conflicts(x))
    for iteration in range(max_iterations):
        for i in range(pop_size):
            # 更新速度和位置
            velocities[i], particles[i] = update_velocity_position(particles[i], pbest[i], gbest, velocities[i], w, c1, c2)
            # 更新个体最优
            if fitness(particles[i]) > fitness(pbest[i]):
                pbest[i] = particles[i]
            # 更新全局最优
            if fitness(particles[i]) > fitness(gbest):
                gbest = particles[i]
        # 检查是否有解
        if calculate_conflicts(gbest) == 0:
            print(f"Solution found in iteration {iteration}: {gbest}")
            return gbest
    print("No solution found.")
    return None

# 运行粒子群优化
pso(pop_size=10, max_iterations=1000, w=0.5, c1=1.5, c2=1.5)

总结

通过书本中的例子，我们可以看到：

  • 遗传算法（GA） 通过模拟生物进化过程，利用选择、交叉和变异等操作优化种群，适用于复杂的多模态优化问题。
  • 粒子群优化（PSO） 通过模拟群体智能，利用个体和群体的历史最优位置更新搜索方向，适用于连续优化问题。

这两种方法在处理复杂、多模态的优化问题时表现出色，尤其适用于传统梯度方法难以处理的情况。通过具体的示例，我们可以更好地理解 GA 和 PSO 的运行过程。

第六章 约束优化

问题定义

目标函数：

f(x)=(x1−1)2+(x2−2)2f(x) = (x_1 - 1)^2 + (x_2 - 2)^2

约束条件：

  • 等式约束：$ h_1(x) = 2x_1 - x_2 = 0 $
  • 不等式约束：$ g_1(x) = x_1 - 5 \leq 0 $

初始点：

x(0)=(10,−5)x^{(0)} = (10, -5)

迭代 1

我们从初始点 x(0)=(10,−5)x^{(0)} = (10, -5) 开始。

1. 梯度计算
1.1 目标函数的梯度

∇f(x)=[2(x1−1)2(x2−2)]\nabla f(x) = \begin{bmatrix} 2(x_1 - 1) \\ 2(x_2 - 2) \end{bmatrix}

在 x(0)=(10,−5)x^{(0)} = (10, -5) 处：

∇f(x(0))=[2(10−1)2(−5−2)]=[18−14]\nabla f(x^{(0)}) = \begin{bmatrix} 2(10 - 1) \\ 2(-5 - 2) \end{bmatrix} = \begin{bmatrix} 18 \\ -14 \end{bmatrix}

1.2 等式约束的梯度

∇h1(x)=[2−1]\nabla h_1(x) = \begin{bmatrix} 2 \\ -1 \end{bmatrix}

1.3 不等式约束的梯度

∇g1(x)=[10]\nabla g_1(x) = \begin{bmatrix} 1 \\ 0 \end{bmatrix}

2. 二次规划 (QP) 子问题的构造

构造二次规划问题：

min⁡p12pTHkp+∇f(x(k))Tp\min_p \quad \frac{1}{2} p^T H_k p + \nabla f(x^{(k)})^T p

subject to:∇h1(x(k))Tp+h1(x(k))=0∇g1(x(k))Tp+g1(x(k))≤0\text{subject to:} \quad \begin{aligned} & \nabla h_1(x^{(k)})^T p + h_1(x^{(k)}) = 0 \\ & \nabla g_1(x^{(k)})^T p + g_1(x^{(k)}) \leq 0 \end{aligned}

2.1 Hessian 矩阵

目标函数的 Hessian 为：

∇2f(x)=[2002]\nabla^2 f(x) = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix}

由于约束是线性的，其 Hessian 为零。因此：

Hk=[2002]H_k = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix}

2.2 线性化约束

在 x(0)=(10,−5)x^{(0)} = (10, -5) 处：

  • 等式约束：

    h1(x(0))=2(10)−(−5)=25h_1(x^{(0)}) = 2(10) - (-5) = 25

    ∇h1(x(0))Tp+h1(x(0))=2p1−p2+25=0\nabla h_1(x^{(0)})^T p + h_1(x^{(0)}) = 2p_1 - p_2 + 25 = 0

  • 不等式约束：

    g1(x(0))=10−5=5g_1(x^{(0)}) = 10 - 5 = 5

    ∇g1(x(0))Tp+g1(x(0))=p1+5≤0\nabla g_1(x^{(0)})^T p + g_1(x^{(0)}) = p_1 + 5 \leq 0

2.3 子问题的形式

将以上代入，我们的二次规划子问题为：

min⁡p12[p1p2]T[2002][p1p2]+[18−14]T[p1p2]\min_p \quad \frac{1}{2} \begin{bmatrix} p_1 \\ p_2 \end{bmatrix}^T \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix} \begin{bmatrix} p_1 \\ p_2 \end{bmatrix} + \begin{bmatrix} 18 \\ -14 \end{bmatrix}^T \begin{bmatrix} p_1 \\ p_2 \end{bmatrix}

subject to:2p1−p2+25=0p1+5≤0\text{subject to:} \begin{aligned} & 2p_1 - p_2 + 25 = 0 \\ & p_1 + 5 \leq 0 \end{aligned}

2.4 求解子问题

  1. 拉格朗日函数：

    L(p,λ,μ)=p12+p22+18p1−14p2+λ(2p1−p2+25)+μ(p1+5)L(p, \lambda, \mu) = p_1^2 + p_2^2 + 18p_1 - 14p_2 + \lambda (2p_1 - p_2 + 25) + \mu (p_1 + 5)

  2. 对 p1,p2p_1, p_2 求偏导：

    ∂L∂p1=2p1+18+2λ+μ=0⇒p1=−9−λ−μ2\frac{\partial L}{\partial p_1} = 2p_1 + 18 + 2\lambda + \mu = 0 \quad \Rightarrow \quad p_1 = -9 - \lambda - \frac{\mu}{2}

    ∂L∂p2=2p2−14−λ=0⇒p2=7+λ2\frac{\partial L}{\partial p_2} = 2p_2 - 14 - \lambda = 0 \quad \Rightarrow \quad p_2 = 7 + \frac{\lambda}{2}

  3. 带入约束：

    2p1−p2+25=02p_1 - p_2 + 25 = 0

    和

    p1+5≤0p_1 + 5 \leq 0

    解得：

    p1=−5,p2=15p_1 = -5, \quad p_2 = 15

3. 更新点

更新公式：

x(1)=x(0)+αpx^{(1)} = x^{(0)} + \alpha p

假设步长 α=1\alpha = 1：

x(1)=(10,−5)+1⋅(−5,15)=(5,10)x^{(1)} = (10, -5) + 1 \cdot (-5, 15) = (5, 10)

迭代 2

以 x(1)=(5,10)x^{(1)} = (5, 10) 继续迭代。

1. 梯度计算

  1. 目标函数梯度：

    ∇f(x(1))=[2(5−1)2(10−2)]=[816]\nabla f(x^{(1)}) = \begin{bmatrix} 2(5 - 1) \\ 2(10 - 2) \end{bmatrix} = \begin{bmatrix} 8 \\ 16 \end{bmatrix}

  2. 等式约束梯度：

    ∇h1(x(1))=[2−1]\nabla h_1(x^{(1)}) = \begin{bmatrix} 2 \\ -1 \end{bmatrix}

  3. 不等式约束梯度：

    ∇g1(x(1))=[10]\nabla g_1(x^{(1)}) = \begin{bmatrix} 1 \\ 0 \end{bmatrix}

2. 构造二次规划子问题

新的二次规划问题如下：

  1. Hessian:

    Hk=[2002]H_k = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix}

  2. 线性化约束：

    • 等式约束：

      h1(x(1))=2(5)−10=0h_1(x^{(1)}) = 2(5) - 10 = 0

      ∇h1(x(1))Tp+h1(x(1))=2p1−p2=0\nabla h_1(x^{(1)})^T p + h_1(x^{(1)}) = 2p_1 - p_2 = 0

    • 不等式约束：

      g1(x(1))=5−5=0g_1(x^{(1)}) = 5 - 5 = 0

      ∇g1(x(1))Tp+g1(x(1))=p1≤0\nabla g_1(x^{(1)})^T p + g_1(x^{(1)}) = p_1 \leq 0

解得：

p1=−4,p2=−8p_1 = -4, \quad p_2 = -8

3. 更新点

x(2)=x(1)+αp=(5,10)+(−4,−8)=(1,2)x^{(2)} = x^{(1)} + \alpha p = (5, 10) + (-4, -8) = (1, 2)

收敛

在 x(2)=(1,2)x^{(2)} = (1, 2) 处：

  1. 梯度为零：

    ∇f(x(2))=[00]\nabla f(x^{(2)}) = \begin{bmatrix} 0 \\ 0 \end{bmatrix}

  2. 满足约束条件：

    h1(x(2))=2(1)−2=0h_1(x^{(2)}) = 2(1) - 2 = 0

    g1(x(2))=1−5≤0g_1(x^{(2)}) = 1 - 5 \leq 0

因此 x∗=(1,2)x^* = (1, 2) 是最终解。

许可协议

本文采用 署名-非商业性使用-相同方式共享 4.0 国际 许可协议，转载请注明出处。

© 2024 oymaster 使用 Stellar 创建

总访问量次