求最优化问题，用matlab

Question

f = x(1)^2 + x(2)^2 - 3*x(1) - x(1)*x(2) + 3;
取x0 = [0; 0]; 用BFGS和DFP算法编写程序，不用函数里的库

Answer 1

引自免费微信小程序：皆我百晓生

为了实现BFGS（Broyden–Fletcher–Goldfarb–Shanno）和DFP（Davidon–Fletcher–Powell）算法，我们需要手动编写相应的迭代过程。这里是一个使用Matlab编写的通用优化器模板，适用于这两种方法。我们将分别实现这两个算法，并应用到给定的问题上。

首先，定义目标函数f(x)和梯度grad_f(x)：

function f = objective_function(x)
    f = x(1)^2 + x(2)^2 - 3*x(1) - x(1)*x(2) + 3;
end

function grad_f = gradient_function(x)
    grad_f = [2*x(1) - 3 - x(2); 2*x(2) - x(1)];
end

接下来，实现BFGS算法：

function [x_opt, f_val] = bfgs_algorithm(f, grad_f, x0, options)
    % 初始化变量
    n = length(x0);
    H = eye(n); % 初始Hessian矩阵为单位矩阵
    rho = ones(n, 1); % 用于计算rho值的向量
    y_k = zeros(n, 1); % 梯度变化向量
    s_k = zeros(n, 1); % 步长向量
    x_k = x0; % 初始解
    f_k = feval(f, x_k); % 初始函数值
    
    % 迭代过程
    for k = 1:options.max_iterations
        % 计算梯度
        g_k = feval(grad_f, x_k);
        
        % 如果满足停止条件，跳出循环
        if norm(g_k) < options.tolerance
            break;
        end
        
        % 更新y_k和s_k
        y_k = g_k - rho(k-1) * s_k;
        rho(k) = 1 / dot(y_k, g_k);
        s_k = - rho(k) * y_k;
        
        % 预测步长
        alpha_k = 1;
        while true
            x_next = x_k + alpha_k * s_k;
            f_next = feval(f, x_next);
            
            % Armijo准则检查
            if f_next <= f_k + options.armijo_factor * alpha_k * dot(s_k, g_k)
                break;
            else
                alpha_k = 0.5 * alpha_k;
            end
        end
        
        % 更新x_k和f_k
        x_k = x_next;
        f_k = f_next;
        
        % BFGS更新Hessian矩阵
        sy = rho(k) * s_k;
        hy = rho(k) * y_k;
        H = H + (sy * sy' - H * s_k * sy') / dot(sy, y_k) - (H * hy * sy' + hy * sy' * H) / (dot(hy, sy));
        H = H + (hy * hy' - hy * sy' * inv(H) * sy * hy') / (dot(hy, sy) - dot(hy, s_k) * dot(sy, y_k)) / dot(sy, y_k);
    end
    
    x_opt = x_k;
    f_val = f_k;
end

然后，实现DFP算法：

function [x_opt, f_val] = dfp_algorithm(f, grad_f, x0, options)
    % 初始化变量
    n = length(x0);
    H = eye(n); % 初始Hessian矩阵为单位矩阵
    y_k = zeros(n, 1); % 梯度变化向量
    s_k = zeros(n, 1); % 步长向量
    x_k = x0; % 初始解
    f_k = feval(f, x_k); % 初始函数值
    
    % 迭代过程
    for k = 1:options.max_iterations
        % 计算梯度
        g_k = feval(grad_f, x_k);
        
        % 如果满足停止条件，跳出循环
        if norm(g_k) < options.tolerance
            break;
        end
        
        % 更新y_k和s_k
        y_k = g_k - H * s_k;
        s_k = - y_k;
        
        % 预测步长
        alpha_k = 1;
        while true
            x_next = x_k + alpha_k * s_k;
            f_next = feval(f, x_next);
            
            % Armijo准则检查
            if f_next <= f_k + options.armijo_factor * alpha_k * dot(s_k, g_k)
                break;
            else
                alpha_k = 0.5 * alpha_k;
            end
        end
        
        % 更新x_k和f_k
        x_k = x_next;
        f_k = f_next;
        
        % DFP更新Hessian矩阵
        sy = y_k - H * s_k;
        H = H + (sy * sy' - H * s_k * sy') / dot(sy, y_k);
    end
    
    x_opt = x_k;
    f_val = f_k;
end

最后，调用上述函数并设置选项：

% 设置初始点和选项
x0 = [0; 0];
options = struct('max_iterations', 1000, 'tolerance', 1e-6, 'armijo_factor', 1e-4);

% 使用BFGS算法
[x_bfgs, f_bfgs] = bfgs_algorithm(@objective_function, @gradient_function, x0, options);

% 使用DFP算法
[x_dfp, f_dfp] = dfp_algorithm(@objective_function, @gradient_function, x0, options);

% 输出结果
fprintf('Optimal solution using BFGS: %f, %f\n', x_bfgs(1), x_bfgs(2))
fprintf('Objective function value at the optimal solution using BFGS: %f\n', f_bfgs)

fprintf('\nOptimal solution using DFP: %f, %f\n', x_dfp(1), x_dfp(2))
fprintf('Objective function value at the optimal solution using DFP: %f\n', f_dfp)

这段代码将分别用BFGS和DFP算法解决给定的优化问题，并输出最优解和对应的函数值。注意，实际运行时可能需要根据实际情况调整选项参数，例如最大迭代次数、收敛阈值等。