Fair reward allocations

This example showcase the capabilities of estimating fair reward allocations within energy communities using concepts from cooperative game theory. By leveraging on TheoryOfGames.jl, the package provides functionalities to compute various allocation methods, such as the Variance (Least) Core, Shapley value, the Nucleolus, and more. These methods help in distributing the benefits of cooperation among the members of the energy community in a fair manner. In the following, we showcase how to set up an energy community optimization problem and compute fair reward allocations based on the results. More information are also available at:

D. Fioriti, G. Bigi, A. Frangioni, M. Passacantando and D. Poli, "Fair Least Core: Efficient, Stable and Unique Game-Theoretic Reward Allocation in Energy Communities by Row-Generation," in IEEE Transactions on Energy Markets, Policy and Regulation, vol. 3, no. 2, pp. 170-181, June 2025, doi: 10.1109/TEMPR.2024.3495237.
D. Fioriti, A. Frangioni, D. Poli, "Optimal sizing of energy communities with fair revenue sharing and exit clauses: Value, role and business model of aggregators and users," in Applied Energy, vol. 299, 2021, 117328,doi: 10.1016/j.apenergy.2021.117328

Initialization of the model

Import the needed packages

using EnergyCommunity, JuMP, HiGHS
using TheoryOfGames
using DataFrames, StatsPlots

Create a base Energy Community example in the data folder; use the default configuration.

folder = joinpath(@__DIR__, "data")
create_example_data(folder, config_name="default")

Input file to load the structure of the energy community based on a yaml file.

input_file = joinpath(@__DIR__, "data/energy_community_model.yml");

define optimizer and options

optimizer = optimizer_with_attributes(
    HiGHS.Optimizer,
    "log_to_console"=>false,  # suppress solver output
    "ipm_optimality_tolerance"=>1e-6,  # set optimality tolerance
)

MathOptInterface.OptimizerWithAttributes(HiGHS.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("log_to_console") => false, MathOptInterface.RawOptimizerAttribute("ipm_optimality_tolerance") => 1.0e-6])

Define the Non Cooperative model

CO_Model = ModelEC(input_file, EnergyCommunity.GroupCO(), optimizer)

An Energy Community Model
Energy Community problem for a Cooperative Model
User set: ["user1", "user2", "user3"]

Enumerative method for reward allocations

Define the enumerative mode for cooperative games

enum_mode = EnumMode(
    CO_Model,
    EnergyCommunity.GroupNC(),  # Base group is the Non Cooperative group
    no_aggregator_type=EnergyCommunity.GroupNC(),  # when the aggregator is not in the community, use Non Cooperative group
)

TheoryOfGames.EnumMode(["EC", "user1", "user2", "user3"], Dict{Set{String}, Float64}(Set(["EC", "user2", "user1"]) => 9122.53161486378, Set(["EC", "user3"]) => 7.566995918750763e-10, Set(["EC", "user1"]) => 0.0, Set(["EC", "user1", "user3"]) => 9404.772014430258, Set(["EC"]) => 0.0, Set(["user3"]) => 0.0, Set(["user2", "user1", "user3"]) => 0.0, Set(["user2", "user1"]) => 0.0, Set(["user1"]) => 0.0, Set() => 0.0…))

Shapley

Calculate fair allocation using Shapley value

dst_sh = shapley_value(enum_mode)
dst_sh

Dict{String, Float64} with 4 entries:
  "EC"    => 6118.81
  "user2" => 2983.88
  "user1" => 4068.05
  "user3" => 3077.96

Nucleolus

Calculate fair allocation using Nucleolus

dst_nuc = nucleolus(enum_mode, optimizer)
dst_nuc

Dict{String, Float64} with 4 entries:
  "EC"    => 4631.83
  "user2" => 3421.97
  "user1" => 4631.83
  "user3" => 3563.09

Variance Core

Calculate fair allocation using Variance Core

dst_vc = var_in_core(enum_mode, optimizer)
dst_vc

Dict{String, Float64} with 4 entries:
  "EC"    => 4062.18
  "user2" => 4062.18
  "user1" => 4062.18
  "user3" => 4062.18

Variance Least Core

Calculatefair allocation using Variance Least Core

dst_vlc = var_least_core(enum_mode, optimizer)
dst_vlc

Dict{String, Float64} with 4 entries:
  "EC"    => 4561.27
  "user2" => 3421.97
  "user1" => 4561.27
  "user3" => 3704.21

Verify stability of the allocations

Check if the Shapley value is in the Core

sh_in_core = verify_in_core(dst_sh, enum_mode, optimizer)
println("Shapley value in Core: ", sh_in_core)

Shapley value in Core: true

Check if the Nucleolus is in the Core

nuc_in_core = verify_in_core(dst_nuc, enum_mode, optimizer)
println("Nucleolus in Core: ", nuc_in_core)

Nucleolus in Core: true

Check if the Variance Core allocation is in the Core

vc_in_core = verify_in_core(dst_vc, enum_mode, optimizer)
println("Variance Core allocation in Core: ", vc_in_core)

Variance Core allocation in Core: true

Check if the Variance Least Core allocation is in the Core

vlc_in_core = verify_in_core(dst_vlc, enum_mode, optimizer)
println("Variance Least Core allocation in Core: ", vlc_in_core)

Variance Least Core allocation in Core: true

Compare reward allocations

Create a DataFrame to compare the different allocations

df = DataFrame(
    Member = collect(keys(dst_sh)),
    Shapley = collect(values(dst_sh)),
    Nucleolus = collect(values(dst_nuc)),
    Variance_Core = collect(values(dst_vc)),
    Variance_Least_Core = collect(values(dst_vlc)),
)
println(df)

