ChainRules integration demo: Relaxed Unit Commitment

In this example, we will demonstrate the integration of DiffOpt with ChainRulesCore.jl, the library allowing the definition of derivatives for functions that can then be used by automatic differentiation systems.

using DiffOpt, Plots, JuMP
using LinearAlgebra
import Clp
using ChainRulesCore

Unit commitment problem

We will consider a unit commitment problem, finding the cost-minimizing activation of generation units in a power network over multiple time periods. The considered constraints include:

  • Demand satisfaction of several loads
  • Ramping constraints
  • Generation limits.

The decisions are:

  • $u_{it} \in \{0,1\}$: activation of the $i$-th unit at time $t$
  • $p_{it}$: power output of the $i$-th unit at time $t$.

DiffOpt handles convex optimization problems only, we therefore relax the domain of the $u_{it}$ variables to $\left[0,1\right]$.

Primal UC problem

ChainRules defines the differentiation of functions. The actual function that is differentiated in the context of DiffOpt is the solution map taking in input the problem parameters and returning the solution.

function unit_commitment(
        load1_demand, load2_demand, gen_costs, noload_costs;
        model = Model(Clp.Optimizer), silent=false)
    MOI.set(model, MOI.Silent(), silent)

    ## Problem data
    unit_codes = [1, 2] # Generator identifiers
    load_names = ["Load1", "Load2"] # Load identifiers
    n_periods = 4 # Number of time periods
    Pmin = Dict(1 => fill(0.5, n_periods), 2 => fill(0.5, n_periods)) # Minimum power output (pu)
    Pmax = Dict(1 => fill(3.0, n_periods), 2 => fill(3.0, n_periods)) # Maximum power output (pu)
    RR = Dict(1 => 0.25, 2 => 0.25) # Ramp rates (pu/min)
    P0 = Dict(1 => 0.0, 2 => 0.0) # Initial power output (pu)
    D = Dict("Load1" => load1_demand, "Load2" => load2_demand) # Demand (pu)
    Cp = Dict(1 => gen_costs[1], 2 => gen_costs[2]) # Generation cost coefficient ($/pu)
    Cnl = Dict(1 => noload_costs[1], 2 => noload_costs[2]) # No-load cost ($)

    ## Variables
    # Note: u represents the activation of generation units.
    # Would be binary in the typical UC problem, relaxed here to u ∈ [0,1]
    # for a linear relaxation.
    @variable(model, 0 <= u[g in unit_codes, t in 1:n_periods] <= 1) # Commitment
    @variable(model, p[g in unit_codes, t in 1:n_periods] >= 0) # Power output (pu)

    ## Constraints

    # Energy balance
    @constraint(
        model,
        energy_balance_cons[t in 1:n_periods],
        sum(p[g, t] for g in unit_codes) == sum(D[l][t] for l in load_names),
    )

    # Generation limits
    @constraint(model, [g in unit_codes, t in 1:n_periods], Pmin[g][t] * u[g, t] <= p[g, t])
    @constraint(model, [g in unit_codes, t in 1:n_periods], p[g, t] <= Pmax[g][t] * u[g, t])

    # Ramp rates
    @constraint(model, [g in unit_codes, t in 2:n_periods], p[g, t] - p[g, t - 1] <= 60 * RR[g])
    @constraint(model, [g in unit_codes], p[g, 1] - P0[g] <= 60 * RR[g])
    @constraint(model, [g in unit_codes, t in 2:n_periods], p[g, t - 1] - p[g, t] <= 60 * RR[g])
    @constraint(model, [g in unit_codes], P0[g] - p[g, 1] <= 60 * RR[g])

    # Objective
    @objective(
        model,
        Min,
        sum((Cp[g] * p[g, t]) + (Cnl[g] * u[g, t]) for g in unit_codes, t in 1:n_periods),
    )

    optimize!(model)
    # asserting finite optimal value
    @assert termination_status(model) == MOI.OPTIMAL
    # converting to dense matrix
    return JuMP.value.(p.data)
