40、优化算法与MATLAB实现-优快云博客

本文链接：https://blog.youkuaiyun.com/perl8/article/details/155221007

优化算法与MATLAB实现

优化是在满足现有约束条件的同时，最小化或最大化目标函数的计算过程。在当今工程领域，优化技术至关重要，它能帮助解决能源成本高、产品价格和质量的全球竞争以及严格的环境法规等问题。在化学工程行业，优化技术在生产、过程控制、调度、库存控制和运输等方面的应用，显著节省了成本并改善了运营环境。对于化学工程师而言，要在工作中应用优化技术，就必须了解基本算法和计算程序。

本文将聚焦于优化算法以及使用MATLAB实现这些优化技术的程序，通过具体的MATLAB示例来展示这些优化技术的应用，让读者更好地理解其工作原理。

无约束优化算法

无约束优化算法主要用于在没有约束条件的情况下，寻找目标函数的极值。下面将介绍几种常见的无约束优化算法及其MATLAB实现。

Fibonacci方法

Fibonacci方法适用于在最少的试验次数或函数评估下，将包含极值点的区间或不确定区间缩小到给定值。Fibonacci数列的生成公式为：
[
\begin{cases}
F_0 = 0 \
F_1 = 1 \
F_i = F_{i - 1} + F_{i - 2}, i = 2,3,\cdots
\end{cases}
]
区间缩减过程可以表示为：
[
\begin{cases}
I_1 = I_2 + I_3 \
I_2 = I_3 + I_4 \
\cdots \
I_{n - 1} = I_n + I_{n + 1}
\end{cases}
]
基于Fibonacci方法的最小化算法步骤如下：
1. 指定两个点 (x_1) 和 (x_4)，表示初始区间 (I_1 = |x_4 - x_1|)。
2. 指定区间缩减的次数 (n)。如果给定了所需的精度 (\epsilon)，则 (n) 可以定义为满足 (\frac{1}{F_n} < \epsilon) 的最小整数。
3. 计算 (v = \frac{\sqrt{5} - 1}{2})，(w = \frac{1 - \sqrt{5}}{1 + \sqrt{5}})，(\alpha = \frac{v(1 - w^n)}{1 - w^{n + 1}})。
4. 确定中间点 (x_3 = \alpha x_4 + (1 - \alpha)x_1)，并计算 (f_3 = f(x_3))。
5. 进行 (n - 1) 次迭代：
- 如果 (i = n - 1)，则 (x_2 = 0.01x_1 + 0.99x_3)；否则 (x_2 = \alpha x_1 + (1 - \alpha)x_4)。
- 计算 (f_2 = f(x_2))。
- 如果 (f_2 < f_3)，则 (x_4 = x_3)，(x_3 = x_2)，(f_3 = f_2)；否则 (x_1 = x_4)，(x_4 = x_2)。
- 更新 (\alpha = \frac{v(1 - w^{n - i})}{1 - w^{n - i + 1}})。

以下是实现Fibonacci搜索算法的MATLAB函数 fbnopt ：

function [x, f, fint] = fbnopt(objfun, a, b, n)
% fbnopt.m: 1-dimensional Fibonacci search method
% inputs:
%    objfun: objective function
%    a, b: initial interval points
%    n: number of reduction
% output:
%    fint: size of final interval
%    x: optimum point
%    f: function value at the optimum point
v = (sqrt(5)-1)/2; w = (1-sqrt(5))/(1+sqrt(5)); x1 = a; x4 = b; alpha = 
v*(1-w^n)/(1-w^(n+1));
x3 = alpha*x4 + (1-alpha)*x1; f3 = objfun(x3);
for k = 1: n - 1
     if k == n - 1
           x2 = 0.01*x1 + 0.99*x3;
     else
           x2 = alpha*x1 + (1-alpha)*x4;
     end
     f2 = objfun(x2);
     if (f2 < f3)
            x4 = x3; x3 = x2; f3 = f2;
     else
            x1 = x4; x4 = x2; f4 = f2;
     end
     alpha = v*(1- w^(n - k)) / (1-w^(n - k + 1));