4×5 DataFrame
 Row │ Member  Shapley  Nucleolus  Variance_Core  Variance_Least_Core
     │ String  Float64  Float64    Float64        Float64
─────┼────────────────────────────────────────────────────────────────
   1 │ EC      6118.81    4631.83        4062.18              4561.27
   2 │ user2   2983.88    3421.97        4062.18              3421.97
   3 │ user1   4068.05    4631.83        4062.18              4561.27
   4 │ user3   3077.96    3563.09        4062.18              3704.21

Plot the comparison

groupedbar(
    df.Member,
    [df.Shapley df.Nucleolus df.Variance_Core df.Variance_Least_Core],
    label = ["Shapley" "Nucleolus" "Variance Core" "Variance Least Core"],
    title = "Comparison of Fair Reward Allocations",
    xlabel = "Community Members",
    ylabel = "Allocated Reward [€]",
    bar_position = :dodge,
)

Iterative method for reward allocations

For large energy communities, the enumerative method may become computationally expensive as it implies the calculation of the utility function for all possible coalitions. In such cases, an iterative method employing row-generation can be employed to estimate fair reward allocations more efficiently. This example showcases how to set up an iterative method for computing the Nucleolus allocation.

Define the iterative mode for cooperative games

iter_mode = IterMode(
    CO_Model,
    EnergyCommunity.GroupNC(),  # Base group is the Non Cooperative group
    no_aggregator_type=EnergyCommunity.GroupNC(),  # when the aggregator is not in the community, use Non Cooperative group
    optimizer=optimizer,
)

TheoryOfGames.IterMode(["EC", "user1", "user2", "user3"], EnergyCommunity.var"#utility_callback_by_subgroup#1210"{EnergyCommunity.var"#objective_callback_by_subgroup#1077"{ModelEC, EnergyCommunity.var"#objective_callback_by_subgroup#945"{JuMP.Containers.DenseAxisArray{Float64, 1, Tuple{Vector{String}}, Tuple{JuMP.Containers._AxisLookup{Dict{String, Int64}}}}}}, EnergyCommunity.var"#objective_callback_by_subgroup#945"{JuMP.Containers.DenseAxisArray{Float64, 1, Tuple{Vector{String}}, Tuple{JuMP.Containers._AxisLookup{Dict{String, Int64}}}}}}(EnergyCommunity.var"#objective_callback_by_subgroup#1077"{ModelEC, EnergyCommunity.var"#objective_callback_by_subgroup#945"{JuMP.Containers.DenseAxisArray{Float64, 1, Tuple{Vector{String}}, Tuple{JuMP.Containers._AxisLookup{Dict{String, Int64}}}}}}(An Energy Community Model
Energy Community problem for a Cooperative Model
User set: ["user1", "user2", "user3"]
, EnergyCommunity.var"#objective_callback_by_subgroup#945"{JuMP.Containers.DenseAxisArray{Float64, 1, Tuple{Vector{String}}, Tuple{JuMP.Containers._AxisLookup{Dict{String, Int64}}}}}(1-dimensional DenseAxisArray{Float64,1,...} with index sets:
    Dimension 1, ["EC", "user1", "user2", "user3"]
And data, a 4-element Vector{Float64}:
       0.0
 -563669.0832118867
 -223241.1205710986
 -435831.69036038255)), EnergyCommunity.var"#objective_callback_by_subgroup#945"{JuMP.Containers.DenseAxisArray{Float64, 1, Tuple{Vector{String}}, Tuple{JuMP.Containers._AxisLookup{Dict{String, Int64}}}}}(1-dimensional DenseAxisArray{Float64,1,...} with index sets:
    Dimension 1, ["EC", "user1", "user2", "user3"]
And data, a 4-element Vector{Float64}:
       0.0
 -563669.0832118867
 -223241.1205710986
 -435831.69036038255)), EnergyCommunity.var"#least_profitable_coalition_callback#1289"{EnergyCommunity.var"#least_profitable_coalition_callback#1288#1290"{Int64, Bool, Float64, Float64, ModelEC, Dict{Any, Any}}}(EnergyCommunity.var"#least_profitable_coalition_callback#1288#1290"{Int64, Bool, Float64, Float64, ModelEC, Dict{Any, Any}}(0, false, 0.0001, 0.0001, Core.Box(nothing), An Energy Community Model
Energy Community problem for a Cooperative Model
User set: ["user1", "user2", "user3"]
, Dict{Any, Any}())))

Calculate Variance Least Core with iterative method

dst_vlc_iter = var_least_core(iter_mode, optimizer)
dst_vlc_iter

Dict{String, Float64} with 4 entries:
  "EC"    => 4561.27
  "user2" => 3421.97
  "user1" => 4561.27
  "user3" => 3704.21

Compare the distributions from the enumerative and iterative methods

println("Difference between enumerative and iterative Variance Least Core allocations:")
for member in keys(dst_vlc)
    diff = dst_vlc[member] - dst_vlc_iter[member]
    println("Member: ", member, ", Absolute Difference [%]: ", 100*abs(diff)/dst_vlc[member])
end

Difference between enumerative and iterative Variance Least Core allocations:
Member: EC, Absolute Difference [%]: 7.95587017568954e-12
Member: user2, Absolute Difference [%]: 1.2318964536748155e-11
Member: user1, Absolute Difference [%]: 7.975809699939387e-12
Member: user3, Absolute Difference [%]: 2.4761724242403935e-11

Info

View this file on Github.

This page was generated using Literate.jl.