无约束优化问题的信赖域求解方法

最新推荐文章于 2026-03-31 09:45:00 发布

原创最新推荐文章于 2026-03-31 09:45:00 发布 · 1k 阅读

28 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

这里写目录标题

一、信赖域方法基本原理
- 1.模型
- 2.步骤
二、matlab仿真
- 1.程序
- 2.仿真结果
三、总结

一、信赖域方法基本原理

1.模型

信赖域方法的基本思想是：首先给定一个所谓的信赖域半径作为位移长度的上界，并以当前迭代点为中心，以此上界为半径确定一个信赖域，通过求解这个区域的信赖域子问题的最优点来确定候选位移。若候选位移能使目标函数值有充分的下降量，则接受该候选位移作为新的位移，并保持或扩大信赖域半径，继续新的迭代。

2.步骤

1）初始化
设置初始迭代点 $x_0$ 和初始信赖域半径，通常取Δ =1或基于问题尺度估计。同时设定收敛阈值 $\epsilon$
2）模型构建
在当前点处构建二次近似模型：
$m_k(p) = f(x_k) + \delta f(x_k)^{T} p+1/2p^{T}B_{k}p$
其中为 $B_{k}$ 为Hessian矩阵或其近似。
3）子问题求解
求解带约束的二次优化子问题：
$min m_k(p) s.t. ||p||<= Δ_k$
4）信赖度评估
$Δf_k$ 为f在第k步的实际下降量， $Δf_k=f(k)-f(x_k+d_k)$ ,
$Δq_k$ 为对应的预测下降量： $Δq_k=q_k(0)-q_k(d_k)$ ,
$r_k = Δf_k/Δq_k$
当 $r_k$ 接近1的时候， $x_{k+1} := x_k+d_k$ 可以作为新的迭代点。
5）迭代更新
k = k+1。

二、matlab仿真

1.程序


%-----function:Unconstrained Nonlinear Optimization------
%-----Remark:Use TrustRegionMethod-----------------------
%-----Data:2025.12.19------------------------------------
%-----Author:Clemence------------------------------------
%-----Reference:<Optimization Calculation Method and Its Matlab Program Implementation>


clc;
clear all;
close all;

%-------------Example-----------------------
x0 = [1,-1]';
epsilon = 1e-6;

[k,x,val]=TrustRegionMethod(x0,epsilon)


%--------------TrustRegion---------------
%   k:Iteration Count
%   x:Optimal Solution
%   val:Objective Function Value
function [k,x,val] = TrustRegionMethod(x0,epsilon)
    
    n = length(x0);
    eta1 = 0.1;
    eta2 = 0.75;
    tau1 = 0.5;
    tau2 = 2.0;
    delta  =1;
    dtabar = 2.0;
    

    x = x0;
   
    bk = hess(x);
    k = 0;
    
    while(k<50)
         fk = fun(x);
         gk = gfun(x);
    if(norm(gk)<epsilon)
        break;
    end
    %Call the Subprogram trustq
    [d,val,lam,i]=trustq(fk,gk,bk,delta);
    deltaq = fk - val;
    deltaf = fun(x) -fun(x+d);
    rk = deltaf/deltaq;
    if(rk<=eta1)
        delta =tau1*delta;
    else if(rk>=eta2 & norm(d)==delta)
            delta = min(tau2*delta,dtabar);
        else
            delta = delta;
        end
    end
    if(rk>eta1)
        x = x +d;
        bk = hess(x);
    end
    k = k+1;
    end
    val = fun(x);
end

%---Function：Solving the TrustRegion Subproblem-----
%   fk:the Objetctive function value at xk,
%   gk:the Gradient at xk
%   bk:the Approximate Hessian Matrix at Iteration k,
%   delta:the Current TrustRegion Radius
%---Output：-----------------------------------------
%   d:the Optimal Point
%   val:the Optimal Value of the Subproblem
%   lam:the Lagrange Multiplier Value
%   i:the Iteration Count
function [d,val,lam,i]=trustq(fk,gk,bk,deltak)
    n = length(gk);
    beta = 0.6;
    sigma = 0.05;
    mu0 = 0.05;
    lam0 = 0.05;
    gamma = 0.05;
    d0 = ones(n,1);
    z0 = [mu0,lam0,d0']';
    zbar = [mu0,zeros(1,n+1)]';
    i = 0;
    z = z0;
    mu = mu0;
    lam = lam0;
    d = d0;
    while(i<150)
        h = dah(mu,lam,d,gk,bk,deltak);
        if(norm(h)<=1e-8)
            break;
        end
        j = jacobih(mu,lam,d,bk,deltak);
        b = psi(mu,lam,d,gk,bk,deltak,gamma)*zbar-h;
        dz = j\b;
        dmu = dz(1);
        dlam = dz(2);
        dd = dz(3:n+2);
        m=0;
        mi=0;
        while(m<20)
            t1=beta^m;
            hnew = dah(mu+t1*dmu,lam+t1*dlam,d+t1*dd,gk,bk,deltak);
            if(norm(hnew)<=(1-sigma*(1-gamma*mu0)*beta^m)*norm(h))
                mi=m;
                break;
            end
            m=m+1;
        end
        alpha = beta^mi;
        mu = mu+alpha*dmu;
        lam = lam + alpha*dlam;
        d = d+alpha*dd;
        i = i+1;
    end
    val = fk+gk'*d+0.5*d'*bk*d;
end

%--------Function phi--------------
function p=phi(mu,a,b)
    p=a+b-sqrt((a-b)^2+4*mu^2);
end

%--------Matrix H--------------
function h = dah(mu,lam,d,gk,bk,deltak)
    n = length(d);
    h = zeros(n+2,1);
    h(1) = mu;
    h(2) = phi(mu,lam,deltak^2-norm(d)^2);
    h(3:n+2)=(bk+lam*eye(n))*d+gk;
end

%--------Matrix H'--------------
function j=jacobih(mu,lam,d,bk,deltak)
    n=length(d);
    j=zeros(n+2,n+2);
    t2 = sqrt((lam+norm(d)^2-deltak^2)^2+4*mu^2);
    pmu=-4*mu/t2;
    thetak=(lam+norm(d)^2-deltak^2)/t2;
    j = [1,0,zeros(1,n);
        pmu,1-thetak,-2*(1+thetak)*d';
        zeros(n,1),d,bk+lam*eye(n)];
end

function si = psi(mu,lam,d,gk,bk,deltak,gamma)
    h=dah(mu,lam,d,gk,bk,deltak);
    si = gamma*norm(h)*min(1,norm(h));
end


%-----Original Function------
function [f]=fun(x)
    f = 100*(x(1)^2-x(2))^2+(x(1)-1)^2;
end

%------First-order Partial Derivative Function-----
function [gf]=gfun(x)
    gf = [400*x(1)*(x(1)^2-x(2))+2*(x(1)-1);-200*(x(1)^2-x(2))];
end

%------Second-order Partial Derivative Function-----
function [h]=hess(x)
    h = [1200*x(1)^2-400*x(2)+2 -400*x(1);-400*x(1) 200];
end