5.5

题目：

上机实验四：共轭梯度法程序设计

1、基本要求

掌握共轭梯度法的基本思想及其迭代步骤；学会运用MATLAB编程实现常用优化算法；能够正确处理实验数据和分析实验结果及调试程序。

2、主要内容

（1）求解无约束优化问题：

（2）终止准则取；

（3）完成FR共轭梯度法的MATLAB编程、调试；

（4）选取几个与实验二实验三中相同的初始点，并给出相关实验结果的对比及分析（从最优解、最优值、收敛速度（迭代次数）等方面进行比较）；

（5）按照模板撰写实验报告，要求规范整洁。

3、操作要点

（1）通过编程实现FR共轭梯度法；

（2）使用MTALAB调试程序，并将实验结果保存到文件中；

（3）撰写实验报告。

4、主要仪器设备

微机及Matlab软件

代码：

1. 算法步骤

初始化

定义目标函数 func，它是一个字符串表达式，描述了优化问题的目标。

初始化起始点 x0，它是一个向量，表示在搜索空间中的起始位置。

设置容忍度 tol，用于确定何时停止迭代（当梯度的范数小于此值时）。

设置最大迭代次数 max_iter，以防止无限循环。

定义变量和函数

使用 length(x0) 确定变量的数量 n。

使用符号变量 vars 定义符号表达式中的变量。

将目标函数字符串 func 转换为符号表达式 funf_sym。

计算目标函数的梯度 gradf_sym 和海森矩阵 hessf_sym（尽管FR方法通常不需要海森矩阵，但这里被计算了）。

将符号表达式转换为MATLAB函数句柄 funf、gradf 和 hessf，以便在数值计算中使用。

迭代过程

初始化当前点 x 为起始点 x0。

计算当前点的梯度 g。

初始化搜索方向 d 为负梯度方向 -g。

初始化迭代计数器 k 为0。

主循环

当梯度的范数大于 tol 且迭代次数小于 max_iter 时，进行以下迭代：

计算搜索方向 d 与海森矩阵 Hd 的乘积（尽管FR方法不使用这个值来计算步长，但这里计算了）。

使用FR公式计算步长 alpha：alpha = -(g' * d) / (d' * Hd)（注意：这里通常使用线搜索来找到最佳步长，而不是直接用这个公式）。

更新当前点 x：x = x + alpha * d。

计算新点的梯度 g_new。

使用FR公式更新搜索方向 d：d = -g_new + beta * d，其中 beta = (g_new' * g_new) / (g' * g)。

更新梯度 g 为 g_new。

递增迭代计数器 k。

输出

迭代完成后，输出最优解 x 和迭代次数 k

2. 代码

function [x_min, f_min] = lab04(f, grad_f, x_initial, tol)

% 定义目标函数和梯度函数

% 注意：这里的f和grad_f应该是作为输入参数传递进来的函数句柄

% 但为了演示，我直接定义了函数内容。在实际应用中，应由调用者提供这些函数句柄。

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

grad_f = @(x) [2*(x(1)+10*x(2)) + 40*(-x(1)+x(4)); ...

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

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

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

% 初始化参数

x = x_initial;

r = -grad_f(x);

p = r;

rsold = r'*r;

tolerance = tol; % 设置容差

maxIter = 1000; % 最大迭代次数

iter = 0;

while (rsold > tolerance) && (iter < maxIter)

    iter = iter + 1;

    Ap = grad_f(x + p); % 计算Ap

    alpha = rsold / (p'*Ap); % 计算步长

    x = x + alpha * p; % 更新x

    r = r - alpha * Ap; % 更新残差

    rsnew = r'*r; % 计算新的rs

    if sqrt(rsnew) < tolerance

        break;

    end

    beta = rsnew / rsold; % Fletcher-Reeves公式计算beta

    p = r + beta * p; % 更新搜索方向

    rsold = rsnew; % 更新rsold

end

% 循环结束后，直接使用最终的x和f(x)作为结果

x_min = x;

f_min = f(x);

end