Thermal Generation Dispatch Example

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 ϕ.

\[\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
import Plots
const MOI = MathOptInterface
MathOptInterface

Define the model that will be construct given a set of parameters.

function generate_model(d::Float64; g_sup::Vector{Float64}, c_g::Vector{Float64}, c_ϕ::Float64)
    # Creation of the Model and Parameters
    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
    set_silent(model)
    I = length(g_sup)

    # Variables
    @variable(model, g[i in  1:I] >= 0.0)
    @variable(model, ϕ >= 0.0)

    # 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
generate_model (generic function with 1 method)

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[:ϕ]]
    I = length(model[:g])

    # Get the primal solution of the model
    vect = MOI.get.(model, MOI.VariablePrimal(), vect_ref)

    # Pass the perturbation to the DiffOpt Framework and set the context to Forward
    constraint_equation = convert(MOI.ScalarAffineFunction{Float64}, ϵ)
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), model[:demand_constraint], constraint_equation)
    DiffOpt.forward_differentiate!(model)

    # Get the derivative of the model
    dvect = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), 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[:ϕ]]
    I = length(model[:g])

    # Get the primal solution of the model
    vect = MOI.get.(model, MOI.VariablePrimal(), vect_ref)

    # Set variables needed for the DiffOpt Backward Framework
    dvect = Array{Float64, 1}(undef, I + 1)
    perturbation = zeros(I + 1)

    # Loop for each primal variable
    for i in 1:I+1
        # Set the perturbation in the Primal Variables and set the context to Backward
        perturbation[i] = ϵ
        MOI.set.(model, DiffOpt.ReverseVariablePrimal(), vect_ref, perturbation)
        DiffOpt.reverse_differentiate!(model)

        # Get the value of the derivative of the model
        dvect[i] = JuMP.constant(MOI.get(model, DiffOpt.ReverseConstraintFunction(), model[:demand_constraint]))
        perturbation[i] = 0.0
    end

    # Return the values as a vector
    return [vect;dvect]
end
diff_reverse (generic function with 2 methods)

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 [-]"
)

This page was generated using Literate.jl.