Edit on GitHub

Example Fitting a Gaussian Process to the profit curve (using AbstractGPs.jl)

This a example on how to use AbstractGPs.jl to fit a Gaussian Process (GP) to the profit curve of a company operating in a defined market (through the OptimalBids API).

  • For this example, we will use a market of type PowerModelsMarkets which represents an energy spot market.
  • Case instance is IEEE 118 Bus Case
  • To ilustrate a situation where the company defines the distribution of volume across it's assets in a predefined way (and reduce the decision dimentionality), the only controlable variable to maximize the company profit is a multiplicative factor of all offers.
  • We will use Optim.jl for the MLE of the tested kernels' hyperparameters.

Include Needed Packages

using OptimalBids
using OptimalBids.PowerModelsMarkets
using HiGHS # Market Clearing Solver
using JuMP: optimizer_with_attributes

using AbstractGPs
using KernelFunctions

using Optim
using ParameterHandling
using Zygote
using ParameterHandling: flatten

using Plots # For some evaluation plots at the end
using Plots.PlotMeasures
using Downloads # To download Test Cases
┌ Warning: Module DistributionsTestExt with build ID ffffffff-ffff-ffff-0000-00698dc5e021 is missing from the cache.
│ This may mean DistributionsTestExt [ffbe0ea5-a612-5ff7-aaf5-cac02eef3019] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:1948
┌ Error: Error during loading of extension DistributionsTestExt of Distributions, use `Base.retry_load_extensions()` to retry.
│   exception =
│    1-element ExceptionStack:
│    Declaring __precompile__(false) is not allowed in files that are being precompiled.
│    Stacktrace:
│      [1] _require(pkg::Base.PkgId, env::Nothing)
│        @ Base ./loading.jl:1952
│      [2] __require_prelocked(uuidkey::Base.PkgId, env::Nothing)
│        @ Base ./loading.jl:1812
│      [3] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│      [4] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│      [5] _require_prelocked
│        @ ./loading.jl:1803 [inlined]
│      [6] _require_prelocked
│        @ ./loading.jl:1802 [inlined]
│      [7] run_extension_callbacks(extid::Base.ExtensionId)
│        @ Base ./loading.jl:1295
│      [8] run_extension_callbacks(pkgid::Base.PkgId)
│        @ Base ./loading.jl:1330
│      [9] run_package_callbacks(modkey::Base.PkgId)
│        @ Base ./loading.jl:1164
│     [10] __require_prelocked(uuidkey::Base.PkgId, env::String)
│        @ Base ./loading.jl:1819
│     [11] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│     [12] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│     [13] _require_prelocked(uuidkey::Base.PkgId, env::String)
│        @ Base ./loading.jl:1803
│     [14] macro expansion
│        @ ./loading.jl:1790 [inlined]
│     [15] macro expansion
│        @ ./lock.jl:267 [inlined]
│     [16] __require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1753
│     [17] #invoke_in_world#3
│        @ ./essentials.jl:926 [inlined]
│     [18] invoke_in_world
│        @ ./essentials.jl:923 [inlined]
│     [19] require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:1746
│     [20] include
│        @ ./Base.jl:495 [inlined]
│     [21] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::String)
│        @ Base ./loading.jl:2222
│     [22] top-level scope
│        @ stdin:3
│     [23] eval
│        @ ./boot.jl:385 [inlined]
│     [24] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
│        @ Base ./loading.jl:2076
│     [25] include_string
│        @ ./loading.jl:2086 [inlined]
│     [26] exec_options(opts::Base.JLOptions)
│        @ Base ./client.jl:316
│     [27] _start()
│        @ Base ./client.jl:552
└ @ Base loading.jl:1301

Case Definition

# Read market data from IEEE 118 bus case
case_name = "pglib_opf_case118_ieee.m"
DATA_DIR = mktempdir()
case_file_path = joinpath(DATA_DIR, case_name)
Downloads.download("https://raw.githubusercontent.com/power-grid-lib/pglib-opf/01681386d084d8bd03b429abcd1ee6966f68b9a3/" * case_name, case_file_path)
network_data = PowerModels.parse_file(case_file_path)