end
x = x3; f = f3; fint = abs(x1 - x4);
end

示例10.1：Fibonacci搜索
使用Fibonacci搜索方法，寻找函数 (f(x) = 2x^2\sin(1.5x) - 4x^2 + 3x - 1) 在区间 ([-8, 8]) 内的最小值，区间缩减次数 (n = 20)。

a = -8; b = 8; n = 20; fobj = @(x) 2*x^2*sin(1.5*x) - 4*x^2 + 3*x-1;
[x, f, fint] = fbnopt(fobj, a, b, n)

运行结果：

x =
   -5.7256
f =
   -197.9705
fint =
      0.0015

黄金分割法

从区间缩减过程可知，当迭代次数增加时，Fibonacci数的比值 (\frac{F_{n - 1}}{F_n}) 趋近于黄金分割比 (\tau = \frac{\sqrt{5} - 1}{2} \approx 0.61803)。

黄金分割法的区间缩减方案遵循Fibonacci算法的步骤，其算法步骤如下：
1. 指定初始点 (x_1) 和步长 (h)，计算 (f_1 = f(x_1))。
2. 令 (x_2 = x_1 + h)，计算 (f_2 = f(x_2))。
3. 如果 (f_2 > f_1)，交换点1和点2，并将 (h) 取反。
4. 令 (h = \tau h)，(x_4 = x_2 + h)，计算 (f_4 = f(x_4))。
5. 如果 (f_4 > f_2)，转到步骤7。
6. 令 (x_1 = x_2)，(x_2 = x_4)，转到步骤4。
7. 令 (x_3 = \tau x_4 + (1 - \tau)x_1)，计算 (f_3 = f(x_3))。
8. 如果 (f_2 < f_3)，则 (x_4 = x_1)，(x_1 = x_3)；否则 (x_1 = x_2)，(x_2 = x_3)，(f_2 = f_3)。
9. 如果不满足停止准则，转到步骤7。

以下是实现黄金分割搜索算法的MATLAB函数 gsopt ：

function [x, f, n] = gsopt(objfun, x1, h, crit)
% gsopt.m: golden section search method
% input:
%    objfun: objective function
%    x1: initial point
%    h: step size
%    crit: stopping criterion
% output:
%    x: optimal point
%    f: function value at the optimal point
%    n: number of function evaluations
tau = (sqrt(5) - 1)/2; n = 0; f1 = objfun(x1); n = n+1;
x2 = x1 + h; f2 = objfun(x2); n = n+1; 
if f2 > f1
      temp = x1; x1 = x2; x2 = temp; temp = f1; f1 = f2; f2 = temp; h = -h;
end
while x1 < 1e50
      h = h/tau; x4 = x2 + h; f4 = objfun(x4); n = n+1;
      if f4 > f2, break; end
      f1 = f2; x1 = x2; f2 = f4; x2 = x4;
end
fold = (f1 + f2 + f4) / 3; ind = 0;
while x1 < 1e50
      if abs(x4-x1) < crit, break; end
      x3 = tau*x4 + (1-tau)*x1; f3 = objfun(x3); n = n+1;
      if f2 < f3
            x4 = x1; x1 = x3; f4 = f1; f1 = f3;
      else
            x1 = x2; x2 = x3; f1 = f2; f2 = f3;
      end
      fpr = (f1 + f2 + f4) / 3;
      if abs(fpr - fold) < crit
           ind = ind + 1;
           if ind == 2, break; end
      else
           ind = 0;
      end
      fold = fpr;
end
x = x2; f = f2;
end

示例10.2：黄金分割法
使用黄金分割搜索方法，寻找函数 (f(x) = 2x^2\sin(1.5x) - 4x^2 + 3x - 1) 的最小值。初始点 (x_1 = -5)，步长 (h = 0.1)，停止准则为 (1\times10^{-6})。

