Add JuMP-like API

Project.toml

@@ -23,5 +23,5 @@ JuMP = "1"
 LazyArrays = "0.21, 0.22, 1"
 MathOptInterface = "1.18"
 MathOptSetDistances = "0.2.9"
-ParametricOptInterface = "0.9.0"
+ParametricOptInterface = "0.12.1"
 julia = "1.6"

README.md

@@ -33,15 +33,12 @@ examples, tutorials, and an API reference.
 
 ### DiffOpt-JuMP API with `Parameters`
 
+Here is an example with a Parametric **Linear Program**:
+
 ```julia
 using JuMP, DiffOpt, HiGHS
 
-model = Model(
-    () -> DiffOpt.diff_optimizer(
-        HiGHS.Optimizer;
-        with_parametric_opt_interface = true,
-    ),
-)
+model = DiffOpt.quadratic_diff_model(HiGHS.Optimizer)
 set_silent(model)
 
 p_val = 4.0
@@ -64,9 +61,9 @@ optimize!(model)
 
 # differentiate w.r.t. p
 direction_p = 3.0
-MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), Parameter(direction_p))
+DiffOpt.set_forward_parameter(model, p, direction_p)
 DiffOpt.forward_differentiate!(model)
-@show MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) == direction_p * 3 / pc_val
+@show DiffOpt.get_forward_variable(model, x) == direction_p * 3 / pc_val
 
 # update p and pc
 p_val = 2.0
@@ -82,45 +79,30 @@ optimize!(model)
 DiffOpt.empty_input_sensitivities!(model)
 # differentiate w.r.t. pc
 direction_pc = 10.0
-MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), Parameter(direction_pc))
+DiffOpt.set_forward_parameter(model, pc, direction_pc)
 DiffOpt.forward_differentiate!(model)
-@show abs(MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) -
+@show abs(DiffOpt.get_forward_variable(model, x) -
     -direction_pc * 3 * p_val / pc_val^2) < 1e-5
 
 # always a good practice to clear previously set sensitivities
 DiffOpt.empty_input_sensitivities!(model)
 # Now, reverse model AD
 direction_x = 10.0
-MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_x)
+DiffOpt.set_reverse_variable(model, x, direction_x)
 DiffOpt.reverse_differentiate!(model)
-@show MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) == MOI.Parameter(direction_x * 3 / pc_val)
-@show abs(MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(pc)).value -
-    -direction_x * 3 * p_val / pc_val^2) < 1e-5
+@show DiffOpt.get_reverse_parameter(model, p) == direction_x * 3 / pc_val
+@show DiffOpt.get_reverse_parameter(model, pc) == -direction_x * 3 * p_val / pc_val^2
 ```
 
-### Low level DiffOpt-JuMP API:
-
-A brief example:
+Available models:
+* `DiffOpt.quadratic_diff_model`: Quadratic Programs (QP) and Linear Programs
+(LP)
+* `DiffOpt.conic_diff_model`: Conic Programs (CP) and Linear Programs (LP)
+* `DiffOpt.nonlinear_diff_model`: Nonlinear Programs (NLP), Quadratic Program
+(QP) and Linear Programs (LP)
+* `DiffOpt.diff_model`: Nonlinear Programs (NLP), Conic Programs (CP),
+Quadratic Programs (QP) and Linear Programs (LP)
 
-```julia
-using JuMP, DiffOpt, HiGHS
-# Create a model using the wrapper
-model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
-# Define your model and solve it
-@variable(model, x)
-@constraint(model, cons, x >= 3)
-@objective(model, Min, 2x)
-optimize!(model)
-# Choose the problem parameters to differentiate with respect to, and set their
-# perturbations.
-MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 1.0)
-# Differentiate the model
-DiffOpt.reverse_differentiate!(model)
-# fetch the gradients
-grad_exp = MOI.get(model, DiffOpt.ReverseConstraintFunction(), cons)  # -3 x - 1
-constant(grad_exp)        # -1
-coefficient(grad_exp, x)  # -3
-```
 
 ## Citing DiffOpt.jl
 

docs/Project.toml

@@ -9,6 +9,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MLDatasets = "eb30cadb-4394-5ae3-aed4-317e484a6458"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

docs/src/examples/Thermal_Generation_Dispatch_Example_new.jl