# Measure maximum load
max_load = sum(load["pd"] for load in values(network_data["load"])) * network_data["baseMVA"]

# Pretend we are a company constructing a new set of generators in the grid.
# Choose a percentage of the total number of buses to install the new generators:
percentage_buses = 0.09

# We need the keys PowerModels uses to reference the appropriate buses in it's network data dictionary.
# First, find out all available keys:
bus_indexes = collect(keys(network_data["bus"]))
# Then, calculate number of buses that consitute the chose percent (`percentage_buses`):
num_buses = length(bus_indexes)
num_strategic_buses = ceil(Int, percentage_buses * num_buses)
# To avoid any biases let's grab some generators in the middle:
bus_indexes = bus_indexes[21:(21 + num_strategic_buses - 1)]
# Finally, add new generators to the network grid data and collect their reference keys.
generator_indexes = [
    add_generator(network_data, parse(Int, bus_idx)) for bus_idx in bus_indexes
]
11-element Vector{String}:
 "55"
 "56"
 "57"
 "58"
 "59"
 "60"
 "61"
 "62"
 "63"
 "64"
 "65"

OptimalBids API


# Define market
market = build_market(
    PowerModelsMarket,
    network_data,
    generator_indexes,
    optimizer_with_attributes(HiGHS.Optimizer, "LogLevel" => 0),
)

# Relative distribution of offers are sometimes predefined and cannot be changed at bidding time.
using Random
rng = MersenneTwister(0)
offer_weights = rand(rng, num_strategic_buses)
offer_weights = offer_weights/ sum(offer_weights)

# However, the decision maker is allowed to increase all bids evenly:
min_total_volume = 0.0
max_total_volume = 655.0
range_mul_factor = min_total_volume:1.0:max_total_volume
bid_range = [offer_weights .* [i] for i in range_mul_factor]
p_curve = profit_curve!(market, bid_range)

Optim API

function bfgs(objective, initial_params; verbose=false)
   training_results = Optim.optimize(
        objective,
        θ -> only(Zygote.gradient(objective, θ)),
        initial_params,
        LBFGS(
            alphaguess = Optim.LineSearches.InitialStatic(scaled=true),
            linesearch = Optim.LineSearches.BackTracking(),
        ),
        Optim.Options(show_trace = verbose);
        inplace=false,
    )

   return (training_results.minimum, training_results.minimizer)
end

function multi_start_hyperopt(target, algorithm, num_iterations, initial_initial_params; update_initials = (θ_last, θ_best, obj_diff) -> rand(length(θ_last)), verbose=false, kwargs...)
    θ_initials = Array{Float64, 2}(undef, num_iterations, length(initial_initial_params));
    vistited_best_fs = Array{Float64}(undef, num_iterations);
    best_θ = fill(1/length(initial_initial_params), length(initial_initial_params))
    best_f = target(best_θ)
    last_θ = best_θ
    last_f = best_f
    for i = 1:num_iterations
        try
            (last_f, last_θ) = algorithm(target, update_initials(last_θ, best_θ, last_f-best_f); verbose=verbose, kwargs...)
        catch
            println("SHIT")
        end
        θ_initials[i,:] .= last_θ
        if last_f < best_f
            best_f = last_f
            best_θ = last_θ
        end
        vistited_best_fs[i] = best_f
    end
    return best_f, best_θ, θ_initials, vistited_best_fs
end
multi_start_hyperopt (generic function with 1 method)

AbstractGPs Interface

Data


# Randomly choose observed data
number_samples = 9
idx = [1; rand(2:size(range_mul_factor,1), number_samples)]
x = collect(range_mul_factor)[idx]
y = p_curve[idx]

plt_comp = plot(collect(range_mul_factor) * 100 / max_load, p_curve,
    label="Range Evaluation",
    ylabel="Profit (USD)",
    size=(900, 600), 
    title="Kernel Comparison", 
    xlabel="Bid Volume (% Market Share)",
    legend=:outertopright,
    color="black",
    width=3,
    left_margin=10mm,
    bottom_margin=10mm,
);
scatter!(plt_comp, x .* 100 ./ max_load, y; label="Data");

plt_comp

Matern32Kernel