crit = 1e-6; x1 = -5; h = 0.1;
fobj = @(x) 2*x^2*sin(1.5*x) - 4*x^2 + 3*x-1;
[x, f, n] = gsopt(fobj, x1, h, crit)

运行结果：

x =
    -5.7248
f =
  -197.9705
n =
     25

需要注意的是，该目标函数有多个局部最小值，不同的初始点会得到不同的结果。例如，当 (x_1 = 1) 时，(x_{min} = 3.7939)，(f_{min} = -63.265)；当 (x_1 = 6) 时，(x_{min} = 7.6665)，(f_{min} = -316.0254)。

x1 = 1; [x, f, n] = gsopt(fobj, x1, h, crit)
x1 = 6; [x, f, n] = gsopt(fobj, x1, h, crit)

Brent二次拟合方法

假设给定三个初始数据点 ((x_1, f_1))，((x_2, f_2))，((x_3, f_3))，其中 (x_1 < x_2 < x_3) 且 (f_2 \leq \min(f_1, f_3))，可以通过这些点拟合一个二次函数 (p(x) = a + bx + cx^2)。使用拉格朗日多项式，函数 (p(x)) 可以表示为：
[
p(x) = \frac{f_1(x - x_2)(x - x_3)}{(x_1 - x_2)(x_1 - x_3)} + \frac{f_2(x - x_1)(x - x_3)}{(x_2 - x_1)(x_2 - x_3)} + \frac{f_3(x - x_1)(x - x_2)}{(x_3 - x_1)(x_3 - x_2)}
]
通过令 (\frac{dp(x)}{dx} = 0)，可以找到最小值点 (x^* = x_4)。

Brent二次拟合算法步骤如下：
1. 指定三个点 (a)，(b)（形成区间）和 (x)（给出最小函数值）。
2. 将 (w) 和 (v) 初始化为 (x)。
3. 如果 (x)，(w) 和 (v) 都不同，则转到步骤5。
4. 使用黄金分割搜索方法，在 (x - a) 或 (x - b) 中较大的区间内找到点 (u)，转到步骤7。
5. 对 (x)，(w) 和 (v) 进行二次拟合，找到最小值点 (u)。
6. 将 (u) 调整到 (x - a) 或 (x - b) 中较大的区间内。
7. 计算 (f_u = f(u))。
8. 找到 (a)，(b)，(x)，(w) 和 (v) 的新值。
9. 检查 (x - a) 或 (x - b) 中较大的区间是否小于指定的收敛准则。如果未达到收敛，转到步骤3。

以下是实现Brent搜索算法的MATLAB函数 brentopt ：

function [xopt, fopt, nf] = brentopt(objfun, x1, h, crit)
% brentopt.m: minimization(1-dimensional search) by Brent's algorithm
% input:
%    objfun: objective function
%    x1: initial point
%    h: step size
%    crit: convergence criterion
% output:
%    x: optimum point
%    f: function value at the optimum point
%    hf: number of function evaluations for search
clarge = 1e40; tau = (sqrt(5) - 1)/2; nf = 1; f1 = objfun(x1);
x2 = x1 + h; nf = nf+1; f2 = objfun(x2);
if f2 > f1
     temp = x1; x1 = x2; x2 = temp; temp = f1; f1 = f2; f2 = temp; h = -h;
end
while x1 < clarge
     h = h/tau; x4 = x2 + h; nf = nf+1; f4 = objfun(x4);
     if f4 > f2, break; end
     f1 = f2; x1 = x2; f2 = f4; x2 = x4;
