最优化方法——Matlab实现黄金分割法一维搜索

黄金分割法一维搜索原理

若保留区间为[x1,b],我们得到的结果是一致的.
该方法称为黄金分割法,实际计算取近似值: x1=a+0.382(b – a), x2=a+0.618(b – a),
所以黄金分割法又称为0.618法.
黄金分割法每次缩小区间的比例是一致的，每次将区间长度缩小到原来的0.618倍.

算法流程：

黄金分割法也称作0.618法，一维指的是只含有一个未知量的情况。

Matlab代码

用matlab实现黄金分割法求解f(x)=x^2-x+2在（-1,3）上的最小值：

clc,clear,close all;
a = -1; b =3;
ep = 0.08*(b-a);
x = a:0.1:b;
f_x = x.^2-x+2;
plot(x, f_x, 'linewidth', 1.5)
axis([-1, 3, 0, 8])
title('f(x)=x^2-x+2')
grid on;
flag = 0;
cnt = 0;
pause(0.5)
while 1
    fprintf('第%d次迭代：\n', cnt)
    if flag==0
        x2 = a + 0.618*(b-a);
        f2 = x2.^2-x2+2;
        x1 = a + b - x2;
        f1 = x1.^2-x1+2;
        fprintf('a = %f, b = %f\n', a, b)
        fprintf('x1 = %f, x2 = %f, f1 = %f, f2 = %f\n', x1, x2, f1, f2)
        hold on
        stem([x2, x1], [f2, f1], '--', 'linewidth', 0.8)
        pause(1)
    elseif flag==1
        x1 = a + b - x2;
        f1 = x1.^2-x1+2;
        fprintf('a = %f, b = %f\n', a, b)
        fprintf('x1 = %f, x2 = %f, f1 = %f, f2 = %f\n', x1, x2, f1, f2)
        stem([x2, x1], [f2, f1], '--', 'linewidth', 0.8)
        pause(1)
    elseif flag==2
        x2 = a + 0.618*(b-a);
        f2 = x2.^2-x2+2;
        fprintf('a = %f, b = %f\n', a, b)
        fprintf('x1 = %f, x2 = %f, f1 = %f, f2 = %f\n', x1, x2, f1, f2)
        stem([x2, x1], [f2, f1], '--', 'linewidth', 0.8)
        pause(1)
    end
    if abs(b-a)<ep
        xb = (a+b)/2;
        disp('最优解为：')
        fprintf('xb = %f, f(xb) = %f\n', xb, xb.^2-xb+2)
        disp('黄金分割法一维搜索完毕.')
        break
    elseif f1<f2
        disp('f1)
        b = x2;
        x2 = x1;
        f2 = f1;
        flag = 1;
    elseif f1==f2
        disp('f1=f2')
        a = x1;
        b = x2;
        flag = 0;
    elseif f1>f2
        disp('f1>f2')
        a = x1;
        x1 = x2;
        f1 = f2;
        flag = 2;
    end
    cnt = cnt + 1;
end
pause(0.5)
stem(xb, xb^2-xb+2, 'r', 'linewidth', 2)
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

代码运行有动态效果，这里就不再保存为GIF动图了，可以复制一键运行尝试：

命令行窗口结果打印：

第0次迭代：
a = -1.000000, b = 3.000000
x1 = 0.528000, x2 = 1.472000, f1 = 1.750784, f2 = 2.694784
f1<f2
第1次迭代：
a = -1.000000, b = 1.472000
x1 = -0.056000, x2 = 0.528000, f1 = 2.059136, f2 = 1.750784
f1>f2
第2次迭代：
a = -0.056000, b = 1.472000
x1 = 0.528000, x2 = 0.888304, f1 = 1.750784, f2 = 1.900780
f1<f2
第3次迭代：
a = -0.056000, b = 0.888304
x1 = 0.304304, x2 = 0.528000, f1 = 1.788297, f2 = 1.750784
f1>f2
第4次迭代：
a = 0.304304, b = 0.888304
x1 = 0.528000, x2 = 0.665216, f1 = 1.750784, f2 = 1.777296
f1<f2
第5次迭代：
a = 0.304304, b = 0.665216
x1 = 0.441520, x2 = 0.528000, f1 = 1.753420, f2 = 1.750784
f1>f2
第6次迭代：
a = 0.441520, b = 0.665216
x1 = 0.528000, x2 = 0.579764, f1 = 1.750784, f2 = 1.756362
最优解为：
xb = 0.553368, f(xb) = 1.752848
黄金分割法一维搜索完毕.
>> 
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

更换匿名函数：

通过更改目标函数 f_x ，对自定义的目标函数进行一维搜索的代码：

clc,clear,close all;
a = -1; b =3;
ep = 0.08*(b-a);
x = a:0.1:b;
f_x = @(x)x.^2-3*x+2;
plot(x, f_x(x), 'linewidth', 1.5)
axis tight
tl = func2str(f_x);
title(tl(5:end))
grid on;
flag = 0;
cnt = 0;
pause(0.5)
while 1
    fprintf('第%d次迭代：\n', cnt)
    if flag==0
        x2 = a + 0.618*(b-a);
        f2 = f_x(x2);
        x1 = a + b - x2;
        f1 = f_x(x1);
        fprintf('a = %f, b = %f\n', a, b)
        fprintf('x1 = %f, x2 = %f, f1 = %f, f2 = %f\n', x1, x2, f1, f2)
        hold on
        stem([x2, x1], [f2, f1], '--', 'linewidth', 0.8)
        pause(1)
    elseif flag==1
        x1 = a + b - x2;
        f1 = f_x(x1);
        fprintf('a = %f, b = %f\n', a, b)
        fprintf('x1 = %f, x2 = %f, f1 = %f, f2 = %f\n', x1, x2, f1, f2)
        stem([x2, x1], [f2, f1], '--', 'linewidth', 0.8)
        pause(1)
    elseif flag==2
        x2 = a + 0.618*(b-a);
        f2 = f_x(x2);
        fprintf('a = %f, b = %f\n', a, b)
        fprintf('x1 = %f, x2 = %f, f1 = %f, f2 = %f\n', x1, x2, f1, f2)
        stem([x2, x1], [f2, f1], '--', 'linewidth', 0.8)
        pause(1)
    end
    if abs(b-a)<ep
        xb = (a+b)/2;
        disp('最优解为：')
        fprintf('xb = %f, f(xb) = %f\n', xb, f_x(xb))
        disp('黄金分割法一维搜索完毕.')
        break
    elseif f1<f2
        disp('f1)
        b = x2;
        x2 = x1;
        f2 = f1;
        flag = 1;
    elseif f1==f2
        disp('f1=f2')
        a = x1;
        b = x2;
        flag = 0;
    elseif f1>f2
        disp('f1>f2')
        a = x1;
        x1 = x2;
        f1 = f2;
        flag = 2;
    end
    cnt = cnt + 1;
end
pause(0.5)
stem(xb, f_x(xb), 'r', 'linewidth', 2)
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

如果不想等待动画加载，ctrl+f, ctrl+r 把pause批量注释即可；

《最优化方法》教材上写成表的答案：

黄金分割法的一些性质

1、x1 = a+b-x2;
2、下一次迭代的区间长度是上一个区间长度的0.618倍；
3、如果f1 同理如果f1>f2，则上一次迭代的x2, f2传给下一次迭代的x1, f1；
4、迭代次数和求解精度取决于终止条件$|b-a|<ϵ$中$ϵ$的大小。