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 as MOI
import Plots

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.