二次优化问题dfp_优化方法·无约束非线性优化

就是求一个函数的最小值问题。并且计算机不适合用解析解。

我们可以从一个初始点出发，一步一步走到期望的最优点，这就要求求出每一步往哪里走，走多长，从而把这种优化问题分解成求方向和求步长两个问题，求步长的时候，可以认为函数只有步长一个自变量，所以又称为一维搜索。

附上我写的如下一段代码：求方向可以选择4种算法之一，求步长则用的是不精确一维搜索方法一

%主函数集成了4种算法，所用的初始点可由用户自定义，用于求解无约束优化问题。
%目标函数定义在程序末尾的fun函数中，其梯度定义在gfun函数中，海塞矩阵（阻尼牛顿法需要）定义在d2fun函数中。
function Ti1
clear,clc;
x0=input('函数和梯度表达式在本程序末尾设置，现在请确定一个迭代初始点，以列向量表示：nx0='); %用户可自定迭代初始点，增强了算法的灵活性
disp('请输入数字选择一种算法：1-最速下降法；2-阻尼牛顿法；3-共轭梯度法；4或其他输入-BFGS算法');
choice=input('选择...'); %choice是特定的选择变量，能够决定程序执行哪个算法
kmax=10000;   %最大迭代次数，避免算法无法收敛而一直进行
lambda=0;   k=0;    epsilon=1e-5;   g=gfun(x0); %不同的方法在迭代前都需要的参数，其中lambda是迭代步长，epsilon是精度
switch choice
    case 1
        disp('选择了最速下降法'); %提示选择的算法，下同
    case 2
        disp('选择了阻尼牛顿法');
        d2f=d2fun(x0); %阻尼牛顿法迭代前需要补充的参数，d2f即为海塞矩阵
    case 3
        disp('选择了共轭梯度法');
        n=length(x0);   s=zeros(n,1);   mu=0; %共轭梯度法迭代前需要补充的参数，mu是确定下个搜索方向最关键的参数μ
    otherwise
        disp('选择了BFGS算法');
        n=length(x0);   H=eye(n);   mu=0;   g=gfun(x0); %BFGS算法迭代前需要补充的参数，H和mu分别为公式中的Hk和μk
end
while(k<kmax)
    if(norm(g)<epsilon) %迭代停止准则
        break;
    end
    switch choice %根据选择的算法确定搜索方向
        case 1
            s=-g; %最速下降法的搜索方向
        case 2
            s=-inv(d2f)*g; %阻尼牛顿法的搜索方向
        case 3
            s=-g+mu*s; %共轭梯度法的搜索方向
        otherwise
            s=-H*g; %BFGS算法的搜索方向
    end
    %求步长采用不精确一维搜索算法一，相比直接求解驻点方程和精确一维搜索更简单，运算也更迅速
    c1=0.1;c2=0.5;a=0;b=inf;t=1;j=0;jmax=100; %一维搜索迭代参数设置
    %其中t为存储步长的变量，最大步数jmax的设置是避免一维搜索无法跳出循环而一直占用计算资源
    while(j<jmax)
        if(fun(x0)-fun(x0+t*s)+c1*t*(gfun(x0))'*s>=0) %采用Wolfe-Powell准则
            if((gfun(x0+t*s)-c2*gfun(x0))'*s>=0)
                break;
            else
                j=j+1;a=t;t=min(2*t,(t+b)/2);
            end
        else
            j=j+1;b=t;t=(t+a)/2;
        end
    end
    %以上为不精确一维搜索算法一过程
    lambda=t; %步长的确定
    x=x0+lambda*s; %下一个迭代点的确定
    %对最速下降法，已经可以进行下一步迭代，但对另外三种方法，需要更新参数
    switch choice
        case 1 %对最速下降法
            x0=x; %更新迭代点的值，下同
            g=gfun(x0); %求解新的梯度，下同
            k=k+1; %进行下一次迭代，下同
        case 2 %对阻尼牛顿法
            x0=x;
            g=gfun(x0);
            k=k+1;
            d2f=d2fun(x0); %阻尼牛顿法中海塞矩阵的更新
        case 3 %对共轭梯度法
            mu=norm(gfun(x))^2/(norm(g)^2); %共轭梯度法中μ的更新
            x0=x;
            g=gfun(x0);
            k=k+1;
        otherwise %对BFGS算法
            dx=x-x0;dg=gfun(x)-g;
            mu=1+(dg'*H*dg)/(dx'*dg); %BFGS算法中μ的更新
            H=H+(mu*(dx*dx')-H*dg*dx'-dx*dg'*H)/(dx'*dg); %BFGS算法中H的更新
            x0=x;
            g=gfun(x0);
            k=k+1;
    end
end
if(k==kmax) %迭代超过最大次数仍未收敛，说明选择的算法收敛性不好