end
unit_commitment (generic function with 1 method)
m = Model(Clp.Optimizer)
@show unit_commitment(
    [1.0, 1.2, 1.4, 1.6], [1.0, 1.2, 1.4, 1.6],
    [1000.0, 1500.0], [500.0, 1000.0],
    model=m, silent=true
)
println(m)
unit_commitment([1.0, 1.2, 1.4, 1.6], [1.0, 1.2, 1.4, 1.6], [1000.0, 1500.0], [500.0, 1000.0], model = m, silent = true) = [2.0 2.4 2.8 3.0; 0.0 0.0 0.0 0.19999999999999996]
Min 1000 p[1,1] + 500 u[1,1] + 1000 p[1,2] + 500 u[1,2] + 1000 p[1,3] + 500 u[1,3] + 1000 p[1,4] + 500 u[1,4] + 1500 p[2,1] + 1000 u[2,1] + 1500 p[2,2] + 1000 u[2,2] + 1500 p[2,3] + 1000 u[2,3] + 1500 p[2,4] + 1000 u[2,4]
Subject to
 energy_balance_cons[1] : p[1,1] + p[2,1] = 2.0
 energy_balance_cons[2] : p[1,2] + p[2,2] = 2.4
 energy_balance_cons[3] : p[1,3] + p[2,3] = 2.8
 energy_balance_cons[4] : p[1,4] + p[2,4] = 3.2
 0.5 u[1,1] - p[1,1] ≤ 0.0
 0.5 u[2,1] - p[2,1] ≤ 0.0
 0.5 u[1,2] - p[1,2] ≤ 0.0
 0.5 u[2,2] - p[2,2] ≤ 0.0
 0.5 u[1,3] - p[1,3] ≤ 0.0
 0.5 u[2,3] - p[2,3] ≤ 0.0
 0.5 u[1,4] - p[1,4] ≤ 0.0
 0.5 u[2,4] - p[2,4] ≤ 0.0
 -3 u[1,1] + p[1,1] ≤ 0.0
 -3 u[2,1] + p[2,1] ≤ 0.0
 -3 u[1,2] + p[1,2] ≤ 0.0
 -3 u[2,2] + p[2,2] ≤ 0.0
 -3 u[1,3] + p[1,3] ≤ 0.0
 -3 u[2,3] + p[2,3] ≤ 0.0
 -3 u[1,4] + p[1,4] ≤ 0.0
 -3 u[2,4] + p[2,4] ≤ 0.0
 -p[1,1] + p[1,2] ≤ 15.0
 -p[2,1] + p[2,2] ≤ 15.0
 -p[1,2] + p[1,3] ≤ 15.0
 -p[2,2] + p[2,3] ≤ 15.0
 -p[1,3] + p[1,4] ≤ 15.0
 -p[2,3] + p[2,4] ≤ 15.0
 p[1,1] ≤ 15.0
 p[2,1] ≤ 15.0
 p[1,1] - p[1,2] ≤ 15.0
 p[2,1] - p[2,2] ≤ 15.0
 p[1,2] - p[1,3] ≤ 15.0
 p[2,2] - p[2,3] ≤ 15.0
 p[1,3] - p[1,4] ≤ 15.0
 p[2,3] - p[2,4] ≤ 15.0
 -p[1,1] ≤ 15.0
 -p[2,1] ≤ 15.0
 u[1,1] ≥ 0.0
 u[2,1] ≥ 0.0
 u[1,2] ≥ 0.0
 u[2,2] ≥ 0.0
 u[1,3] ≥ 0.0
 u[2,3] ≥ 0.0
 u[1,4] ≥ 0.0
 u[2,4] ≥ 0.0
 p[1,1] ≥ 0.0
 p[2,1] ≥ 0.0
 p[1,2] ≥ 0.0
 p[2,2] ≥ 0.0
 p[1,3] ≥ 0.0
 p[2,3] ≥ 0.0
 p[1,4] ≥ 0.0
 p[2,4] ≥ 0.0
 u[1,1] ≤ 1.0
 u[2,1] ≤ 1.0
 u[1,2] ≤ 1.0
 u[2,2] ≤ 1.0
 u[1,3] ≤ 1.0
 u[2,3] ≤ 1.0
 u[1,4] ≤ 1.0
 u[2,4] ≤ 1.0

Perturbation of a single input parameter

Let us vary the demand at the second time frame on both loads:

demand_values = 0.05:0.05:3.0
pvalues = map(demand_values) do di
    unit_commitment(
        [1.0, di, 1.4, 1.6], [1.0, di, 1.4, 1.6],
        [1000.0, 1500.0], [500.0, 1000.0],
        silent=true,
    )
end
pflat = [getindex.(pvalues, i) for i in eachindex(pvalues[1])]
8-element Vector{Vector{Float64}}:
 [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0  …  2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0  …  3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0]
 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  2.0999999999999996, 2.2, 2.3, 2.4000000000000004, 2.5, 2.5999999999999996, 2.7, 2.8, 2.9000000000000004, 3.0]
 [2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8  …  2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8]
 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 [3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0  …  3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0]
 [0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996  …  0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996, 0.19999999999999996]

The influence of this variation of the demand is piecewise linear on the generation at different time frames:

scatter(demand_values, pflat)
title!("Generation at different time frames and generators for a single variation")
xlims!(0.0, 3.5)

Forward Differentiation

Forward differentiation rule for the solution map of the unit commitment problem. It takes as arguments:

  1. the perturbations on the input parameters
  2. the differentiated function
  3. the primal values of the input parameters,

and returns a tuple (primal_output, perturbations), the main primal result and the perturbation propagated to this result:

