共轭梯度法程序设计

共轭梯度法程序设计

一、实验目的

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

二、实验内容

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

（2）终止准则取；

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

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

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

三、算法步骤、代码、及结果

1. 算法步骤

初始化参数后，开始迭代，循环迭代，直到达到最大迭代次数或满足终止条件。计算梯度g=feval(gfun,x0);计算搜索方向、当前迭代步的周期索引。根据FR公式计算新的搜索方向。如果是周期的第一步，则方向为负梯度方向，否则计算共轭方向。检查终止条件使用Armijo准则确定步长。更新变量，最后返回结果。

2. 代码

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

    maxk = 5000;

    rho = 1e-6;  % 步长缩减因子

    sigma = 0.01;

    k = 0;

    epsilon = 1e-6;

    n = length(x0);

    g0 = feval(gfun, x0);

    d0 = -g0;

    while (k < maxk)

        g = feval(gfun, x0);

       

        % 限制梯度值避免过大

        g = max(min(g, 1e5), -1e5);

       

        itern = k - (n + 1) * floor(k / (n + 1));

        itern = itern + 1;

        if (itern == 1)

            d = -g;

        else

            beta = (g' * g) / (g0' * g0);

            d = -g + beta * d0;

            gd = g' * d;

            if (gd >= 0.0)

                d = -g;

            end

        end

        if (norm(g) < epsilon), break; end

        m = 0; mk = 0;

        while (m < 20)

            if (feval(fun, x0 + rho^m * d) < feval(fun, x0) + sigma * rho^m * g' * d)

                mk = m; break;

            end

            m = m + 1;

        end

        x0 = x0 + rho^mk * d;

        val = feval(fun, x0);

        g0 = g;

        d0 = d;

        k = k + 1;

       

        % 调试信息输出

        fprintf('Iteration: %d, x: [%s], f(x): %.6f, grad: [%s], step size: %.6f\n', ...

                k, num2str(x0'), val, num2str(g'), rho^mk);

       

        % 数值检查

        if any(isnan(x0)) || any(isinf(x0))

            fprintf('Numerical issue detected at iteration %d\n', k);

            break;

        end

    end

    x = x0;

    val = feval(fun, x);

end

 

% 目标函数

function f = fun(x)

    f = (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

 

 

x0 = [1; 1; 1; 1];

[x_opt, fval_opt, iter] = frcg(@fun, @gfun, x0);

fprintf('Optimal point: [%s]\n', num2str(x_opt'));

fprintf('Optimal value: %.6f\n', fval_opt);

fprintf('Iterations: %d\n', iter);

 

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*((3) - x(4)) + 8*(x(2) - 2*x(3))^3;

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

    ];

end

 

function f = fun(x)

    f = (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