问题描述：

$min f(x)=x_1^2+x_2^2$

$s.t. x_1+2x_2-2=0$

PH算法执行过程：

Step 1：选定初始点$x_0$、初始乘子向量$\lambda_1$，初始罚因子$\sigma_1$、及其放大系数 $c>1$、控制误差 $\epsilon>0$与常数 $\theta\in(0,1)$，令$k=1$

Step 2：以$x_{k-1}$为初始点求解无约束问题$min M(x,\lambda_k,\sigma_k)=f(x)-\lambda_k^\mathrm{T}c(x)+\frac{\sigma_k}{2}c(x)^\mathrm{T}x(x)$得最优解$x_k$

Step 3：当 $\Vert c(x_k)\Vert<\epsilon$时，$x_k$为所求最优解，停。否则转Step 4

Step 4：当$\Vert c(x_k)\Vert/\Vert c(x_{k-1})\Vert\leq\theta$ 时,转 Step 5，否则令$\sigma_{k-1}=c\sigma_k$，转 Step 5.

Step 5：令$\lambda_{k+1}=\lambda_k-\sigma_kc(x_k),k=k+1$，转 Step 2

Matlab代码实现：

function x_optimal = ph_algorithm_no_toolbox()
    % Step 1: Initialization
    x = [0; 0];          % Initial point
    lambda = zeros(1, 1); % Initial Lagrange multiplier
    sigma = 1;           % Initial penalty parameter
    c = 2;               % Penalty amplification factor
    epsilon = 1e-6;      % Tolerance for stopping criterion
    theta = 0.1;         % Threshold for penalty adjustment
    k = 1;               % Iteration counter

    while true
        % Step 2: Solve unconstrained problem
        x_k = solve_unconstrained_no_toolbox(x, lambda, sigma);

        % Step 3: Check stopping criterion
        if norm(constraint_function(x_k)) < epsilon
            x_optimal = x_k; % Optimal solution found
            break;
        end

        % Step 4: Check for penalty adjustment
        if norm(constraint_function(x_k)) / norm(constraint_function(x)) <= theta
            % Skip penalty adjustment, proceed to Step 5
        else
            sigma = c * sigma; % Adjust penalty parameter
        end

        % Step 5: Update Lagrange multiplier
        lambda = lambda - sigma * constraint_function(x_k);

        % Prepare for the next iteration
        x = x_k;
        k = k + 1;
    end
end

function x_optimal = solve_unconstrained_no_toolbox(x, lambda, sigma)
    % Objective function
    objective_function = @(x) x(1)^2 + x(2)^2;

    % Lagrangian function with penalty term
    lagrangian_function = @(x) objective_function(x) - lambda * constraint_function(x) + (sigma / 2) * norm(constraint_function(x))^2;

    % Solve unconstrained problem using fminsearch (no toolbox required)
    x_optimal = fminsearch(lagrangian_function, x);
end

function c_x = constraint_function(x)
    % Equality constraint
    c_x = x(1) + 2 * x(2) - 2;
end

本文由mdnice多平台发布