最优化算法中三点二次插值法求最优解matlab程序

function [xopt, fopt, iter] = quadratic_interpolation(f, x0, eps)
% f: 目标函数句柄
% x0: 初始点
% eps: 迭代停止条件
% xopt: 最优解
% fopt: 最优值
% iter: 迭代次数

x1 = x0;
f1 = feval(f, x1);
iter = 0;

while true
iter = iter + 1;
% 计算梯度
g = gradient(f, x1);
% 计算海森矩阵
H = hessian(f, x1);
% 求解方程 H d = -g
d = -H \ g;
% 一维搜索
[a, b] = bracket(f, x1, d);
x2 = golden_section(f, a, b);
% 判断是否收敛
if norm(x2 - x1) < eps
xopt = x2;
fopt = feval(f, xopt);
break;
end
% 更新x1和f1
f2 = feval(f, x2);
if f2 < f1
x1 = x2;
f1 = f2;
else
% 二次插值
alpha = quadratic_interpolation_1d(f, x1, f1, x2, f2);
x3 = x1 + alpha * d;
% 更新x1和f1
x2 = x1;
x1 = x3;
f1 = feval(f, x1);
end
end

end

function alpha = quadratic_interpolation_1d(f, x1, f1, x2, f2)
% 一维二次插值
% f: 目标函数句柄
% x1: 已知点1
% f1: 已知点1的函数值
% x2: 已知点2
% f2: 已知点2的函数值
% alpha: 最优步长

% 计算一维向量
d = x2 - x1;
% 计算一维向量的长度
L = norm(d);
% 计算一维向量的单位向量
u = d / L;
% 计算一维向量的中点
x0 = (x1 + x2) / 2;
% 计算一维向量的中点函数值
f0 = feval(f, x0);
% 计算二次插值系数a,b,c
a = (f1 - 2*f0 + f2) / L^2;
b = (f2 - f1) / L - a*L;
c = f1;
% 求解极小点
alpha = -b / (2*a);
% 限制步长范围
alpha = max(alpha, 0);
alpha = min(alpha, L);

end

function [a, b] = bracket(f, x, d)
% 一维搜索区间搜索
% f: 目标函数句柄
% x: 当前点
% d: 搜索方向
% a: 搜索区间左端点
% b: 搜索区间右端点

% 初始步长
alpha = 0.1;
% 初始函数值
f1 = feval(f, x);
% 向搜索方向移动一定距离
x2 = x + alpha*d;
% 向反方向移动一定距离
x0 = x - alpha*d;
% 判断函数值变化情况
f2 = feval(f, x2);
if f2 > f1
% 如果函数值增加，则继续向搜索方向移动
while true
alpha = alpha * 2;
x2 = x + alpha*d;
f2 = feval(f, x2);
if f2 < f1
break;
end
end
a = x;
b = x2;
elseif f2 < f1
% 如果函数值减少，则继续向反方向移动
while true
alpha = alpha * 2;
x0 = x - alpha*d;
f0 = feval(f, x0);
if f0 < f1
break;
end
end
a = x0;
b = x;
else
% 如果函数值没有变化，则向两个方向移动
while true
alpha = alpha * 2;
x2 = x + alpha*d;
f2 = feval(f, x2);
x0 = x - alpha*d;
f0 = feval(f, x0);
if f0 < f1 || f2 < f1
break;
end
end
a = x0;
b = x2;
end

end

function xopt = golden_section(f, a, b)
% 一维搜索的黄金分割法
% f: 目标函数句柄
% a: 左端点
% b: 右端点
% xopt: 最优解

% 黄金分割比例
alpha = (sqrt(5) - 1) / 2;
% 初始区间长度
L = b - a;
% 初始内点
x1 = a + alpha * L;
x2 = b - alpha * L;
% 初始函数值
f1 = feval(f, x1);
f2 = feval(f, x2);
% 迭代
while L > eps
if f1 > f2
% 右侧区间更优
a = x1;
x1 = x2;
f1 = f2;
L = b - a;
x2 = b - alpha * L;
f2 = feval(f, x2);
else
% 左侧区间更优
b = x2;
x2 = x1;
f2 = f1;
L = b - a;
x1 = a + alpha * L;
f1 = feval(f, x1);
end
end

% 返回内点
xopt = (a + b) / 2;

end

function g = gradient(f, x)
% 计算目标函数的梯度
% f: 目标函数句柄
% x: 自变量向量
% g: 梯度向量

h = 1e-6;
n = length(x);
g = zeros(n, 1);
for i = 1:n
dx = zeros(n, 1);
dx(i) = h;
g(i) = (feval(f, x+dx) - feval(f, x-dx)) / (2*h);
end

end

function H = hessian(f, x)
% 计算目标函数的海森矩阵
% f: 目标函数句柄
% x: 自变量向量
% H: 海森矩阵

h = 1e-6;
n = length(x);
H = zeros(n, n);
for i = 1:n
dx = zeros(n, 1);
dx(i) = h;
for j = i:n
dy = zeros(n, 1);
dy(j) = h;
H(i,j) = (feval(f, x+dx+dy) - feval(f, x+dx-dy) - feval(f, x-dx+dy) + feval(f, x-dx-dy)) / (4*h^2);
H(j,i) = H(i,j);
end
end

end