【MATLAB学习笔记】数值方法——黄金分割法（求极小值）

前言

  黄金分割法是一种优化算法，用于求解函数的极小值。它基于黄金比例，通过不断缩小搜索范围来逼近极小值点。该方法可以在一定程度上避免了传统的梯度下降法和牛顿法可能陷入局部最优解的问题。

具体步骤如下：

1. 确定初始搜索范围$[a,b]$，使得函数在这个范围内有极小值点。
2. 计算两个内点${\alpha }_1$和${\alpha }_2$，使得${\alpha }_1$和${\alpha }_2$分别与$a$和$b$的距离满足黄金比例。
3. 计算函数在${\alpha }_1$和${\alpha }_2$的取值，选择较小的值所对应的点作为新的搜索范围的端点。
4. 不断重复步骤2和3，直到搜索范围足够小，可以认为找到了极小值点。

通过黄金分割法，可以逐步逼近函数的极小值点，是一种有效的优化算法。

算法流程图

算法代码

function Hjfg_minValue(f_sym,x0,eps,N)
    % f_sym为所求函数, 符号函数
    % x0为初始点，默认为0
    % eps为允许误差, 默认为1e-4
    % N为最大迭代次数, 默认为500
    if nargin < 2, x0 = 0; end
    if nargin < 3, eps = 1e-4; end
    if nargin < 4, N = 500; end
    
    f = matlabFunction(f_sym);
    lamda = (sqrt(5)-1)/2;
    a = x0(1);
    b = x0(2);
    x1 = b - lamda*(b-a);
    x2 = a + lamda*(b-a);
    y1 = f(x1);
    y2 = f(x2);

    k = 0;
    fprintf('函数f=%s迭代开始:\n',f_sym)
    while true
        k = k + 1;

        fprintf('k=%2.d, a=%.4f, x1=%.4f, x2=%.4f, b=%.4f\n',k,a,x1,x2,b);
        
        if y1 >= y2
            a = x1;
            x1 = x2;
            y1 = y2;

            x2 = a + lamda*(b-a);
            y2 = f(x2);
        else
            b = x2;
            x2 = x1;
            y2 = y1;

            x1 = b - lamda*(b-a);
            y1 = f(x1);
        end


        if abs((b-a)/(b+1e-2)) < eps && abs((y2-y1)/(y2+1e-2)) < eps
            x_star = (1/2)*(a+b);
            fprintf('已收敛!结果如下：\nk=%2.d, x1=%.4f, x2=%.4f, x_star=%.4f, f_star=%.4f',...
                k+1,x1,x2,x_star,f(x_star));
            break
        end 


        if k > N
            fprintf('已经超出最大迭代次数%d！\n',N)
            break
        end
        
    end
    fprintf('\n\n')

end

测试

测试函数如下：

syms x
f1_sym = 2*x.^2 - x - 1;
f2_sym = x.^2;
f3_sym = 3*x.^2 + 1;

解析解：

f1 = diff(f1_sym);
x1_min = double(vpa(solve(f1,x)));
f2 = diff(f2_sym);
x2_min = double(vpa(solve(f2,x)));
f3 = diff(f3_sym);
x3_min = double(vpa(solve(f3,x)));

x1_min	x2_min	x3_min
0.25	0	0

算法解：

x01 = [-1,1];
x02 = [0,1];
x03 = [-1,1];

eps1 = 0.16;
eps2 = 0.2;
eps3 = 0.04;

N = 100;
Hjfg_minValue(f1_sym,x01,eps1,N)
Hjfg_minValue(f2_sym,x02,eps2,N)
Hjfg_minValue(f3_sym,x03,eps3,N)

迭代过程如下：

函数f=2*x^2 - x - 1迭代开始:
k= 1, a=-1.0000, x1=-0.2361, x2=0.2361, b=1.0000
k= 2, a=-0.2361, x1=0.2361, x2=0.5279, b=1.0000
k= 3, a=-0.2361, x1=0.0557, x2=0.2361, b=0.5279
k= 4, a=0.0557, x1=0.2361, x2=0.3475, b=0.5279
k= 5, a=0.0557, x1=0.1672, x2=0.2361, b=0.3475
k= 6, a=0.1672, x1=0.2361, x2=0.2786, b=0.3475
k= 7, a=0.1672, x1=0.2098, x2=0.2361, b=0.2786
k= 8, a=0.2098, x1=0.2361, x2=0.2523, b=0.2786
已收敛!结果如下：
k= 9, x1=0.2523, x2=0.2624, x_star=0.2574, f_star=-1.1249

函数f=x^2迭代开始:
k= 1, a=0.0000, x1=0.3820, x2=0.6180, b=1.0000
k= 2, a=0.0000, x1=0.2361, x2=0.3820, b=0.6180
k= 3, a=0.0000, x1=0.1459, x2=0.2361, b=0.3820
k= 4, a=0.0000, x1=0.0902, x2=0.1459, b=0.2361
k= 5, a=0.0000, x1=0.0557, x2=0.0902, b=0.1459
k= 6, a=0.0000, x1=0.0344, x2=0.0557, b=0.0902
k= 7, a=0.0000, x1=0.0213, x2=0.0344, b=0.0557
k= 8, a=0.0000, x1=0.0132, x2=0.0213, b=0.0344
k= 9, a=0.0000, x1=0.0081, x2=0.0132, b=0.0213
k=10, a=0.0000, x1=0.0050, x2=0.0081, b=0.0132
k=11, a=0.0000, x1=0.0031, x2=0.0050, b=0.0081
k=12, a=0.0000, x1=0.0019, x2=0.0031, b=0.0050
k=13, a=0.0000, x1=0.0012, x2=0.0019, b=0.0031
已收敛!结果如下：
k=14, x1=0.0007, x2=0.0012, x_star=0.0010, f_star=0.0000

