MATLAB实现BFGS算法，用wolfe选择步长

BFGS 算法简介

BFGS（Broyden-Fletcher-Goldfarb-Shanno）算法是一种用于无约束优化问题的准牛顿法。它通过构建和更新近似的海森矩阵（Hessian Matrix）来加速优化过程。BFGS 是一种常用且有效的优化算法，尤其适用于大规模问题。

步长选择：Wolfe 条件

为了选择合适的步长，我们可以使用 Wolfe 条件。这些条件包括：

Armijo 条件（或称为充分下降条件）：确保步长不会过大，从而使得目标函数值显著减少。
曲率条件：确保步长不会过小，使得导数的减少足够大。

软约束处理方法

在处理有约束的优化问题时，可以使用软约束方法，例如罚函数法。通过在目标函数中添加一个惩罚项来处理约束条件。

MATLAB 代码示例

以下是一个可以在 MATLAB 中直接运行的代码示例，展示如何使用 BFGS 算法解决实际问题，包括 Wolfe 条件的步长选择和软约束处理。

function bfgs_optimization_example()
    % 初始值
    x0 = [1.2; 1.2];
    
    % 目标函数
    objective = @(x) (1 - x(1))^2 + 100 * (x(2) - x(1)^2)^2;
    
    % 约束函数
    constraints = @(x) x(1) + x(2) - 1.5;
    
    % 罚函数系数
    penalty_coeff = 1000;
    
    % 优化
    [x_opt, f_opt] = bfgs(objective, constraints, x0, penalty_coeff);
    
    % 结果显示
    disp('Optimal solution:');
    disp(x_opt);
    disp('Optimal objective function value:');
    disp(f_opt);
end

function [x_opt, f_opt] = bfgs(objective, constraints, x0, penalty_coeff)
    % 初始值
    x = x0;
    n = length(x);
    H = eye(n);  % 初始近似海森矩阵
    
    % 迭代参数
    max_iter = 100;
    tol = 1e-6;
    
    % 记录目标函数值
    f_values = zeros(max_iter, 1);
    
    % 创建绘图窗口
    figure;
    hold on;
    xlabel('Iteration');
    ylabel('Objective Function Value');
    title('Objective Function Value over Iterations');
    
    for iter = 1:max_iter
        % 计算目标函数和梯度
        [f, grad] = penalized_objective(objective, constraints, x, penalty_coeff);
        
        % 记录当前目标函数值
        f_values(iter) = f;
        
        % 实时更新绘图
        plot(1:iter, f_values(1:iter), 'b-');
        drawnow;
        
        % 计算搜索方向
        p = -H * grad;
        
        % Wolfe 条件选择步长
        alpha = line_search_wolfe(objective, constraints, x, p, grad, penalty_coeff);
        
        % 更新
        x_new = x + alpha * p;
        
        % 计算新的梯度
        [~, grad_new] = penalized_objective(objective, constraints, x_new, penalty_coeff);
        
        % BFGS 公式更新 H
        s = x_new - x;
        y = grad_new - grad;
        rho = 1 / (y' * s);
        H = (eye(n) - rho * s * y') * H * (eye(n) - rho * y * s') + rho * s * s';
        
        % 检查收敛
        if norm(x_new - x) < tol
            break;
        end
        
        x = x_new;
    end
    
    x_opt = x;
    f_opt = objective(x);
end

function [f, grad] = penalized_objective(objective, constraints, x, penalty_coeff)
    % 原始目标函数值和梯度
    f = objective(x);
    grad = numerical_gradient(objective, x);
    
    % 约束惩罚
    c = constraints(x);
    f = f + penalty_coeff * max(0, c)^2;
    
    % 约束梯度惩罚
    if c > 0
        grad = grad + 2 * penalty_coeff * c * numerical_gradient(@(x) constraints(x), x);
    end
end

function grad = numerical_gradient(f, x)
    % 数值梯度计算
    epsilon = 1e-8;
    grad = zeros(length(x), 1);
    
    for i = 1:length(x)
        x1 = x;
        x1(i) = x1(i) + epsilon;
        grad(i) = (f(x1) - f(x)) / epsilon;
    end
end

function alpha = line_search_wolfe(objective, constraints, x, p, grad, penalty_coeff)
    % Wolfe 条件的步长选择
    alpha = 1;
    c1 = 1e-4;
    c2 = 0.9;
    
    phi = @(a) penalized_objective(objective, constraints, x + a * p, penalty_coeff);
    [f0, ~] = phi(0);
    [f_alpha, grad_alpha] = phi(alpha);
    
    while f_alpha > f0 + c1 * alpha * grad' * p || grad_alpha' * p < c2 * grad' * p
        alpha = alpha / 2;
        [f_alpha, grad_alpha] = phi(alpha);
    end
end