Run UC - ED simulation with RH and lmp exports

market_sim.jl

using PowerSystems
const PSY = PowerSystems
using PowerSystems: UtilsData
using PowerSimulations
const PSI = PowerSimulations
using Dates
using TimeSeries
using Random
using Cbc
using GLPK
using Logging
using PowerGraphics
using CSV
using DataFrames
using DataStructures
plotlyjs()

const initial_time = DateTime("2020-09-01")
const DISPATCH_INCREASE = 2.0
const FIX_DECREASE = 0.3

cbc_solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.05)
glpk_solver = optimizer_with_attributes(GLPK.Optimizer)

rts_dir = download("https://github.com/GridMod/RTS-GMLC", "master", pwd())
rts_src_dir = joinpath(rts_dir, "RTS_Data", "SourceData")
rts_siip_dir = joinpath(rts_dir, "RTS_Data", "FormattedData", "SIIP");

rawsys = PSY.PowerSystemTableData(
    rts_src_dir,
    100.0,
    joinpath(rts_siip_dir, "user_descriptors.yaml"),
    timeseries_metadata_file = joinpath(rts_siip_dir, "timeseries_pointers.json"),
    generator_mapping_file = joinpath(rts_siip_dir, "generator_mapping.yaml"),
)

sys_DA = System(rawsys; time_series_resolution = Dates.Hour(1))
sys_RT = System(rawsys; time_series_resolution = Dates.Minute(5))

# Add area renewable energy forecasts for RT model
area_mapping = get_aggregation_topology_mapping(Area, sys_RT)
for (k, buses_in_area) in area_mapping
    k == "1" && continue
    remove_component!(sys_RT, get_component(Area, sys_RT, k))
    for b in buses_in_area
        PSY.set_area!(b, get_component(Area, sys_RT, "1"))
    end
end

for sys in [sys_DA, sys_RT]
    # Adjust Reserve Provisions
    # Remove Flex Reserves
    set_units_base_system!(sys, "SYSTEM_BASE")
    res_up = get_component(VariableReserve{ReserveUp}, sys, "Flex_Up")
    remove_component!(sys, res_up)
    res_dn = get_component(VariableReserve{ReserveDown}, sys, "Flex_Down")
    remove_component!(sys, res_dn)
    reg_reserve_up = get_component(VariableReserve, sys, "Reg_Up")
    mult = (sys == sys_DA) ? 1.5 : 1.0
    set_requirement!(reg_reserve_up, mult*get_requirement(reg_reserve_up))
    reg_reserve_dn = get_component(VariableReserve, sys, "Reg_Down")
    mult = (sys == sys_DA) ? 1.5 : 1.0
    set_requirement!(reg_reserve_dn, mult*get_requirement(reg_reserve_dn))
    spin_reserve_R1 = get_component(VariableReserve, sys, "Spin_Up_R1")
    spin_reserve_R2 = get_component(VariableReserve, sys, "Spin_Up_R2")
    spin_reserve_R3 = get_component(VariableReserve, sys, "Spin_Up_R3")

    for g in get_components(ThermalStandard, sys, x -> get_prime_mover(x) in [PrimeMovers.CT, PrimeMovers.CC])
        if get_fuel(g) == ThermalFuels.DISTILLATE_FUEL_OIL
            remove_component!(sys, g)
            continue
        end
        g.operation_cost.shut_down = g.operation_cost.start_up/2.0

        if PSY.get_base_power(g) > 3
            #remove_service!(g, reg_reserve_dn)
            #remove_service!(g, reg_reserve_up)
            continue
        end
        clear_services!(g)
        add_service!(g, reg_reserve_dn)
        add_service!(g, reg_reserve_up)

        if get_prime_mover(g) == PrimeMovers.CT
            set_status!(g, false)
            set_active_power!(g, 0.0)
            old_pwl_array = get_variable(get_operation_cost(g)) |> get_cost
            new_pwl_array = similar(old_pwl_array)
            for (ix, tup) in enumerate(old_pwl_array)
                if ix ∈ [1, length(old_pwl_array)]
                    cost_noise = 50.0*rand()
                    new_pwl_array[ix] = ((tup[1] + cost_noise), tup[2])
                else
                    try_again = true
                    while try_again
                        cost_noise = 50.0*rand()
                        power_noise = 0.01*rand()
                        slope_previous = ((tup[1] + cost_noise) - old_pwl_array[ix - 1][1])/((tup[2] - power_noise) - old_pwl_array[ix - 1][2])
                        slope_next = (- (tup[1] + cost_noise) + old_pwl_array[ix + 1][1])/(-(tup[2] - power_noise) + old_pwl_array[ix + 1][2])
                        new_pwl_array[ix] = ((tup[1] + cost_noise), (tup[2] - power_noise))
                        try_again = slope_previous > slope_next
                    end
                end
            end
            get_variable(get_operation_cost(g)).cost = new_pwl_array
        end

    end
    #Remove units that make no sense to include
    names = [
        "114_SYNC_COND_1",
        "314_SYNC_COND_1",
        "313_STORAGE_1",
        "214_SYNC_COND_1",
        "212_CSP_1",
    ]
    for d in get_components(Generator, sys, x -> x.name ∈ names)
        remove_component!(sys, d)
    end
    for br in get_components(DCBranch, sys)
        remove_component!(sys, br)
    end
    for d in get_components(Storage, sys)
        remove_component!(sys, d)
    end

    # Remove large Coal and Nuclear from reserves
    for d in
        get_components(ThermalStandard, sys, x -> (occursin(r"STEAM|NUCLEAR", get_name(x))))
        get_fuel(d) == ThermalFuels.COAL && (get_ramp_limits(d) = (up = 0.001, down = 0.001))
        if get_fuel(d) == ThermalFuels.DISTILLATE_FUEL_OIL
            remove_component!(sys, d)
            continue
        end
        get_operation_cost(d).shut_down = get_operation_cost(d).start_up/2.0
        if get_rating(d) < 3
            set_status!(d, false)
            #clear_services!(d)
            #reserve_hydro && add_service!(d, reg_reserve_up)
            #reserve_hydro && add_service!(d, reg_reserve_dn)
            #add_service!(d, spin_reserve_R1)
            set_status!(d, false)
            set_active_power!(d, 0.0)
            continue
        end
        clear_services!(d)
        get_operation_cost(d).shut_down = get_operation_cost(d).start_up/2.0
        if get_fuel(d) == ThermalFuels.NUCLEAR
            get_ramp_limits(d) = (up = 0.0, down = 0.0)
            get_time_limits(d) = (up = 4380.0, down = 4380.0)
        end

    end

    for d in get_components(RenewableDispatch, sys)
        clear_services!(d)
    end

    # Add Hydro to regulation reserves
    for d in get_components(HydroEnergyReservoir, sys)
       remove_component!(sys, d)
    end

    for d in get_components(HydroDispatch, sys)
        clear_services!(d)
    end

    for g in get_components(RenewableDispatch, sys, x -> get_prime_mover(x) == PrimeMovers.PVe)
        rat_ = get_rating(g)
        set_rating!(g, DISPATCH_INCREASE*rat_)
    end

    for g in get_components(RenewableFix, sys, x -> get_prime_mover(x) == PrimeMovers.PVe)
        rat_ = get_rating(g)
        set_rating!(g, FIX_DECREASE*rat_)
    end