@@ -0,0 +1,191 @@
+# # Thermal Generation Dispatch Example
+
+#md # [![](https://img.shields.io/badge/GitHub-100000?style=for-the-badge&logo=github&logoColor=white)](@__REPO_ROOT_URL__/examples/Thermal_Generation_Dispatch_Example.jl)
+
+# This example illustrates the sensitivity analysis of thermal generation dispatch problem.
+
+# This problem can be described as the choice of thermal generation `g` given a demand `d`, a price for thermal generation `c` and a penalty price `c_{ϕ}` for any demand not attended ϕ.
+
+# ```math
+# \begin{split}
+# \begin{array} {ll}
+# \mbox{minimize} & \sum_{i=1}^{N} c_{i} g_{i} + c_{\phi} \phi \\
+# \mbox{s.t.} & g_{i} \ge 0 \quad i=1..N  \\
+#             & g_{i} \le G_{i} \quad i=1..N  \\
+#             & \sum_{i=1}^{N} g_{i} + \phi = d\\
+# \end{array}
+# \end{split}
+# ```
+# where
+# - `G_{i}` is the maximum possible generation for a thermal generator `i`
+
+# ## Define and solve the Thermal Dispatch Problem
+
+# First, import the libraries.
+
+using Test
+using JuMP
+import DiffOpt
+import LinearAlgebra: dot
+import HiGHS
+import MathOptInterface as MOI
+import Plots
+
+# Define the model that will be construct given a set of parameters.
+
+function generate_model(
+    d_data::Float64;
+    g_sup::Vector{Float64},
+    c_g::Vector{Float64},
+    c_ϕ::Float64,
+)
+    ## Creation of the Model and Parameters
+    model = DiffOpt.quadratic_diff_model(HiGHS.Optimizer)
+    set_silent(model)
+    I = length(g_sup)
+
+    ## Variables
+    @variable(model, g[i in 1:I] >= 0.0)
+    @variable(model, ϕ >= 0.0)
+
+    ## Parameters
+    @variable(model, d in Parameter(d_data))
+
+    ## Constraints
+    @constraint(model, limit_constraints_sup[i in 1:I], g[i] <= g_sup[i])
+    @constraint(model, demand_constraint, sum(g) + ϕ == d)
+
+    ## Objectives
+    @objective(model, Min, dot(c_g, g) + c_ϕ * ϕ)
+
+    ## Solve the model
+    optimize!(model)
+
+    ## Return the solved model
+    return model
+end
+
+# Define the functions that will get the primal values `g` and `\phi` and sensitivity analysis of the demand `dg/dd` and `d\phi/dd` from a optimized model.
+
+function diff_forward(model::Model, ϵ::Float64 = 1.0)
+    ## Initialization of parameters and references to simplify the notation
+    vect_ref = [model[:g]; model[:ϕ]]
+
+    ## Get the primal solution of the model
+    vect = value.(vect_ref)
+
+    ## Reset the sensitivities of the model
+    DiffOpt.empty_input_sensitivities!(model)
+
+    ## Pass the perturbation to the DiffOpt Framework
+    DiffOpt.set_forward_parameter(model, model[:d], ϵ)
+
+    ## Compute the derivatives with the Forward Mode
+    DiffOpt.forward_differentiate!(model)
+
+    ## Get the derivative of the model
+    dvect = DiffOpt.get_forward_variable.(model, vect_ref)
+
+    ## Return the values as a vector
+    return [vect; dvect]
+end
+
+function diff_reverse(model::Model, ϵ::Float64 = 1.0)
+    ## Initialization of parameters and references to simplify the notation
+    vect_ref = [model[:g]; model[:ϕ]]
+
+    ## Get the primal solution of the model
+    vect = value.(vect_ref)
+
+    ## Set variables needed for the DiffOpt Backward Framework
+    dvect = Array{Float64,1}(undef, length(vect_ref))
+
+    ## Loop for each primal variable
+    for i in 1:I+1
+        ## Reset the sensitivities of the model
+        DiffOpt.empty_input_sensitivities!(model)
+
+        ## Pass the perturbation to the DiffOpt Framework
+        DiffOpt.set_reverse_variable.(model, vect_ref[i], ϵ)
+
+        ## Compute the derivatives with the Forward Mode
+        DiffOpt.reverse_differentiate!(model)
+
+        ## Get the derivative of the model
+        dvect[i] = DiffOpt.get_reverse_parameter(model, model[:d])
+    end
+
+    ## Return the values as a vector
+    return [vect; dvect]
+end
+
+# Initialize of Parameters
+
+g_sup = [10.0, 20.0, 30.0]
+I = length(g_sup)
+d = 0.0:0.1:80
+d_size = length(d)
+c_g = [1.0, 3.0, 5.0]
+c_ϕ = 10.0;
+
+# Generate models for each demand `d`
+models = generate_model.(d; g_sup = g_sup, c_g = c_g, c_ϕ = c_ϕ);
+
+# Get the results of models with the DiffOpt Forward and Backward context
+
+result_forward = diff_forward.(models)
+optimize!.(models)
+result_reverse = diff_reverse.(models);
+
+# Organization of results to plot
+# Initialize data_results that will contain every result
+data_results = Array{Float64,3}(undef, 2, d_size, 2 * (I + 1));
+
+# Populate the data_results array
+for k in 1:d_size
+    data_results[1, k, :] = result_forward[k]
+    data_results[2, k, :] = result_reverse[k]
+end
+
+# ## Results with Plot graphs
+# ### Results for the forward context
+# Result Primal Values:
+Plots.plot(
+    d,
+    data_results[1, :, 1:I+1];
+    title = "Generation by Demand",
+    label = ["Thermal Generation 1" "Thermal Generation 2" "Thermal Generation 3" "Generation Deficit"],
+    xlabel = "Demand [unit]",
+    ylabel = "Generation [unit]",
+)
+
+# Result Sensitivity Analysis:
+Plots.plot(
+    d,
+    data_results[1, :, I+2:2*(I+1)];
+    title = "Sensitivity of Generation by Demand",
+    label = ["T. Gen. 1 Sensitivity" "T. Gen. 2 Sensitivity" "T. Gen. 3 Sensitivity" "Gen. Deficit Sensitivity"],
+    xlabel = "Demand [unit]",
+    ylabel = "Sensitivity [-]",
+)
+
+# ### Results for the reverse context
+# Result Primal Values:
+Plots.plot(
+    d,
+    data_results[2, :, 1:I+1];
+    title = "Generation by Demand",
+    label = ["Thermal Generation 1" "Thermal Generation 2" "Thermal Generation 3" "Generation Deficit"],
+    xlabel = "Demand [unit]",
+    ylabel = "Generation [unit]",
+)
+
+# Result Sensitivity Analysis:
+Plots.plot(
+    d,
+    data_results[2, :, I+2:2*(I+1)];
+    title = "Sensitivity of Generation by Demand",
+    label = ["T. Gen. 1 Sensitivity" "T. Gen. 2 Sensitivity" "T. Gen. 3 Sensitivity" "Gen. Deficit Sensitivity"],
+    xlabel = "Demand [unit]",
+    ylabel = "Sensitivity [-]",
+)