end
x3 = x4; f3 = f4; a = x1; b = x3; fa = f1; fb = f3; ha = a; hb = b;
if (b < a), ha = b; hb = a; end
x = x2; fx = fb; w = x; v = w; ev = 0; fx = f2; fw = fx; fv = fw;
while 1
    hm = (ha + hb)/2;
    if (abs(x - hm) <= crit - (hb - ha)/2), x2 = x; f2 = fx; return; end
    etemp = 0; p = 0; q = 0;
    if (abs(ev) > crit)
         r = (x - w)*(fx - fv); q = (x - v)*(fx - fw);
         p = (x - v)*q - (x - w)*r; q = 2*(q - r);
         if (q > 0), p = -p;
         else, q = -q; end
         etemp = ev; ev = d; 
    end
    ind1 = 1;
    if ((q * (ha - x) - p) < 0), ind1 = -1; end
    ind2 = 1;
    if ((q * (hb - x) - p) < 0), ind2 = -1; end
    if ((abs(p) >= abs(q*etemp/2)) | (ind1 == ind2))
          if (x < hm), ev = hb - x;
          else, ev = ha - x; end
          d = (1 - tau)*ev;
    else
        d = p/q; u = x + d;
        if (((u - ha) < crit) | ((hb - u) < crit))
            if (x < hm), d = crit;
            else, d = -crit; end
        end
    end
    if (abs(d) >= crit)
        u = x + d;
    else
        if (d > 0), u = x + crit; else, u = x - crit; end
    end
    nf = nf+1; fu = objfun(u);
    if (fu <= fx)
        if (u < x), hb = x; fb = fx;
        else, ha = x; fa = fx; end
        v = w; fv = fw; w = x; fw = fx; x = u; fx = fu;
    else
        if (u < x), ha = u; fa = fu;
        else, hb = u; fb = fu; end
        if ((fu <= fw) | (w == x))
            v = w; fv = fw; w = u; fw = fu;
        elseif ((fu <= fv) | (v == x) | (v == w))
            v = u; fv = fu;
        end
    end
    if ((fa - fx) + (fb - fx) < crit), x2 = x; f2 = fx; return; end
    xopt = x2; fopt = f2;
end

示例10.3：Brent算法
使用Brent搜索方法，寻找函数 (f(x) = \exp(x) - 3x + \frac{0.02}{x} - \frac{0.00004}{x}) 的最小值。初始点 (x_1 = 0.01)，步长 (h = 0.2)，停止准则为 (1\times10^{-6})。

x1 = 0.01; h = 0.2; crit = 1e-6;
fun = @(x) exp(x) -3*x + 0.02/x - 0.00004/x;
[xopt, fopt, nf] = brentopt(fun, x1, h, crit)

运行结果：

xopt =
     1.0572
fopt =
    -0.2744
nf =
     12

Shubert - Piyavskii方法

考虑一个简单的最大化问题：
[
\begin{cases}
\max f(x) \
a \leq x \leq b
\end{cases}
]
假设函数 (f(x)) 是Lipschitz连续的，即存在一个已知的常数 (C)，使得对于任意的 (x, y \in [a, b])，有 (|f(x) - f(y)| \leq C|x - y|)。

Shubert - Piyavskii算法步骤如下：
1. 选择区间 ([a, b]) 的中点作为第一个采样点 (x_0 = \frac{a + b}{2})。
2. 在 (x_0) 处画一对斜率为 (C) 的直线 (y(x) = y_0 + C|x - x_0|)，这对直线与区间 ([a, b]) 的端点的交点给出一个由两个点组成的初始锯齿形。
3. 构造锯齿形顶点序列 ([(t_1, z_1), (t_2, z_2), \cdots, (t_n, z_n)])（(z_1 \leq z_2 \leq \cdots \leq z_n)）。
4. 确定新的锯齿形覆盖 ([(t_1, z_1), (t_2, z_2), \cdots, (t_{n - 1}, z_{n - 1}), (t_l, z_l), (t_r, z_r)])。
5. 当 (|z_n - y_n| < \epsilon) 时，停止该过程。

以下是实现Shubert - Piyavskii算法的MATLAB函数 shubertopt ：

