From 2f5199b443dbb39c1b239549fed528d98222b9b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sat, 1 May 2021 17:36:09 +0200 Subject: [PATCH 1/4] added ChainRules examples with forward and backward --- Project.toml | 1 + examples/chainrules.jl | 124 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 125 insertions(+) create mode 100644 examples/chainrules.jl diff --git a/Project.toml b/Project.toml index 5f8213ed6..02188968a 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.1.0" [deps] BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" +ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/examples/chainrules.jl b/examples/chainrules.jl new file mode 100644 index 000000000..da2840794 --- /dev/null +++ b/examples/chainrules.jl @@ -0,0 +1,124 @@ +using JuMP +using Clp +using DiffOpt +using Test +using ChainRulesCore + +# This script creates the JuMP problem for a small unit commitment instance + +ATOL=1e-4 +RTOL=1e-4 + +function unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer))) + ## 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 + 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 + @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 + + ## 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) + @assert termination_status(model) == MOI.OPTIMAL + # converting to dense matrix + return JuMP.value.(p.data) +end + +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]) + + +load1_demand = [1.0, 1.2, 1.4, 1.6] +load2_demand = [1.0, 1.2, 1.4, 1.6] +gen_costs = [1000.0, 1500.0] +noload_costs = [500.0, 1000.0] + +function ChainRulesCore.frule((_, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs), ::typeof(unit_commitment), load1_demand, load2_demand, gen_costs, noload_costs) + model = Model(() -> diff_optimizer(Clp.Optimizer)) + pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model) + energy_balance_cons = model[:energy_balance_cons] + MOI.set.(model, DiffOpt.ForwardIn{DiffOpt.ConstraintConstant}(), energy_balance_cons, [d1 + d2 for (d1, d2) in zip(Δload1_demand, Δload1_demand)]) + + p = model[:p] + u = model[:u] + + for t in size(p, 2) + MOI.set(model, DiffOpt.ForwardIn{DiffOpt.LinearObjective}(), p[1,t], Δgen_costs[1]) + MOI.set(model, DiffOpt.ForwardIn{DiffOpt.LinearObjective}(), p[2,t], Δgen_costs[2]) + MOI.set(model, DiffOpt.ForwardIn{DiffOpt.LinearObjective}(), u[1,t], Δnoload_costs[1]) + MOI.set(model, DiffOpt.ForwardIn{DiffOpt.LinearObjective}(), u[2,t], Δnoload_costs[2]) + end + DiffOpt.forward(JuMP.backend(model)) + Δp = MOI.get.(model, DiffOpt.ForwardOut{MOI.VariablePrimal}(), p) + return (pv, Δp.data) +end + + +Δload1_demand = 0 * load1_demand .+ 0.1 +Δload2_demand = 0 * load2_demand .+ 0.2 +Δgen_costs = 0 * gen_costs .+ 0.1 +Δnoload_costs = 0 * noload_costs .+ 0.4 +(pv, Δpv) = ChainRulesCore.frule((nothing, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs), unit_commitment, load1_demand, load2_demand, gen_costs, noload_costs) + + +function ChainRulesCore.rrule(::typeof(unit_commitment), load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer))) + + pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model) + function pullback_unit_commitment(pb) + p = model[:p] + u = model[:u] + energy_balance_cons = model[:energy_balance_cons] + + MOI.set.(model, DiffOpt.BackwardIn{MOI.VariablePrimal}(), p, pb) + DiffOpt.backward(JuMP.backend(model)) + + dgen_costs = similar(gen_costs) + dgen_costs[1] = sum(MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.LinearObjective}(), p[1,:])) + dgen_costs[2] = sum(MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.LinearObjective}(), p[2,:])) + + dnoload_costs = similar(noload_costs) + dnoload_costs[1] = sum(MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.LinearObjective}(), u[1,:])) + dnoload_costs[2] = sum(MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.LinearObjective}(), u[2,:])) + + dload1_demand = MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.ConstraintConstant}(), energy_balance_cons) + dload2_demand = copy(dload1_demand) + return (dload1_demand, dload2_demand, dgen_costs, dnoload_costs) + end + return (pv, pullback_unit_commitment) +end + +(pv, pullback_unit_commitment) = ChainRulesCore.rrule(unit_commitment, load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer))) +pullback_unit_commitment(ones(size(pv))) From dbdfbeeba5bd461e060b715fad47c931c8c59558 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Wed, 5 May 2021 23:03:45 +0200 Subject: [PATCH 2/4] example chainrule in tests --- docs/Project.toml | 1 + docs/src/index.md | 9 +++------ examples/chainrules.jl | 38 ++++++++++++++++++++++++++------------ test/runtests.jl | 1 + 4 files changed, 31 insertions(+), 18 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index ffbec23bc..f2b96fd7e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/docs/src/index.md b/docs/src/index.md index 4a2d95864..b7647925f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,10 +1,6 @@ # DiffOpt.jl -[![Build Status](https://travis-ci.org/AKS1996/DiffOpt.jl.svg?branch=master)](https://travis-ci.org/AKS1996/DiffOpt.jl) -[![Coverage Status](https://coveralls.io/repos/github/AKS1996/DiffOpt.jl/badge.svg?branch=master)](https://coveralls.io/github/AKS1996/DiffOpt.jl?branch=master) -[![AppVeyor Status](https://ci.appveyor.com/api/projects/status/github/AKS1996/DiffOpt.jl?branch=master&svg=true)](https://ci.appveyor.com/project/AKS1996/diffopt-jl) -[![Docs status](https://img.shields.io/badge/docs-dev-blue.svg)](https://aks1996.github.io/DiffOpt.jl/dev/) -[DiffOpt](https://github.com/AKS1996/JuMP.jl) is a package for differentiating convex optimization program ([JuMP.jl](https://github.com/jump-dev/JuMP.jl) or [MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) models) with respect to program parameters. Note that this package does not contain any solver. This package has two major backends, available via `backward` and `forward` methods, to differentiate models (quadratic or conic) with optimal solutions. +[DiffOpt.jl](https://github.com/jump-dev/JuMP.jl) is a package for differentiating convex optimization program ([JuMP.jl](https://github.com/jump-dev/JuMP.jl) or [MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) models) with respect to program parameters. Note that this package does not contain any solver. This package has two major backends, available via `backward` and `forward` methods, to differentiate models (quadratic or conic) with optimal solutions. !!! note Currently supports *linear programs* (LP), *convex quadratic programs* (QP) and *convex conic programs* (SDP, SOCP constraints only). @@ -17,7 +13,8 @@ DiffOpt can be installed through the Julia package manager: ``` ## Why are Differentiable optimization problems important? -Differentiable optimization is a promising field of convex optimization and has many potential applications in game theory, control theory and machine learning (specifically deep learning - refer [this video](https://www.youtube.com/watch?v=NrcaNnEXkT8) for more). Recent work has shown how to differentiate specific subclasses of convex optimization problems. But several applications remain unexplored (refer section 8 of this [really good thesis](https://github.com/bamos/thesis)). With the help of automatic differentiation, differentiable optimization can a significant impact on creating end-to-end systems for modelling a neural network, stochastic process, or a game. +Differentiable optimization is a promising field of convex optimization and has many potential applications in game theory, control theory and machine learning (specifically deep learning - refer [this video](https://www.youtube.com/watch?v=NrcaNnEXkT8) for more). +Recent work has shown how to differentiate specific subclasses of convex optimization problems. But several applications remain unexplored (refer section 8 of this [really good thesis](https://github.com/bamos/thesis)). With the help of automatic differentiation, differentiable optimization can have a significant impact on creating end-to-end differentiable systems to model neural networks, stochastic processes, or a game. ## Contributing diff --git a/examples/chainrules.jl b/examples/chainrules.jl index da2840794..3913f4ce1 100644 --- a/examples/chainrules.jl +++ b/examples/chainrules.jl @@ -4,11 +4,23 @@ using DiffOpt using Test using ChainRulesCore -# This script creates the JuMP problem for a small unit commitment instance +# This script creates the JuMP model for a small unit commitment instance +# represented in a solution map function taking parameters as arguments and returning +# the optimal solution as output. + +# The derivatives of this solution map can then be expressed in ChainRules semantics +# and implemented using DiffOpt ATOL=1e-4 RTOL=1e-4 +""" +Solution map of the problem using parameters: +- `load1_demand, load2_demand` for the demand of two nodes +- `gen_costs` is the vector of generator costs +- `noload_costs` is the vector of fixed activation costs of the generators, +and returning the optimal output power `p`. +""" function unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer))) ## Problem data unit_codes = [1, 2] # Generator identifiers @@ -58,14 +70,10 @@ function unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs; mo return JuMP.value.(p.data) end -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]) - - -load1_demand = [1.0, 1.2, 1.4, 1.6] -load2_demand = [1.0, 1.2, 1.4, 1.6] -gen_costs = [1000.0, 1500.0] -noload_costs = [500.0, 1000.0] +@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]) +# Forward differentiation rule for the solution map of the unit commitment problem +# taking in input perturbations on the input parameters and returning perturbations propagated to the result function ChainRulesCore.frule((_, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs), ::typeof(unit_commitment), load1_demand, load2_demand, gen_costs, noload_costs) model = Model(() -> diff_optimizer(Clp.Optimizer)) pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model) @@ -87,15 +95,21 @@ function ChainRulesCore.frule((_, Δload1_demand, Δload2_demand, Δgen_costs, end +load1_demand = [1.0, 1.2, 1.4, 1.6] +load2_demand = [1.0, 1.2, 1.4, 1.6] +gen_costs = [1000.0, 1500.0] +noload_costs = [500.0, 1000.0] + Δload1_demand = 0 * load1_demand .+ 0.1 Δload2_demand = 0 * load2_demand .+ 0.2 Δgen_costs = 0 * gen_costs .+ 0.1 Δnoload_costs = 0 * noload_costs .+ 0.4 -(pv, Δpv) = ChainRulesCore.frule((nothing, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs), unit_commitment, load1_demand, load2_demand, gen_costs, noload_costs) - +@show (pv, Δpv) = ChainRulesCore.frule((nothing, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs), unit_commitment, load1_demand, load2_demand, gen_costs, noload_costs) +# Reverse-mode differentiation of the solution map +# The computed pullback takes a seed for the optimal solution `̄p` and returns +# derivatives wrt each input parameter. function ChainRulesCore.rrule(::typeof(unit_commitment), load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer))) - pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model) function pullback_unit_commitment(pb) p = model[:p] @@ -121,4 +135,4 @@ function ChainRulesCore.rrule(::typeof(unit_commitment), load1_demand, load2_dem end (pv, pullback_unit_commitment) = ChainRulesCore.rrule(unit_commitment, load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer))) -pullback_unit_commitment(ones(size(pv))) +@show pullback_unit_commitment(ones(size(pv))) diff --git a/test/runtests.jl b/test/runtests.jl index da4274998..c0ce05672 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,6 +27,7 @@ const RTOL = 1e-4 include(joinpath(@__DIR__, "../examples/solve-QP.jl")) include(joinpath(@__DIR__, "../examples/unit_example.jl")) include(joinpath(@__DIR__, "../examples/sensitivity-SVM.jl")) + include(joinpath(@__DIR__, "../examples/chainrules.jl")) end @testset "Generate random problems" begin From b5c56d01a00592f4d692cd1777f5390e9b1b5287 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sun, 9 May 2021 16:24:53 +0200 Subject: [PATCH 3/4] Update examples/chainrules.jl Co-authored-by: Raphael Saavedra --- examples/chainrules.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/chainrules.jl b/examples/chainrules.jl index 3913f4ce1..e2ac215c7 100644 --- a/examples/chainrules.jl +++ b/examples/chainrules.jl @@ -30,7 +30,7 @@ function unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs; mo 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 + 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 ($) From 9ce17c48eeb07e85f9d40790353ebdf572796aa6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sun, 9 May 2021 16:49:17 +0200 Subject: [PATCH 4/4] detailed comments on example --- examples/chainrules.jl | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/examples/chainrules.jl b/examples/chainrules.jl index e2ac215c7..f4c7e708e 100644 --- a/examples/chainrules.jl +++ b/examples/chainrules.jl @@ -35,8 +35,11 @@ function unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs; mo 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 + @variable(model, p[g in unit_codes, t in 1:n_periods] >= 0) # Power output (pu) ## Constraints @@ -75,14 +78,25 @@ end # Forward differentiation rule for the solution map of the unit commitment problem # taking in input perturbations on the input parameters and returning perturbations propagated to the result function ChainRulesCore.frule((_, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs), ::typeof(unit_commitment), load1_demand, load2_demand, gen_costs, noload_costs) + # creating the UC model with a DiffOpt optimizer wrapper around Clp model = Model(() -> diff_optimizer(Clp.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] - MOI.set.(model, DiffOpt.ForwardIn{DiffOpt.ConstraintConstant}(), energy_balance_cons, [d1 + d2 for (d1, d2) in zip(Δload1_demand, Δload1_demand)]) + + # Setting some perturbation of the right-hand side of the energy balance constraints + # the RHS is equal to the sum of load demands at each period. + # the corresponding perturbation are set accordingly as the set of perturbations of the two loads + MOI.set.( + model, + DiffOpt.ForwardIn{DiffOpt.ConstraintConstant}(), energy_balance_cons, + [d1 + d2 for (d1, d2) in zip(Δload1_demand, Δload1_demand)], + ) p = model[:p] u = model[:u] + # setting the perturbation of the linear objective for t in size(p, 2) MOI.set(model, DiffOpt.ForwardIn{DiffOpt.LinearObjective}(), p[1,t], Δgen_costs[1]) MOI.set(model, DiffOpt.ForwardIn{DiffOpt.LinearObjective}(), p[2,t], Δgen_costs[2]) @@ -90,6 +104,7 @@ function ChainRulesCore.frule((_, Δload1_demand, Δload2_demand, Δgen_costs, MOI.set(model, DiffOpt.ForwardIn{DiffOpt.LinearObjective}(), u[2,t], Δnoload_costs[2]) end DiffOpt.forward(JuMP.backend(model)) + # querying the corresponding perturbation of the decision Δp = MOI.get.(model, DiffOpt.ForwardOut{MOI.VariablePrimal}(), p) return (pv, Δp.data) end @@ -110,6 +125,7 @@ noload_costs = [500.0, 1000.0] # The computed pullback takes a seed for the optimal solution `̄p` and returns # derivatives wrt each input parameter. function ChainRulesCore.rrule(::typeof(unit_commitment), load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer))) + # solve the forward UC problem pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model) function pullback_unit_commitment(pb) p = model[:p] @@ -119,6 +135,7 @@ function ChainRulesCore.rrule(::typeof(unit_commitment), load1_demand, load2_dem MOI.set.(model, DiffOpt.BackwardIn{MOI.VariablePrimal}(), p, pb) DiffOpt.backward(JuMP.backend(model)) + # computing derivative wrt linear objective costs dgen_costs = similar(gen_costs) dgen_costs[1] = sum(MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.LinearObjective}(), p[1,:])) dgen_costs[2] = sum(MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.LinearObjective}(), p[2,:])) @@ -127,6 +144,7 @@ function ChainRulesCore.rrule(::typeof(unit_commitment), load1_demand, load2_dem dnoload_costs[1] = sum(MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.LinearObjective}(), u[1,:])) dnoload_costs[2] = sum(MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.LinearObjective}(), u[2,:])) + # computing derivative wrt constraint constant dload1_demand = MOI.get.(model, DiffOpt.BackwardOut{DiffOpt.ConstraintConstant}(), energy_balance_cons) dload2_demand = copy(dload1_demand) return (dload1_demand, dload2_demand, dgen_costs, dnoload_costs)