数值计算器fminunc处理约束问题的惩罚/障碍函数法

 

在数学领域约束优化中，障碍函数是一个连续函数，其中点的值随着点到达优化问题的可行区域的边界而增加到无穷大。这些函数用于通过更容易处理的目标函数中的惩罚项来代替不等式约束。

两种最常见的屏障功能类型是反向屏障功能和对数屏障功能。 恢复对数屏障功能是由于它们与原始双重内点方法的联系。

罚函数法又称乘子法，是指将有约束最优化问题转化为求解无约束最优化问题：其中M为足够大的正数， 起惩罚作用， 称之为罚因子，F(x, M )称为罚函数。内部罚函数法也称为障碍罚函数法。

数值求解器迭代次数不够往往求得的解不是实际的解，慎用。

function [x_opt, fval] = Barrier_Function(f, g, x0, mu, tol, iter)
    % f: 目标函数句柄，例如 f = @(x) x(1) * x(2);
    % g: 约束函数句柄，例如 g = @(x) -2 * x(1) + x(2) + 3;
    % x0: 初始点
    % mu: 障碍函数的系数
    % tol: 容忍度
    % iter: 最大迭代次数

    % 定义障碍函数
    barrier_func = @(x) f(x) - mu * (1 / g(x));

    % 初始化
    options = optimoptions('fminunc', 'Algorithm', 'quasi-newton', 'Display', 'iter', 'TolFun', tol);
    x = x0; % 初始点
    k = 0; % 迭代计数器

    while k < iter
        % 使用 fminunc 求解无约束问题
        [x, fval] = fminunc(barrier_func, x, options);

        % 输出当前迭代信息
        fprintf('迭代次数: %d, 当前解: [%f, %f], 目标函数值: %f, 约束函数值: %f\n', ...
            k, x(1), x(2), fval, g(x));

        % 检查是否满足约束条件
        if g(x) >= tol
            disp('找到满足约束条件的解');
            break;
        end

        % 增大障碍系数 mu，使得解远离约束边界
        mu = mu * 0.1;
        barrier_func = @(x) f(x) - mu * (1 / g(x));

        % 增加迭代次数
        k = k + 1;
    end

    % 如果达到最大迭代次数，输出警告
    if k == iter
        warning('达到最大迭代次数，可能未找到满足约束的解。');
    end

    % 输出最终结果
    x_opt = x;
end

function [x_opt, fval] = Penalty_Function_Numerical(f, G, x0, rho, tol, iter)
    % f: 目标函数句柄，例如 f = @(x) -x(1) * x(2);
    % G: 约束函数句柄的单元数组，例如 G = {@(x) x(1)^2 + x(2)^2 - 1, @(x) -x(1) - x(2)};
    % x0: 初始点
    % rho: 初始惩罚参数
    % tol: 容忍度
    % iter: 最大迭代次数
    
    % 定义惩罚函数的数值形式
    penalty_func = @(x) f(x) + rho * sum(cellfun(@(g) max(0, g(x))^2, G));

    % 初始化
    options = optimoptions('fminunc', 'Algorithm', 'quasi-newton', 'Display', 'iter');
    x = x0; % 初始点
    k = 0; % 迭代计数器

    while k < iter
        % 使用 fminunc 求解无约束问题
        [x, fval] = fminunc(penalty_func, x, options);

        % 计算当前点处的约束值
        constraint_vals = cellfun(@(g) g(x), G);

        % 输出当前迭代信息
        fprintf('迭代次数: %d, 当前解: [%.2f, %.2f], 目标函数值: %.2f, 约束值: [%.2f, %.2f]\n', ...
            k, x(1), x(2), fval, constraint_vals(1), constraint_vals(2));

        % 检查是否满足约束条件
        if all(constraint_vals <= tol)
            disp('找到满足约束条件的解');
            break;
        end

        % 增大惩罚参数 rho
        rho = rho * 10;
        penalty_func = @(x) f(x) + rho * sum(cellfun(@(g) max(0, g(x))^2, G));

        % 增加迭代次数
        k = k + 1;
    end

    % 如果达到最大迭代次数，输出警告
    if k == iter
        warning('达到最大迭代次数，可能未找到满足约束的解。');
    end

    % 将解和目标函数值保留两位有效数字
    x_opt = round(x, 2);  % 保留两位有效数字
    fval = round(fval, 2); % 保留两位有效数字

    % 输出最终结果
    fprintf('最优解: [%.2f, %.2f]\n', x_opt(1), x_opt(2));
    fprintf('最优目标函数值: %.2f\n', fval);
end


%{
f = @(x) -x(1) * x(2); % 目标函数
G = {@(x) x(1)^2 + x(2)^2 - 1, @(x) -x(1) - x(2)}; % 约束函数
x0 = [0.5; -0.5]; % 初始点
rho = 10; % 初始惩罚参数
tol = 1e-6; % 约束容忍度
iter = 10; % 最大迭代次数

[x_opt, fval] = Penalty_Function_Numerical(f, G, x0, rho, tol, iter);
%}