matlab diff函数用法_工程优化设计与Matlab实现——无约束问题的导数解法（变尺度法）...

尺有所短，寸有所长。 ——屈原《楚辞·卜居》

取其精华，去其糟粕。 ——鲁迅《拿来主义》

变尺度法

如果说我们想要在距离极值点比较远的起点处，像最速下降法那样小计算量且速度不算太慢迭代，在距离极值点比较近时像牛顿法那样快速收敛，该怎么办呢？

让我们回顾一下这两种方法的迭代思路：

最速下降法迭代：

阻尼牛顿法迭代：

可以看出，迭代的差别在于方向的不同，阻尼牛顿方法的方向比最速下降法多了一个Hessian矩阵的逆。

那我们是否可以构造一个矩阵

来取代阻尼牛顿法中Hessian矩阵的位置，并且这个矩阵可以随着迭代点位置作出改变，在起始点时

，在极值点时

。

答案当然是肯定的，这就是像“初一到十五的月亮”一样的变尺度法，其迭代公式如下：

那么问题就来了，这个矩阵

怎么构造呢？别急，看Point 1 。

Point 1 变尺度矩阵

的建立

一方面，构造近似矩阵

的迭代格式：

其中

与

是n阶对称阵，则修正矩阵

也一定是n阶对称阵

所以

的格式可以设为

故有

另一方面，对目标函数在

点泰勒二阶展开，有：

求一阶导数有：

将迭代点

代入上式，有：

移项得：

此时，令

则有

两个方面都准备好后，我们用

代替

，就有：

代入

，并移项得：

观察上式，令

，得：

比较左右两边系数，可得：

将

、

、

、

代入

就得到变尺度矩阵：

Point 2 变尺度法的步长

变尺度矩阵构造完成，搜索方向也就确定了，然后又是一维搜索求最优步长了，此处省略N个字。如有疑问，请往前参考：工程优化设计与Matlab实现——一维搜索方法（总结）中提到的一维搜索解决多维问题，工程优化设计与Matlab实现——无约束问题的导数解法（一）中最速下降法求最优步长，工程优化设计与Matlab实现——无约束问题的导数解法（三）中阻尼牛顿法求最优步长。

变尺度法流程

Point 3 栗子

利用变尺度法求目标函数

，在初始值为

时的极小值点和极小值。

主程序：

clc

目标函数定义：

function

变尺度法函数定义：

使用求导的方式，获得步长

%-----------变尺度法---主函数-----------%
function [BianChiDuFa_x, BianChiDuFa_xf, BianChiDuFa_n] = BianChiDuFa(e, x0)
    syms t1 t2 t3 t4 a
    x = [t1;t2];
    f = func(x);
    
    %求x0处jacobian矩阵的值
    yf = jacobian(f)';       %求jacobian矩阵的转置（梯度）

    n = 0;                     %迭代次数统计
    A0 = 1;
    k = 1;                      %迭代开
    while k == 1
        n = n + 1;                     %迭代次数统计     
        yf_x0 = subs(yf, x, x0);     %将x0代入，得x0处的梯度
        d0 = -A0 * yf_x0;                                   %x0迭代方向

        %求最优步长a
        x1 = x0 + a * d0;                 %设最优步长为未知数a，求得x1
        f_x1 = subs(f, x, x1);            %将x1代入原函数
        diff_f_a =diff(f_x1, a);          %将x1代入原函数后得到的含有a的函数，对a求偏导 
        a_n = solve(diff_f_a, a);         %偏导数等于0，求得使函数取极值的a

        %将求得的a代回，求得具体的x1，f（x1）,和x1处的梯度
        x1 = subs(x1, a, a_n);             %x1的具体值
        f_x1 = subs(f, x, x1);             %x1的函数值
        yf_x1 = subs(yf, x, x1);           %x1的梯度值
        
        f_x0 = subs(f, x, x0);            %将x0代入原函数
        
        %判断迭代是否终止
        if abs(f_x1-f_x0) <= e
            break
        else
            Dx0 = x1 - x0;
            Dyf0 = yf_x1 - yf_x0;
            DA0 = (Dx0 * Dx0') / (Dx0' * Dyf0) - (A0 * Dyf0 * Dyf0' * A0) / (Dyf0' * A0 * Dyf0);
            A1 = A0 + DA0;
            x0 = x1;
            A0 = A1;
        end
    end

    BianChiDuFa_x = vpa(x1);
    BianChiDuFa_xf = vpa(f_x1);
    BianChiDuFa_n = n;

使用一维搜索方法，获得步长

%-----------变尺度法---主函数-----------%
function [BianChiDuFa_x, BianChiDuFa_xf, BianChiDuFa_n] = BianChiDuFa2(e, x0)
    syms t1 t2 t3 t4 a
    x = [t1;t2];
    f = func(x);
    
    %求x0处jacobian矩阵的值
    yf = jacobian(f)';       %求jacobian矩阵的转置（梯度）

    n = 0;                     %迭代次数统计
    A0 = 1;
    k = 1;                      %迭代开
    while k == 1
        n = n + 1;                     %迭代次数统计     
        %fprintf('迭代次数：%fn',n)
        
        f_x0 = subs(f, x, x0);       %将x0代入，得x0处的函数值
        yf_x0 = subs(yf, x, x0);     %将x0代入，得x0处的梯度
        d0 = -A0 * yf_x0;            %x0迭代方向

        %一维搜索方法求最优步长（求x1）
        [JT_a, JT_b, JT_c] = JinTuiFa(x0, 0.1, d0);    %进退法初始搜索步长，设为0.1    
        [x1, f_x1, HJFGF_n] = HuangJinFenGeFa(JT_a, JT_c, e);   %黄金分割法
        yf_x1 = subs(yf, x, x1);           %x1的梯度值
        
        if abs(f_x1-f_x0) <= e
            break
        else
            Dx0 = x1 - x0;
            Dyf0 = yf_x1 - yf_x0;
            DA0 = (Dx0 * Dx0') / (Dx0' * Dyf0) - (A0 * Dyf0 * Dyf0' * A0) / (Dyf0' * A0 * Dyf0);
            A1 = A0 + DA0;
            x0 = x1;
            A0 = A1;
        end
    end

    BianChiDuFa_x = vpa(x1);
    BianChiDuFa_xf = vpa(f_x1);
    BianChiDuFa_n = n;

    

进退法程序：工程优化设计与Matlab实现——一维搜索方法（一）

黄金分割法程序：工程优化设计与Matlab实现——一维搜索方法（二）