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 updates 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.4869e-13 3.0000e+00
30 5.2580e-13 3.0000e+00
40 2.3448e-12 3.0000e+00
50 1.8645e-11 3.0000e+00
60 3.0295e-10 2.4900e+00
70 3.5567e-10 1.4700e+00
80 1.3594e-10 4.5000e-01
90 1.0039e-10 4.5000e-01
100 3.4697e-01 4.1145e-01
110 5.8841e-10 6.0000e-02
120 6.6430e-08 6.0000e-02
130 3.8446e-02 3.5288e-02
140 1.7896e-02 2.1351e-02
150 7.6684e-03 1.2838e-02
160 2.7903e-03 7.6736e-03
170 6.1194e-04 4.5598e-03
180 2.5229e-04 2.6939e-03
190 5.1047e-04 1.5824e-03
200 5.1314e-04 9.2406e-04
210 4.2628e-04 5.3642e-04
220 3.2306e-04 3.0947e-04
230 2.3179e-04 1.7739e-04
240 1.6031e-04 1.0099e-04The consensus first-stage decision is:
x̄212.96199680106278Progressive 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 updates 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.4859e-10 3.0000e+00
20 1.0441e+00 6.0000e-02
30 1.4831e+00 4.2482e-02
40 7.6746e-01 6.4743e-02
50 3.9929e-01 4.6819e-02
60 2.5251e-01 6.1416e-02
70 1.6810e-01 4.2108e-02
80 1.1027e-01 3.7627e-02
90 7.1751e-02 2.5459e-02
100 5.8669e-02 1.7451e-02
110 3.6128e-02 1.2071e-02
120 2.1899e-02 8.4360e-03
130 1.3410e-02 4.2488e-03
140 1.0677e-02 2.4693e-03
150 8.3335e-03 1.5271e-03
160 4.9783e-03 8.9629e-04
170 2.9734e-03 3.1547e-04
180 1.7656e-03 5.8957e-05
190 1.3543e-03 5.9853e-05
200 1.0346e-03 9.3503e-05
210 6.0475e-04 9.9162e-05
220 3.5150e-04 8.4215e-05
230 2.6324e-04 8.5537e-05
240 1.9818e-04 7.4850e-05
250 1.1317e-04 5.2520e-05The consensus first-stage decision is:
x̄212.9620043913626Try tuning the values of τ and μ. Can you get the algorithm to converge in fewer iterations?