优化算法：信赖域法（含matlab代码）

Williamlliw

于 2025-01-30 16:07:23 发布

阅读量324

点赞数 11

文章标签：算法 matlab 开发语言

本文链接：https://blog.csdn.net/m0_62150241/article/details/144184846

版权

1.简介

在求解优化问题时，信赖域法不要求目标函数的凸性。其基本思想是在迭代点附近（信赖域内）用二次型逼近目标函数，并求其极小值作为下一个迭代点。当迭代后函数下降不充分时，调整信赖域的大小。

2.算法介绍

以下内容引用自文献[1]。

3.matlab代码

3.1 信赖域法matlab函数

%f：目标函数
%H0：初始的海森矩阵
%x0：初始的迭代点 为列向量
%k：迭代次数
%X存储每次迭代的x，F为目标函数，G为每次的梯度，H为海森阵，HN为海森矩阵的逆

function [X, F, G, H] = trust_region(x0, f, my_gradient, my_hessian)
%% 变量初始化
max_step = 1e3;
eps = 1e-5; 
m = length(x0);
X = zeros(m, max_step+1); % x1、x2的迭代值
H = zeros(m, m, max_step+1); % hessian的迭代值
F = zeros(max_step+1, 1); % function的迭代值
G = zeros(m, max_step+1); % 梯度的迭代值
% 赋初始值
X(:,1) = x0; % X(x1, x2)初始化
F(1, 1) = f(x0);
% 信赖域问题的参数
exp_rou = 0.8; % 预期下降量
gamma = 0.7; % 信赖域半径调整参数
delta = 2; % 初始信赖域半径
     
%% 迭代
for k = 1: max_step
    G(:, k) = my_gradient(X(:, k)); % 求梯度
    H(:, :, k) = my_hessian(X(:, k)); % 求hessian阵

    % 求信赖域子问题
    A=[];  Aeq=[];  b=[];  beq=[]; lb=[]; ub=[];
    trust_f = @(s) sub_f(s, G(:, k), H(:, :, k));
    trust_con = @(s) con_sub_f(s, delta);
    [trust_step, trust_fval] = fmincon(trust_f, X(:, k), A,b,Aeq,beq,lb,ub,trust_con);
    rou = (f(X(:,k)) - f(X(:,k) + trust_step))/(trust_f([0;0]) - trust_fval);
    % x迭代
    if rou > exp_rou && isreal(rou)
        X(:,k+1) = X(:,k) + trust_step;
    else
        X(:,k+1) = X(:,k);
        % 调整信赖域半径
        delta = gamma*delta;
    end

    F(k) = f(X(:, k));
    % 梯度消失后结束程序
    if norm(G(:, k)) < eps
        break;
    end
end
% 输出
X = X(:, 1:k);
H = H(:, :, 1:k);
F = F(1: k);
G = G(:, 1:k);
end

% 信赖域子问题的约束函数
function [c, ceq] = con_sub_f(s, delta)
c =  s'*s - delta^2;
ceq = [];
end

% 信赖域子问题的目标函数
function out = sub_f(s, G, H)
if length(s(1, :)) > 1
    s = s';
end
out = G'*s + 1/2*s'*H*s;
end

3.2 一个简单的测试

f = @(x) x(1).^2 + 3*x(2).^2 + exp(x(1));
grd = @(x) [2*x(1) + exp(x(1)); 6*x(2)];
hes = @(x) [2 + exp(x(1)), 0;
            0, 6];
x0 = [3; 3];
[X, F, G, H] = trust_region(x0, f, grd, hes);

3.3 非凸算例

考虑如下的问题，其目标函数非凸，在最优解附近凸。该问题来自于北航最优化理论与方法I2024秋季的大作业。

计算该函数的函数值，梯度，hessian阵。

function [cost, gradient, my_hessian] = target_function(input_x1, input_x2, flag)
%% 初始化
cost = 0;
gradient = 0;
my_hessian = 0;
%% 符号函数运算
syms x1 x2
a = [0.3; 0.6; 0.2];
b = [5; 26; 3];
c = [40; 1; 10];
u = x1 - 0.8;
v = x2 - a'*[1; u.^2.*sqrt(1-u); -u];
alpha = b'*[-1; u.^2.*sqrt(1+u); u];
beta = c(1)*v.^2.*(1-c(2)*v)./(1+c(3)*u.^2);
f = alpha.*exp(-beta);
%% 数值运算
% 目标函数
if ~(flag == 2 || flag == 3)
    cost_fun = matlabFunction(f);
    cost = cost_fun(input_x1, input_x2);
end
% 梯度函数
if ~(flag == 1 || flag == 3)
    grad_f(1) = diff(f, x1);
    grad_f(2,1) = diff(f, x2);
    grad_f_fun = matlabFunction(grad_f);
    gradient = grad_f_fun(input_x1, input_x2);
end
% Hessian矩阵
if ~(flag == 1 || flag == 2)
    grad_f(1) = diff(f, x1);
    grad_f(2,1) = diff(f, x2);
    H = jacobian(grad_f, [x1;x2]);
    H_fun = matlabFunction(H);
    my_hessian = H_fun(input_x1, input_x2);
end
end

信赖域法求解该问题：

% 信赖域法求解
clear; clc; close all;
%% 信赖域法求解
% 第一个点
x0_1 = [0.8; 0.3];
[x_rec1, f_rec1, G_rec1, H_rec1] = trust_region(x0_1, @cal_cost, @cal_gradient, @cal_hessian);
% 第二个点
x0_2 = [1; 0.5];
[x_rec2, f_rec2, G_rec2, H_rec2] = trust_region(x0_2, @cal_cost, @cal_gradient, @cal_hessian);
%% 画图
figure(1);
plot_equal_height;
hold on;
plot(x_rec1(1, :), x_rec1(2, :), 'o', 'LineWidth', 1.5, 'Color', [0.7, 0, 0]);
plot(x_rec1(1, :), x_rec1(2, :), 'LineWidth', 1.5, 'Color', [0.7, 0, 0]);
title('信赖域法迭代点，x_0 = (0.8, 0.3)');

figure(2);
plot_equal_height;
hold on;
plot(x_rec2(1, :), x_rec2(2, :), 'o', 'LineWidth', 1.5, 'Color', [0.7, 0, 0]);
plot(x_rec2(1, :), x_rec2(2, :), 'LineWidth', 1.5, 'Color', [0.7, 0, 0]);
title('信赖域法迭代点，x_0 = (1, 0.5)');

%% 计算目标函数
function out = cal_cost(x)
[out, ~, ~] = target_function(x(1), x(2), 1);
end

% 计算目标函数梯度
function out = cal_gradient(x)
[~, out, ~] = target_function(x(1), x(2), 2);
end

% 计算目标函数hessian阵
function out = cal_hessian(x)
[~, ~, out] = target_function(x(1), x(2), 3);
end

求解结果如下：