信赖域方法求解最优化问题代码MATLAB

信赖域方法是一类求解无约束最优化问题的优化算法。下面是一个MATLAB代码示例，用于实现信赖域方法：

function [x, fval, exitflag] = trustregion(fun, x0, options)
% TRUSTREGION Trust region method for unconstrained optimization
%   [X, FVAL, EXITFLAG] = TRUSTREGION(FUN, X0, OPTIONS) finds the minimum
%   point X of the function FUN starting from the point X0 using trust
%   region method. FUN is a function handle to the objective function that
%   accepts a vector input X and returns a scalar objective function value
%   F. X0 is the initial point and OPTIONS is a structure created using the
%   OPTIMSET function. The field names in the OPTIONS structure correspond
%   to parameter names for the TRUSTREGION function. Default values for
%   parameters not specified in OPTIONS are:
%       OPTIONS.tol = 1e-6;
%       OPTIONS.maxIter = 1000;
%       OPTIONS.maxFunEvals = Inf;
%       OPTIONS.radius = 1;
%       OPTIONS.eta = 0.1;
%       OPTIONS.verbose = false;
%
%   X is the solution found by TRUSTREGION, FVAL is the objective function
%   value at X, and EXITFLAG is an integer identifying the reason for
%   termination:
%       1 - Function converged to a solution X satisfying the tolerance
%           criterion.
%       0 - Maximum number of iterations or function evaluations reached.
%
%   Example:
%       fun = @(x) (x(1)-1)^2 + (x(2)-2.5)^2;
%       x0 = [0 0];
%       options = optimset('tol', 1e-8, 'maxIter', 500);
%       [x, fval, exitflag] = trustregion(fun, x0, options);
%
%   References:
%       [1] Nocedal, J., & Wright, S. J. (2006). Numerical optimization
%           (2nd ed.). Springer.

% Set default options if not specified by user
defaultOptions.tol = 1e-6;
defaultOptions.maxIter = 1000;
defaultOptions.maxFunEvals = Inf;
defaultOptions.radius = 1;
defaultOptions.eta = 0.1;
defaultOptions.verbose = false;
if nargin < 3
    options = defaultOptions;
else
    % Fill in missing values of options with defaults
    optionNames = fieldnames(defaultOptions);
    for i = 1:length(optionNames)
        if ~isfield(options, optionNames{i})
            options.(optionNames{i}) = defaultOptions.(optionNames{i});
        end
    end
end

% Initialize variables
x = x0;
fval = fun(x);
grad = gradient(fun, x);
Hess = hessian(fun, x);
iter = 0;
funEvals = 1;

% Main loop
while true
    % Compute the Cauchy point
    [pk, mk] = cauchy(fun, x, grad, Hess, options.radius);
    
    % Compute the actual reduction
    actualReduction = fval - fun(x + pk);
    
    % Compute the predicted reduction
    predictedReduction = -mk + fun(x);
    
    % Compute the ratio of actual to predicted reduction
    rho = actualReduction / predictedReduction;
    
    % Update the trust region radius
    if rho < 0.25
        options.radius = options.radius / 4;
    elseif rho > 0.75 && abs(norm(pk) - options.radius) < options.tol
        options.radius = min(2 * options.radius, options.maxRadius);
    end
    
    % Update x and fval if the predicted reduction is good
    if rho > options.eta
        x = x + pk;
        fval = fun(x);
        grad = gradient(fun, x);
        Hess = hessian(fun, x);
        funEvals = funEvals + 1;
    end
    
    % Check termination criteria
    iter = iter + 1;
    if norm(grad) < options.tol || iter >= options.maxIter || ...
            funEvals >= options.maxFunEvals
        break;
    end
    
    % Output iteration information if verbose option is true
    if options.verbose
        fprintf('Iteration %d: fval = %g, radius = %g, norm(grad) = %g\n', ...
            iter, fval, options.radius, norm(grad));
    end
end

% Set exit flag based on termination reason
if norm(grad) < options.tol
    exitflag = 1;
else
    exitflag = 0;
end
end

function [pk, mk] = cauchy(fun, x, grad, Hess, radius)
% Compute the Cauchy point for trust region method
g = grad(:);
H = Hess(:);
alpha = norm(g)^3 / (radius * g' * H * g);
pk = -alpha * min(radius, norm(g)) * g;
mk = fun(x) + g' * pk + 0.5 * pk' * Hess * pk;
end

function grad = gradient(fun, x)
% Compute the gradient of the objective function using central difference
n = length(x);
grad = zeros(n, 1);
h = 1e-6;
for i = 1:n
    xplus = x;
    xplus(i) = x(i) + h;
    xminus = x;
    xminus(i) = x(i) - h;
    grad(i) = (fun(xplus) - fun(xminus)) / (2 * h);
end
end

function Hess = hessian(fun, x)
% Compute the Hessian of the objective function using central difference
n = length(x);
Hess = zeros(n, n);
h = 1e-6;
for i = 1:n
    xplus = x;
    xplus(i) = x(i) + h;
    xminus = x;
    xminus(i) = x(i) - h;
    gradplus = gradient(fun, xplus);
    gradminus = gradient(fun, xminus);
    Hess(:, i) = (gradplus - gradminus) / (2 * h);
end
end

上述代码实现了信赖域方法的主循环，以及计算梯度、海森矩阵和 Cauchy 点的函数。要使用此代码，您需要提供一个函数句柄 `fun`，该函数计算目标函数值，一个初始点 `x0`，以及一个选项结构 `options`，该结构包含有关算法的各种参数设置。请注意，此代码使用中心差分来估计梯度和海森矩阵，因此可能不适用于计算资源受限的问题。