## Define GP prior with Matern32Kernel
flat_initial_params, unflatten = flatten((;
    var_kernel = bounded(10.0, 1.0, 1e8),
    λ = bounded(0.005, 0.0, 0.01),
))

# Construct a function to unpack flattened parameters and pull out the raw values.
unpack = ParameterHandling.value ∘ unflatten
params = unpack(flat_initial_params)

function build_gp(params)
    return GP(params.var_kernel * Matern52Kernel() ∘ ScaleTransform(params.λ))
end

# Define MLE objective
function objective(params)
    f = build_gp(unpack(params))
    return -logpdf(f(x, 0.001), y)
end

# Optimise using Optim
(best_f, best_θ, θ_initials, vistited_best_fs) = multi_start_hyperopt(objective, bfgs, 100, flat_initial_params)

# Extracting the optimal values of the parameters
optimal_params = unpack(best_θ)

# Final GP
f = build_gp(optimal_params)

# Finite projection of `f` at inputs `x`.
# Added Gaussian noise with variance 0.001.
fx = f(x, 0.001)

# Exact posterior given `y`. This is another GP.
p_fx = posterior(fx, y)

# Plot posterior.
plot!(plt_comp, collect(range_mul_factor) * 100 / max_load, p_fx(collect(range_mul_factor)); label="Posterior - " * string(Matern32Kernel))

RationalKernel


## Define GP prior with RationalKernel
flat_initial_params, unflatten = flatten((;
    var_kernel = bounded(10.0, 6.0, 1e8),
    λ = bounded(0.005, 0.0, 0.01),
    α = positive(0.05),
))

# Construct a function to unpack flattened parameters and pull out the raw values.
unpack = ParameterHandling.value ∘ unflatten
params = unpack(flat_initial_params)

function build_gp(params)
    return GP(params.var_kernel * RationalKernel(α=params.α) ∘ ScaleTransform(params.λ))
end

# Define MLE objective
function objective(params)
    f = build_gp(unpack(params))
    return -logpdf(f(x, 0.001), y)
end

# Optimise using Optim
(best_f, best_θ, θ_initials, vistited_best_fs) = multi_start_hyperopt(objective, bfgs, 100, flat_initial_params)

# Extracting the optimal values of the parameters
optimal_params = unpack(best_θ)

# Final GP
f = build_gp(optimal_params)

# Finite projection of `f` at inputs `x`.
# Added Gaussian noise with variance 0.001.
fx = f(x, 0.001)

# Exact posterior given `y`. This is another GP.
p_fx = posterior(fx, y)

# Plot posterior.
plot!(plt_comp, collect(range_mul_factor) * 100 / max_load, p_fx(collect(range_mul_factor)); label="Posterior - " * string(RationalKernel))

FBMKernel


## Define GP prior with FBMKernel
flat_initial_params, unflatten = flatten((;
    var_kernel = bounded(10.0, 6.0, 1e8),
    λ = bounded(0.005, 0.0, 0.01),
    h = bounded(0.05, 0.0, 1.0),
))

# Construct a function to unpack flattened parameters and pull out the raw values.
unpack = ParameterHandling.value ∘ unflatten
params = unpack(flat_initial_params)

function build_gp(params)
    return GP(params.var_kernel * FBMKernel(h=params.h) ∘ ScaleTransform(params.λ))
end

# Define MLE objective
function objective(params)
    f = build_gp(unpack(params))
    return -logpdf(f(x, 0.001), y)
end

# Optimise using Optim
(best_f, best_θ, θ_initials, vistited_best_fs) = multi_start_hyperopt(objective, bfgs, 100, flat_initial_params)

# Extracting the optimal values of the parameters
optimal_params = unpack(best_θ)

# Final GP
f = build_gp(optimal_params)

# Finite projection of `f` at inputs `x`.
# Added Gaussian noise with variance 0.001.
fx = f(x, 0.001)

# Exact posterior given `y`. This is another GP.
p_fx = posterior(fx, y)

# Plot posterior.
plot!(plt_comp, collect(range_mul_factor) * 100 / max_load, p_fx(collect(range_mul_factor)); label="Posterior - " * string(FBMKernel))