function [xopt, fopt, nf] = shubertopt(objfun, a, b, C, crit, nfmax)
% shubertopt.m: optimization by Shubert-Piyavskii algorithm
%                   1-D maximization of Lipschitz functions
% inputs:
%    objfun: objective function
%    a,b: end points of initial interval
%    C: Lipschitz constant
%    crit: tolerance
%    nfmax: maximum number of function evaluations
% outputs:
%    xopt: optimum point
%    fopt: function value at the optimum point
%    nf: number of function evaluations
nf = 0; x0 = (a + b)/2; nf = nf + 1; y0 = objfun(x0); ymax = y0; xmax = 
x0; fmax = y0 + C*(b - a)/2;
T(1) = b; Z(1) = y0 + C*(b - a)/2; T(2) = a; Z(2) = y0 + C*(b - a)/2; n = 2;
while ((fmax - ymax) > crit & nf <= nfmax)
    tn = T(n); zn = Z(n); nf = nf + 1; yn = objfun(tn);
    if (yn > ymax), ymax = yn; xmax = tn; end
    zL = (zn + yn)/2; zR = zL;
    tL = tn - (zn - yn)/2/C; tR = tn + (zn - yn)/2/C;
    ind1 = 0; ind2 = 0;
    if (tL >= a & tL <= b), ind1 = 1; end
    if (tR >= a & tR <= b), ind2 = 1; end
    if (ind1 == 1 & ind2 == 0)
          T(n) = tL; Z(n) = zL;
    elseif (ind1 == 0 & ind2 == 1)
          T(n) = tR; Z(n) = zR;
    elseif (ind1 == 1 & ind2 == 1)
          T(n) = tL; Z(n) = zL; T(n+1) = tR; Z(n+1) = zR; n = n+1;
    end
    [Z, Indx] = sort(Z); T = T(Indx); fmax = Z(n);
end
xopt = xmax; fopt = ymax;
end

总结

本文介绍了几种常见的无约束优化算法，包括Fibonacci方法、黄金分割法、Brent二次拟合方法和Shubert - Piyavskii方法，并给出了相应的MATLAB实现。这些算法在不同的场景下具有各自的优势，读者可以根据具体问题选择合适的算法。同时，通过MATLAB的实现，读者可以方便地将这些算法应用到实际问题中，提高工作效率。

算法名称	适用场景	特点
Fibonacci方法	需在最少试验次数或函数评估下缩小区间	利用Fibonacci数列进行区间缩减
黄金分割法	通用的无约束优化	基于黄金分割比进行区间缩减
Brent二次拟合方法	有一定数据点可用于二次拟合的情况	通过二次拟合寻找最小值点
Shubert - Piyavskii方法	目标函数Lipschitz连续的最大化问题	构造锯齿形覆盖寻找全局最大值

mermaid流程图示例（以Fibonacci方法为例）：

graph TD;
    A[指定初始区间x1,x4] --> B[指定区间缩减次数n];
    B --> C[计算v,w,α];
    C --> D[确定中间点x3并计算f3];
    D --> E{i = 1: n - 1};
    E -- 是 --> F{是否i = n - 1};
    F -- 是 --> G[x2 = 0.01x1 + 0.99x3];
    F -- 否 --> H[x2 = αx1 + (1 - α)x4];
    G --> I[计算f2];
    H --> I;
    I --> J{f2 < f3};
    J -- 是 --> K[x4 = x3, x3 = x2, f3 = f2];
    J -- 否 --> L[x1 = x4, x4 = x2];
    K --> M[更新α];
    L --> M;
    M --> E;
    E -- 否 --> N[输出结果x,f,fint];

优化算法与MATLAB实现

各算法复杂度和性能对比

为了更清晰地了解上述几种无约束优化算法的特点，我们对它们的复杂度和性能进行简单对比。

算法名称	时间复杂度	空间复杂度	性能特点
Fibonacci方法	与区间缩减次数 (n) 相关，整体复杂度适中	主要存储几个关键变量，空间复杂度低	适合对区间缩减次数有要求的场景，能在较少的试验次数内达到一定精度
黄金分割法	收敛速度相对稳定，复杂度与迭代次数有关	存储少量变量，空间开销小	通用性强，在大多数无约束优化问题中都能较好地工作，但可能需要较多的迭代次数
Brent二次拟合方法	复杂度受二次拟合和区间调整影响	存储多个中间变量，空间复杂度稍高	利用二次拟合能更快地逼近最小值点，在有一定数据点可供拟合时效率较高
Shubert - Piyavskii方法	复杂度与锯齿形覆盖的构造和更新次数有关	存储锯齿形顶点序列，空间开销较大	专门用于目标函数Lipschitz连续的最大化问题，能有效寻找全局最大值