函数f=3*x^2 + 1迭代开始:
k= 1, a=-1.0000, x1=-0.2361, x2=0.2361, b=1.0000
k= 2, a=-0.2361, x1=0.2361, x2=0.5279, b=1.0000
k= 3, a=-0.2361, x1=0.0557, x2=0.2361, b=0.5279
k= 4, a=-0.2361, x1=-0.0557, x2=0.0557, b=0.2361
k= 5, a=-0.0557, x1=0.0557, x2=0.1246, b=0.2361
k= 6, a=-0.0557, x1=0.0132, x2=0.0557, b=0.1246
k= 7, a=-0.0557, x1=-0.0132, x2=0.0132, b=0.0557
k= 8, a=-0.0132, x1=0.0132, x2=0.0294, b=0.0557
k= 9, a=-0.0132, x1=0.0031, x2=0.0132, b=0.0294
k=10, a=-0.0132, x1=-0.0031, x2=0.0031, b=0.0132
k=11, a=-0.0031, x1=0.0031, x2=0.0069, b=0.0132
k=12, a=-0.0031, x1=0.0007, x2=0.0031, b=0.0069
k=13, a=-0.0031, x1=-0.0007, x2=0.0007, b=0.0031
k=14, a=-0.0007, x1=0.0007, x2=0.0016, b=0.0031
k=15, a=-0.0007, x1=0.0002, x2=0.0007, b=0.0016
k=16, a=-0.0007, x1=-0.0002, x2=0.0002, b=0.0007
k=17, a=-0.0002, x1=0.0002, x2=0.0004, b=0.0007
k=18, a=-0.0002, x1=0.0000, x2=0.0002, b=0.0004
已收敛!结果如下：
k=19, x1=-0.0000, x2=0.0000, x_star=-0.0000, f_star=1.0000

总代码

clc;clear;close all
%% 
syms x
f1_sym = 2*x.^2 - x - 1;
f2_sym = x.^2;
f3_sym = 3*x.^2 + 1;
% 解析解
f1 = diff(f1_sym);
x1_min = double(vpa(solve(f1,x)));
f2 = diff(f2_sym);
x2_min = double(vpa(solve(f2,x)));
f3 = diff(f3_sym);
x3_min = double(vpa(solve(f3,x)));

x01 = [-1,1];
x02 = [0,1];
x03 = [-1,1];

eps1 = 0.16;
eps2 = 0.2;
eps3 = 0.04;

N = 100;
Hjfg_minValue(f1_sym,x01,eps1,N)
Hjfg_minValue(f2_sym,x02,eps2,N)
Hjfg_minValue(f3_sym,x03,eps3,N)

function Hjfg_minValue(f_sym,x0,eps,N)
    % f_sym为所求函数, 符号函数
    % x0为初始点，默认为0
    % eps为允许误差, 默认为1e-4
    % N为最大迭代次数, 默认为500
    if nargin < 2, x0 = 0; end
    if nargin < 3, eps = 1e-4; end
    if nargin < 4, N = 500; end
    
    f = matlabFunction(f_sym);
    lamda = (sqrt(5)-1)/2;
    a = x0(1);
    b = x0(2);
    x1 = b - lamda*(b-a);
    x2 = a + lamda*(b-a);
    y1 = f(x1);
    y2 = f(x2);

    k = 0;
    fprintf('函数f=%s迭代开始:\n',f_sym)
    while true
        k = k + 1;

        fprintf('k=%2.d, a=%.4f, x1=%.4f, x2=%.4f, b=%.4f\n',k,a,x1,x2,b);
        
        if y1 >= y2
            a = x1;
            x1 = x2;
            y1 = y2;

            x2 = a + lamda*(b-a);
            y2 = f(x2);
        else
            b = x2;
            x2 = x1;
            y2 = y1;

            x1 = b - lamda*(b-a);
            y1 = f(x1);
        end


        if abs((b-a)/(b+1e-2)) < eps && abs((y2-y1)/(y2+1e-2)) < eps
            x_star = (1/2)*(a+b);
            fprintf('已收敛!结果如下：\nk=%2.d, x1=%.4f, x2=%.4f, x_star=%.4f, f_star=%.4f',...
                k+1,x1,x2,x_star,f(x_star));
            break
        end 


        if k > N
            fprintf('已经超出最大迭代次数%d！\n',N)
            break
        end
        
    end
    fprintf('\n\n')

end

总结

  这只是一个基础的示例，实际中还会有更具体的、更细致的要求，这就需要再做额外调整；另外本人也仍在学习中，这只是个人的学习笔记，可能还有一些不足之处，欢迎指正。

参考文献

主编白清顺. 机械优化设计.第6版 [M]. 2017.