数学建模竞赛算法全集附matlab实现（建议收藏）

5.无约束优化算法

原始的牛顿迭代法的缺点：

1）初值点选择（敏感）：基本牛顿法初始点需要足够靠近极小点，否则有可能导致算法不收敛。

2）Hessen矩阵正定

代码：

unction [x, val, k] = niudun1(x0)
    %基本牛顿法求函数最小值
    %x是最小值点，val是最小值，k是迭代次数
    epsilon = 1e-5; %终止误差值
    k = 0;  %迭代次数
    %如果梯度的范数小于终止误差值则停止迭代
    while norm(g(x0)) > epsilon
        d = -G(x0)\g(x0);
        x0 = x0 + d';
        k = k+1;
    end
    x = x0;
    val = f(x0);
end
 
%%原函数
function y = f(x)
    y = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;
end
 
%%原函数的梯度,要弄成列向量的形式
function y = g(x)
    y(1,1) = 400*x(1)*(x(1)^2 - x(2)) + 2*(x(1) - 1);
    y(2,1) = -200*(x(1)^2 - x(2));
end
 
%%原函数的Hessen矩阵，注意矩阵每个元素的值
function y = G(x)
    y(1, 1) = 1200*x(1)^2 - 400*x(2) + 2;
    y(1, 2) = -400*x(1);
    y(2, 1) = -400*x(1);
    y(2, 2) = 200;
end

（2）阻尼牛顿法（

优点：阻尼牛顿法在基本牛顿法的基础上增加对步长的线搜索，克服了基本牛顿法对于初始点敏感的问题。

算法：

代码：（对于x0要输入列向量）

function [x,val,k]=damp_Newton(fun,gfun, Hess,x0)
%功能: 用阻尼牛顿法求解无约束问题:  min f(x)
%输入: x0是初始点, fun, gfun, Hess 分别是求
%         目标函数,梯度,Hesse 阵的函数
%输出:  x, val分别是近似最优点和最优值,  k是迭代次数.
maxk=10000;   %给出最大迭代次数
rho=0.55;sigma=0.5;
k=0;  epsilon=1e-5;
tic
while(k<maxk)
    gk=feval(gfun,x0); %计算梯度
    Gk=feval(Hess,x0);  %计算Hesse阵
    dk=-Gk\gk; %解方程组Gk*dk=-gk, 计算搜索方向
    if(norm(gk)<epsilon), break; end  %检验终止准则
    m=0; mk=0;
    while(m<20)   % 用Armijo搜索求步长 
        if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk)
            mk=m; break;
        end
        m=m+1;
    end
    x0=x0+rho^mk*dk;
    k=k+1
end
x=x0;
val=feval(fun,x); 
toc
%gval=norm(gfun(x));
%%%% 目标函数 %%%%%%%%%
%function f=fun(x)
%f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;
%%%% 梯度 %%%%%%%%%%%%%%%%%%%
%function g=gfun(x)
%g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1), -200*(x(1)^2-x(2))]';
%%%% Hesse 阵 %%%%%%%%%%%%%%%%%%%
%function He=Hess(x)
%n=length(x);
%He=zeros(n,n);
%He=[1200*x(1)^2-400*x(2)+2, -400*x(1); 
 %        -400*x(1),                         200        ];
 

3）最速下降法

优点：最速下降法只需要求函数的梯度，并不需要求函数的Hessen矩阵，而且最速梯度法对于初始值的选取并不敏感，全局收敛。

缺点：迭代速度较慢。

代码：

function [x, val, k] = zuisuxiajiang(x0)
    %使用最速下降法求解最小值
    %x是迭代后的最小值点，val是最小值，k是迭代次数
    %x0是初始值，最速下降法对于初始值不敏感
    k = 0;
    epsilon = 1e-10;    %终止误差值
    %rho和sigma是线搜索法的两个参数
    rho = 0.55;
    sigma = 0.5;
    %当函数梯度的范数小于等于终止误差值时停止迭代
    while norm(g(x0)) > epsilon
        d = -g(x0);
        m = 0;
        while (f(x0 + rho^m*d) > (f(x0)+sigma*rho^m*g(x0)'*d))
            m = m + 1;
        end
        x0 = x0 + rho^m*d;
        k = k + 1;
    end
    x = x0;
    val = f(x0);
end
 
%%原函数
function y = f(x)
    y = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;
end
 
%%原函数的梯度,要弄成列向量的形式
function y = g(x)
    y(1,1) = 400*x(1)*(x(1)^2 - x(2)) + 2*(x(1) - 1);
    y(2,1) = -200*(x(1)^2 - x(2));
end
 
%%原函数的Hessen矩阵，注意矩阵每个元素的值
function y = G(x)
    y(1, 1) = 1200*x(1)^2 - 400*x(2) + 2;
    y(1, 2) = -400*x(1);
    y(2, 1) = -400*x(1);
    y(2, 2) = 200;
end

（4）信赖域方法
优点：信赖域方法是一种保证全局收敛的优化算法，并且可以直接确定位移，产生新的迭代点。
思想：（求出位移，并根据位移的相关信息调整信赖域半径）

算法：

代码：

function [xk,val,k]=trustm(x0)
%功能: 牛顿型信赖域方法求解无约束优化问题 min f(x)
%输入: x0是初始迭代点
%输出: xk是近似极小点, val是近似极小值, k是迭代次数
n=length(x0);  x=x0; dta=1;
eta1=0.15; eta2=0.75;  dtabar=2.0;
tau1=0.5; tau2=2.0; epsilon=1e-6;
k=0;  Bk=Hess(x);  %Bk=eye(n);  
while(k<150)
    gk=gfun(x);  %%%% 梯度值
    if(norm(gk)<epsilon) %终止条件
        break;
    end
    [d,val,lam,ik]=trustq(gk,Bk,dta);%%%%%信頼域子问题
    deltaq=-qk(x,d); %%%%预测下降
    deltaf=fun(x)-fun(x+d);%%%实际下降
    rk=deltaf/deltaq;%%%%%二者比值
    if(rk<=eta1)
        dta=tau1*dta;%%%%缩小信頼域半径，只有这种情况需要重新计算子问题
    else if (rk>=eta2&norm(d)==dta)
            dta=min(tau2*dta,dtabar); %%%%放大
        else
            dta=dta;  %%%维持不动
        end
    end
    if(rk>eta1)
        x0=x;     x=x+d;    %%%产生新的迭代点
       % sk=x-x0;  yk=gfun(x)-gfun(x0);  
        %vk=sqrt(yk'*Bk*yk)*(sk/(sk'*yk)-Bk*yk/(yk'*Bk*yk));
        %Bk=Bk-Bk*yk*yk'*Bk/(yk'*Bk*yk)+sk*sk'/(sk'*yk)+vk*vk'
        %pause
        Bk=Hess(x);
    end
    k=k+1;
end
xk=x;
val=fun(xk);
%%% 目标函数  %%%%%%%%%%%%%%%
function f=fun(x)
 f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;
 %%% 子问题目标函数 %%%%%%%%%%%%%
function qd=qk(x,d)
gk=gfun(x);  Bk=Hess(x);
qd=gk'*d+0.5*d'*Bk*d;
%%% 目标函数的梯度 %%%%%%%%%%%%%%
function gf=gfun(x)
gf=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1), -200*(x(1)^2-x(2))]';
%%% 目标函数的Hesse阵 %%%%%%%%%%%%%%
function He=Hess(x)
n=length(x);
He=zeros(n,n);
He=[1200*x(1)^2-400*x(2)+2, -400*x(1); 
         -400*x(1),                         200        ];