2024.5.22（周三）

function [x, val, k] = dampnm(fun, gfun, Hess, x0)

% 功能: 用阻尼牛顿法求解无约束优化问题：min f(x)

% 输入: x0 是初始点, fun, gfun, Hess 分别是目标函数和梯度 Hess 阵函数

% 输出: x, val 分别是近似最优解和近似最优值, k 是迭代次数

maxk = 5000;

rho = 0.5;

sigma = 0.4;

k = 0;

epsilon = 1e-6;

while (k < maxk)

gk = feval(gfun, x0);

Gk = feval(Hess, x0);

% 使用更为稳定的解法

dk = -Gk \ gk;

if (norm(dk) < epsilon)

break;

end

m = 0;

while (m < 20)

if (feval(fun, x0 + rho^m * dk) <= feval(fun, x0) + sigma * rho^m * gk' * dk)

break;

end

m = m + 1;

end

x0 = x0 + rho^m * dk;

k = k + 1;

end

x = x0;

val = feval(fun, x);

end

% 定义目标函数

function y = fun(x)

y = (x(1) + 10 * x(2))^2 + 5 * (x(3) - x(4))^2 + (x(2) - 2 * x(3))^4 + 10 * (x(1) - x(4))^4;

end

% 定义梯度函数

function g = gfun(x)

g = [2 * (x(1) + 10 * x(2)) + 40 * (x(1) - x(4))^3;

20 * (x(1) + 10 * x(2)) + 4 * (x(2) - 2 * x(3))^3;

10 * (x(3) - x(4)) - 8 * (x(2) - 2 * x(3))^3;

-10 * (x(3) - x(4)) - 40 * (x(1) - x(4))^3];

end

% 定义 Hessian 矩阵函数

function He = Hess(x)

He = [2 + 120 * (x(1) - x(4))^2, 20, 0, -120 * (x(1) - x(4))^2;

20, 200 + 12 * (x(2) - 2 * x(3))^2, -24 * (x(2) - 2 * x(3))^2, 0;

0, -24 * (x(2) - 2 * x(3))^2, 10 + 48 * (x(2) - 2 * x(3))^2, -10;

-120 * (x(1) - x(4))^2, 0, -10, 10 + 120 * (x(1) - x(4))^2];

end

% 调用示例

x0 = [0; 1; 2; 2]; % 初始点

[x, val, k] = dampnm(@fun, @gfun, @Hess, x0);

% 输出结果

disp('近似最优解:');

disp(x);

disp('近似最优值:');

disp(val);

disp('迭代次数:');

disp(k);