前言
黄金分割法是一种优化算法,用于求解函数的极小值。它基于黄金比例,通过不断缩小搜索范围来逼近极小值点。该方法可以在一定程度上避免了传统的梯度下降法和牛顿法可能陷入局部最优解的问题。
具体步骤如下:
- 确定初始搜索范围 [ a , b ] [a, b] [a,b],使得函数在这个范围内有极小值点。
- 计算两个内点 α 1 \alpha_1 α1和 α 2 \alpha_2 α2,使得 α 1 \alpha_1 α1和 α 2 \alpha_2 α2分别与 a a a和 b b b的距离满足黄金比例。
- 计算函数在 α 1 \alpha_1 α1和 α 2 \alpha_2 α2的取值,选择较小的值所对应的点作为新的搜索范围的端点。
- 不断重复步骤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.