function ChainRulesCore.frule(
        (_, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs),
        ::typeof(unit_commitment),
        load1_demand, load2_demand, gen_costs, noload_costs;
        optimizer=Clp.Optimizer,
        )
    # creating the UC model with a DiffOpt optimizer wrapper around Clp
    model = Model(() -> diff_optimizer(optimizer))
    # building and solving the main model
    pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model)
    energy_balance_cons = model[:energy_balance_cons]

    # Setting some perturbation of the energy balance constraints
    # Perturbations are set as MOI functions
    Δenergy_balance = [
        convert(MOI.ScalarAffineFunction{Float64}, d1 + d2)
        for (d1, d2) in zip(Δload1_demand, Δload2_demand)
    ]
    MOI.set.(
        model,
        DiffOpt.ForwardInConstraint(), energy_balance_cons,
        Δenergy_balance,
    )

    p = model[:p]
    u = model[:u]

    # setting the perturbation of the linear objective
    Δobj = sum(Δgen_costs ⋅ p[:,t] + Δnoload_costs ⋅ u[:,t] for t in size(p, 2))
    MOI.set(model, DiffOpt.ForwardInObjective(), Δobj)

    DiffOpt.forward(JuMP.backend(model))
    # querying the corresponding perturbation of the decision
    Δp = MOI.get.(model, DiffOpt.ForwardOutVariablePrimal(), p)
    return (pv, Δp.data)
end

We can now compute the perturbation of the output powers Δpv for a perturbation of the first load demand at time 2:

load1_demand = [1.0, 1.0, 1.4, 1.6]
load2_demand = [1.0, 1.0, 1.4, 1.6]
gen_costs = [1000.0, 1500.0]
noload_costs = [500.0, 1000.0]

# all input perturbations are 0
# except first load at time 2
Δload1_demand = 0 * load1_demand
Δload1_demand[2] = 1.0
Δload2_demand = 0 * load2_demand
Δgen_costs = 0 * gen_costs
Δnoload_costs = 0 * noload_costs
(pv, Δpv) = ChainRulesCore.frule(
    (nothing, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs),
    unit_commitment,
    load1_demand, load2_demand, gen_costs, noload_costs,
)
([2.0 2.0 2.8 3.0; 0.0 0.0 0.0 0.19999999999999996], [0.0 -0.9999999999999998 0.0 0.0; 0.0 1.870045459630154e-16 0.0 0.0])
Δpv
2×4 Matrix{Float64}:
 0.0  -1.0          0.0  0.0
 0.0   1.87005e-16  0.0  0.0

The result matches what we observe in the previous figure: the generation of the first generator at the second time frame (third element on the plot).

Reverse-mode differentiation of the solution map

The rrule returns the primal and a pullback. The pullback takes a seed for the optimal solution ̄p and returns derivatives with respect to each input parameter of the function.

function ChainRulesCore.rrule(
        ::typeof(unit_commitment),
        load1_demand, load2_demand, gen_costs, noload_costs;
        optimizer=Clp.Optimizer,
        silent=false)
    model = Model(() -> diff_optimizer(optimizer))
    # solve the forward UC problem
    pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model, silent=silent)
    function pullback_unit_commitment(pb)
        p = model[:p]
        u = model[:u]
        energy_balance_cons = model[:energy_balance_cons]

        MOI.set.(model, DiffOpt.BackwardInVariablePrimal(), p, pb)
        DiffOpt.backward(JuMP.backend(model))

        obj = MOI.get(model, DiffOpt.BackwardOutObjective())

        # computing derivative wrt linear objective costs
        dgen_costs = similar(gen_costs)
        dgen_costs[1] = sum(JuMP.coefficient.(obj, p[1,:]))
        dgen_costs[2] = sum(JuMP.coefficient.(obj, p[2,:]))

        dnoload_costs = similar(noload_costs)
        dnoload_costs[1] = sum(JuMP.coefficient.(obj, u[1,:]))
        dnoload_costs[2] = sum(JuMP.coefficient.(obj, u[2,:]))

        # computing derivative wrt constraint constant
        dload1_demand = JuMP.constant.(MOI.get.(model, DiffOpt.BackwardOutConstraint(), energy_balance_cons))
        dload2_demand = copy(dload1_demand)
        return (dload1_demand, dload2_demand, dgen_costs, dnoload_costs)
    end
    return (pv, pullback_unit_commitment)
end

We can set a seed of one on the power of the first generator at the second time frame and zero for all other parts of the solution:

(pv, pullback_unit_commitment) = ChainRulesCore.rrule(
    unit_commitment,
    load1_demand, load2_demand, gen_costs, noload_costs,
    optimizer=Clp.Optimizer,
    silent=true,
)
dpv = 0 * pv
dpv[1,2] = 1
dargs = pullback_unit_commitment(dpv)
(dload1_demand, dload2_demand, dgen_costs, dnoload_costs) = dargs

The sensitivities with respect to the load demands are:

dload1_demand
4-element Vector{Float64}:
 -0.0
 -0.9999999999999998
 -0.0
 -0.0
dload2_demand
4-element Vector{Float64}:
 -0.0
 -0.9999999999999998
 -0.0
 -0.0

The sensitivity of the generation is propagated to the sensitivity of both loads at the second time frame.

This example integrating ChainRules was designed with support from Invenia Technical Computing.