Progressive Hedging
This tutorial was generated using Literate.jl. Download the source as a .jl file.
This tutorial demonstrates the Progressive Hedging (PH) algorithm, a decomposition method for stochastic programming that solves scenario subproblems iteratively. It may be helpful to read Two-stage stochastic programs first.
Learning intentions:
- Understand how Progressive Hedging decomposes a stochastic program into per-scenario subproblems linked by a quadratic penalty toward the current consensus first-stage decision
- Use
Parametervariables andset_parameter_valueto update penalty coefficients and dual prices between iterations without rebuilding the model - Extend the basic algorithm with an adaptive penalty rule that balances primal and dual residuals to speed convergence
Required packages
This tutorial uses the following packages:
using JuMPimport Distributionsimport Ipoptimport PrintfBackground
Progressive Hedging (PH) is a popular decomposition algorithm for stochastic programming. It decomposes a stochastic problem into scenario subproblems that are solved iteratively, with penalty terms driving solutions toward consensus.
In Progressive Hedging, each scenario subproblem includes a quadratic penalty term:
\[\min\limits_{x_s}: f_s(x_s) + \frac{\rho}{2} ||x_s - \bar{x}||^2 + w_s^\top x_s\]
where:
- $x_s$ is the primal variable in scenario $s$
- $f_s(x)$ is the original scenario objective
- $\rho$ is the penalty parameter
- $\bar{x}$ is the current consensus (average) solution
- $w_s$ is the dual price (Lagrangian multiplier) in scenario $s$
Progressive Hedging is an iterative algorithm. In each iteration, it solves all the penalized scenario subproblems, then it applies two updates:
- $\bar{x} = \mathbb{E}_s[x_s]$
- $w_s = w_s + \rho (x_s - \bar{x})$
The algorithm terminates if $|\bar{x} - x_s| \le \varepsilon$ for all scenarios (the primal residual), and $\bar{x}$ has not changed by much between iterations (the dual residual).
$\rho$ can be optionally updated between iterations. How to do so is an open question. There is a large literature on different update strategies.
In this tutorial we use parameters for $\rho$, $w$, and $\bar{x}$ to efficiently modify each scenario's subproblem between PH iterations.
Building a single scenario
The building block of Progressive Hedging is a separate JuMP model for each scenario. Here's an example, using the problem from Two-stage stochastic programs:
function build_subproblem(; demand::Float64) model = Model(Ipopt.Optimizer) set_silent(model) @variable(model, x >= 0) @variable(model, 0 <= y <= demand) @constraint(model, y <= x) @variable(model, ρ in Parameter(1)) @variable(model, x̄ in Parameter(0)) @variable(model, w in Parameter(0)) @expression(model, f_s, 2 * x - 5 * y + 0.1 * (x - y)) @objective(model, Min, f_s + ρ / 2 * (x - x̄)^2 + w * x) return modelendbuild_subproblem (generic function with 1 method)Using the build_subproblem function, we can create one JuMP model for each scenario:
N = 10demands = rand(Distributions.TriangularDist(150.0, 250.0, 200.0), N);subproblems = map(demands) do demand return (; model = build_subproblem(; demand), probability = 1 / N)end;The Progressive Hedging loop
We're almost ready for our optimization loop, but first, here's a helpful function for logging:
function print_iteration(iter, args...) if mod(iter, 10) == 0 f(x) = Printf.@sprintf("%15.4e", x) println(lpad(iter, 9), " ", join(f.(args), " ")) end returnendprint_iteration (generic function with 1 method)Now we can implement our algorithm:
function solve_progressive_hedging( subproblems; iteration_limit::Int = 400, atol::Float64 = 1e-4, ρ::Float64 = 1.0,) x̄_old, x̄ = 0.0, 0.0 x, w = zeros(length(subproblems)), zeros(length(subproblems)) println("iteration primal_residual dual_residual") # For each iteration... for iter in 1:iteration_limit # For each subproblem... for (i, data) in enumerate(subproblems) # Update the parameters set_parameter_value(data.model[:ρ], ρ) set_parameter_value(data.model[:x̄], x̄) set_parameter_value(data.model[:w], w[i]) # Solve the subproblem optimize!(data.model) assert_is_solved_and_feasible(data.model) # Store the primal solution x[i] = value(data.model[:x]) end # Compute the consensus solution for the first-stage variables x̄ = sum(s.probability * x_s for (s, x_s) in zip(subproblems, x)) # Compute the primal and dual residuals primal_residual = maximum(abs, x_s - x̄ for x_s in x) dual_residual = ρ * abs(x̄ - x̄_old) print_iteration(iter, primal_residual, dual_residual) # Check for convergence if primal_residual < atol && dual_residual < atol break end # Update x̄_old = x̄ w .+= ρ .* (x .- x̄) end return x̄endx̄ = solve_progressive_hedging(subproblems);iteration primal_residual dual_residual
10 1.2434e-13 3.0000e+00
20 2.4158e-13 3.0000e+00
30 5.2580e-13 3.0000e+00
40 2.3590e-12 3.0000e+00
50 1.8247e-11 3.0000e+00
60 6.5765e-10 2.4900e+00
70 1.5246e+00 1.8106e+00
80 9.1392e-02 6.0000e-02
90 1.2808e-08 6.0000e-02
100 8.5467e-02 2.8523e-02
110 4.6490e-02 1.8070e-02
120 2.4945e-02 1.1333e-02
130 1.3165e-02 7.0452e-03
140 6.8049e-03 4.3447e-03
150 3.4229e-03 2.6597e-03
160 1.6581e-03 1.6171e-03
170 7.5900e-04 9.7687e-04
180 3.1565e-04 5.8650e-04
190 1.0705e-04 3.5004e-04
200 1.6006e-05 2.0770e-04
210 1.8480e-05 1.2253e-04The consensus first-stage decision is:
x̄206.71538424630526Progressive Hedging with an adaptive penalty parameter
You can also make the penalty parameter $\rho$ adaptive. How to do so is an open question. There is a large literature on different update strategies. One approach is to increase $\rho$ if the primal residual is much larger than the dual residual, and to decrease $\rho$ if the dual residual is much larger than the primal residual.
function solve_adaptive_progressive_hedging( subproblems; iteration_limit::Int = 400, atol::Float64 = 1e-4, ρ::Float64 = 1.0, τ::Float64 = 1.3, μ::Float64 = 15.0,) x̄_old, x̄ = 0.0, 0.0 x, w = zeros(length(subproblems)), zeros(length(subproblems)) println("iteration primal_residual dual_residual") for iter in 1:iteration_limit for (i, data) in enumerate(subproblems) set_parameter_value(data.model[:ρ], ρ) set_parameter_value(data.model[:x̄], x̄) set_parameter_value(data.model[:w], w[i]) optimize!(data.model) assert_is_solved_and_feasible(data.model) x[i] = value(data.model[:x]) end x̄ = sum(s.probability * x_s for (s, x_s) in zip(subproblems, x)) primal_residual = maximum(abs, x_s - x̄ for x_s in x) dual_residual = ρ * abs(x̄ - x̄_old) print_iteration(iter, primal_residual, dual_residual) if primal_residual < atol && dual_residual < atol break end w .+= ρ .* (x .- x̄) x̄_old = x̄ # Adaptive ρ update if primal_residual > μ * dual_residual ρ *= τ elseif dual_residual > μ * primal_residual ρ /= τ end end return x̄endx̄ = solve_adaptive_progressive_hedging(subproblems);iteration primal_residual dual_residual
10 1.5011e-10 3.0000e+00
20 3.8219e-10 6.0000e-02
30 7.3621e-01 4.5341e-02
40 3.5820e-01 1.6047e-02
50 2.0691e-01 1.6670e-02
60 1.1946e-01 1.5301e-02
70 6.8947e-02 1.3164e-02
80 3.8299e-02 1.2519e-02
90 2.1403e-02 1.0986e-02
100 1.5644e-02 9.1913e-03
110 8.6489e-03 5.9248e-03
120 4.7301e-03 3.7721e-03
130 2.5063e-03 3.1206e-03
140 1.7403e-03 2.5188e-03
150 1.2466e-03 1.9808e-03
160 6.4078e-04 1.2194e-03
170 3.1992e-04 7.4531e-04
180 1.5018e-04 5.8933e-04
190 9.2911e-05 4.5994e-04
200 5.4523e-05 2.7475e-04
210 1.9178e-05 1.6410e-04The consensus first-stage decision is:
x̄206.71550019405603Try tuning the values of τ and μ. Can you get the algorithm to converge in fewer iterations?