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/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
new file mode 100644
index 000000000..f4c7e708e
--- /dev/null
+++ b/examples/chainrules.jl
@@ -0,0 +1,156 @@
+using JuMP
+using Clp
+using DiffOpt
+using Test
+using ChainRulesCore
+
+# 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
+    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)
+    @assert termination_status(model) == MOI.OPTIMAL
+    # converting to dense matrix
+    return JuMP.value.(p.data)
+end
+
+@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)
+    # 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]
+
+    # 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])
+        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))
+    # querying the corresponding perturbation of the decision
+    Δp = MOI.get.(model, DiffOpt.ForwardOut{MOI.VariablePrimal}(), p)
+    return (pv, Δp.data)
+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
+@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)))
+    # 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]
+        u = model[:u]
+        energy_balance_cons = model[:energy_balance_cons]
+
+        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,:]))
+
+        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,:]))
+        
+        # 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)
+    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)))
+@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