end

transform_single_time_series!(sys_DA, 48, Hour(24))
transform_single_time_series!(sys_RT, 12, Minute(15))

to_json(sys_RT, joinpath(pwd(), "RTS_5min.json"), force = true)
to_json(sys_DA, joinpath(pwd(), "RTS_1hr.json"), force = true)

sys_RT = System("RTS_5min.json")
sys_DA = System("RTS_1hr.json")

rt_ptdf = PTDF(sys_RT);
devices_uc = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :HydroROR => DeviceModel(HydroDispatch, HydroDispatchRunOfRiver),
    :Hydro => DeviceModel(HydroEnergyReservoir, FixedOutput),
    :RenFx => DeviceModel(RenewableFix, FixedOutput),
)
template_uc = template_unit_commitment(devices = devices_uc)
devices_ed = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalRampLimited),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :HydroROR => DeviceModel(HydroDispatch, HydroDispatchRunOfRiver),
    :Hydro => DeviceModel(HydroEnergyReservoir, FixedOutput),
    :RenFx => DeviceModel(RenewableFix, FixedOutput),
)
template_ed =
    template_economic_dispatch(network = StandardPTDFModel,
    devices = devices_ed,
    services = template_uc.services)

stages_definition = Dict(
    "UC" => Stage(template_uc, sys_DA, cbc_solver; system_to_file = false),
    "ED" => Stage(
        template_ed,
        sys_RT,
        GLPK_solver;
        PTDF = rt_ptdf,
        system_to_file = false,
        services_slack_variables = true,
        balance_slack_variables = true,
        constraint_duals = [:CopperPlateBalance, :network_flow]
    ),
)

feedforward_chronologies = Dict(
    ("UC" => "ED") => Synchronize(periods = 24),
)

feedforward = Dict(
    ("ED", :devices, :Generators) =>
        SemiContinuousFF(binary_source_stage = ON, affected_variables = [ACTIVE_POWER]),
    )

order = Dict(
    1 => "UC",
    2 => "ED"
    )

horizons = Dict(
                "UC" => 48,
                "ED" => 12
                )
intervals = Dict(
    "UC" => (Hour(24), Consecutive()),
    "ED" => (Minute(15), RecedingHorizon()),
)

DA_RT_sequence = SimulationSequence(
    step_resolution = Hour(24),
    order = order,
    horizons = horizons,
    intervals = intervals,
    ini_cond_chronology = InterStageChronology(),
    feedforward_chronologies = feedforward_chronologies,
    feedforward = feedforward,
)


sim = Simulation(
    name = "rts-test",
    steps = 1,
    stages = stages_definition,
    stages_sequence = DA_RT_sequence,
    simulation_folder = pwd(),
    initial_time = initial_time,
)

build!(sim, console_level = Logging.Warn, recorders = [:simulation])

sim_results = execute!(sim)

uc_results = load_simulation_results("rts-test/1", "UC")
fuel_plot(uc_results, sys_DA, stair = true, load=true, curtailment = true)

ed_results = load_simulation_results("rts-test/1", "ED")
fuel_plot(ed_results, sys_RT, stair = true, load = true, curtailment = true)