从性能上看，Fibonacci方法和黄金分割法更侧重于区间的缩减，而Brent二次拟合方法则通过二次函数拟合来加速收敛。Shubert - Piyavskii方法则针对特定类型的最大化问题，利用锯齿形覆盖逐步逼近全局最优解。

实际应用中的选择建议

在实际应用中，选择合适的优化算法至关重要。以下是一些选择建议：
1. 精度要求和试验次数 ：如果对精度要求较高，且希望在较少的试验次数内达到目标，Fibonacci方法是一个不错的选择。例如，在一些资源有限的实验场景中，需要在有限的测试次数内找到最优解。
2. 通用性和简单性 ：黄金分割法具有较强的通用性，实现简单，适用于大多数无约束优化问题。当对问题的特性了解较少时，可以优先考虑使用黄金分割法。
3. 数据拟合和快速收敛 ：如果有一定的初始数据点，并且希望能够快速收敛到最小值点，Brent二次拟合方法更为合适。在一些已知部分数据特征的工程问题中，该方法可以发挥出更好的性能。
4. 特定类型的最大化问题 ：对于目标函数Lipschitz连续的最大化问题，Shubert - Piyavskii方法是专门为这类问题设计的，能够有效地找到全局最大值。

优化算法的拓展和改进

虽然上述几种算法在无约束优化问题中表现出色，但在实际应用中，可能会遇到更复杂的情况。以下是一些对优化算法的拓展和改进思路：
1. 结合多种算法 ：可以将不同的优化算法结合起来使用，充分发挥它们的优势。例如，先使用黄金分割法进行初步的区间搜索，然后再使用Brent二次拟合方法进行精确的局部优化。
2. 自适应调整参数 ：在算法运行过程中，根据问题的特点和当前的搜索状态，自适应地调整算法的参数，以提高算法的性能。例如，在Fibonacci方法中，根据区间的变化情况动态调整区间缩减次数。
3. 并行计算 ：对于大规模的优化问题，可以利用并行计算技术，同时进行多个搜索过程，加快搜索速度。例如，将不同的初始点分配给不同的处理器，并行地进行优化计算。

总结与展望

本文详细介绍了几种常见的无约束优化算法，包括Fibonacci方法、黄金分割法、Brent二次拟合方法和Shubert - Piyavskii方法，并给出了相应的MATLAB实现。这些算法在不同的场景下具有各自的优势，读者可以根据具体问题选择合适的算法。

在未来的优化问题中，随着问题规模的不断增大和复杂度的不断提高，对优化算法的性能和效率提出了更高的要求。一方面，需要进一步研究和改进现有的优化算法，提高其收敛速度和精度；另一方面，需要探索新的优化算法和方法，以应对更加复杂的实际问题。同时，结合人工智能和机器学习技术，将优化算法应用到更广泛的领域，如数据分析、机器学习模型训练等，也是未来的发展方向之一。

mermaid流程图示例（以结合黄金分割法和Brent二次拟合方法为例）：

graph TD;
    A[开始] --> B[使用黄金分割法进行初步搜索];
    B --> C{是否满足初步精度要求};
    C -- 是 --> D[使用Brent二次拟合方法进行精确优化];
    C -- 否 --> B;
    D --> E{是否满足最终精度要求};
    E -- 是 --> F[输出结果];
    E -- 否 --> D;

通过以上的介绍和分析，希望读者能够对无约束优化算法有更深入的理解，并能够在实际应用中灵活运用这些算法，解决各种优化问题。同时，也鼓励读者进一步探索和研究优化算法的相关理论和技术，为优化领域的发展做出贡献。