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 JuMP
import Distributions
import Ipopt
import 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 model
endbuild_subproblem (generic function with 1 method)Using the build_subproblem function, we can create one JuMP model for each scenario:
N = 10
demands = 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
return
endprint_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̄
end
x̄ = solve_progressive_hedging(subproblems);iteration primal_residual dual_residual
10 1.1369e-13 3.0000e+00
20 2.0606e-13 3.0000e+00
30 3.8369e-13 3.0000e+00
40 1.7053e-12 3.0000e+00
50 1.3642e-11 3.0000e+00
60 2.6369e-09 1.9800e+00
70 1.4259e+00 2.3455e-01
80 1.3795e-09 6.0000e-02
90 1.2303e-01 3.0765e-02
100 7.6570e-02 1.6277e-02
110 4.7271e-02 8.4409e-03
120 2.8966e-02 4.2643e-03
130 1.7628e-02 2.0784e-03
140 1.0658e-02 9.6058e-04
150 6.4046e-03 4.0648e-04
160 3.8256e-03 1.4372e-04
170 2.2718e-03 2.7516e-05
180 1.3413e-03 1.7711e-05
190 7.8736e-04 3.0450e-05
200 4.5950e-04 2.9682e-05
210 2.6656e-04 2.4336e-05The consensus first-stage decision is:
x̄191.6017963845735Progressive 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̄
end
x̄ = solve_adaptive_progressive_hedging(subproblems);iteration primal_residual dual_residual
10 1.0938e-10 3.0000e+00
20 6.2224e-01 1.2090e-02
30 1.5750e-01 7.1293e-03
40 9.3203e-02 4.1380e-03
50 5.5305e-02 6.2773e-04
60 4.2403e-02 9.5003e-04
70 2.4891e-02 1.6291e-03
80 1.4556e-02 1.7755e-03
90 8.4105e-03 2.1262e-03
100 6.3039e-03 2.1000e-03
110 4.7481e-03 1.8173e-03
120 2.7102e-03 1.2725e-03
130 1.5359e-03 8.6481e-04
140 8.5876e-04 7.5518e-04
150 6.2855e-04 6.2898e-04
160 4.5822e-04 3.9732e-04
170 2.5110e-04 2.5347e-04
180 1.3592e-04 1.5995e-04
190 9.3357e-05 1.3047e-04
200 6.6387e-05 1.0320e-04The consensus first-stage decision is:
x̄191.60183084576389Try tuning the values of τ and μ. Can you get the algorithm to converge in fewer iterations?