Skip to content

Commit

Permalink
yc_test_scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
Yonghong Chen committed Oct 11, 2024
1 parent 0d654ac commit 2b67507
Show file tree
Hide file tree
Showing 5 changed files with 1,239 additions and 0 deletions.
23 changes: 23 additions & 0 deletions Multi-Stage/Project.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
HydroPowerSimulations = "fc1677e0-6ad7-4515-bf3a-bd6bf20a0b1b"
InfrastructureSystems = "2cd47ed4-ca9b-11e9-27f2-ab636a7671f1"
Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PowerGraphics = "5f7eddb3-86b1-49e8-a95a-ddc0f45eea41"
PowerNetworkMatrices = "bed98974-b02a-5e2f-9fe0-a103f5c450dd"
PowerSimulations = "e690365d-45e2-57bb-ac84-44ba829e73c4"
PowerSimulationsDecomposition = "bed98974-b02a-5e2f-9ee0-a103f5c450dd"
PowerSystemCaseBuilder = "f00506e0-b84f-492a-93c2-c0a9afc4364e"
PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
ShiftedArrays = "1277b4bf-5013-50f5-be3d-901d8477a67a"
StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
StorageSystemsSimulations = "e2f1a126-19d0-4674-9252-42b2384f8e3c"
TimeSeries = "9e3dc215-6440-5c97-bce1-76c03772f85e"
Xpress = "9e70acf3-d6c9-5be6-b5bd-4e2c73e3e054"
Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
268 changes: 268 additions & 0 deletions Multi-Stage/YC_areaptdf_DA-areasplitptdf_DA.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,268 @@
using Pkg
Pkg.activate(".")
using Revise

ENV["XPRESSDIR"] = "C:\\xpressmp"
ENV["XPAUTH_PATH"] = "C:\\xpressmp\\bin"

using HydroPowerSimulations
using PowerSimulations
using PowerSystemCaseBuilder
using InfrastructureSystems
using StorageSystemsSimulations
using PowerSimulationsDecomposition
using PowerSystems
# using PowerGraphics

using Xpress
using JuMP
using Logging
using Dates
using DataFrames

#include("coordination_EnergyOnly.jl")
#include("GlobalM2M.jl")
#import StatsPlots
#import Plots
include("YC_test_function.jl")
#using Plots
#gr()
const PSY = PowerSystems
const IF = InfrastructureSystems
const PSB = PowerSystemCaseBuilder

NT=5
selected_line=["CA-1", "CB-1", "AB1"]
monitoredlined_line=["A28"]
limit=20
use_monitoredline=1

sys = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
sys2 = build_system(PSISystems, "modified_RTS_GMLC_DA_sys")
transform_single_time_series!(sys, Hour(NT), Hour(NT))
transform_single_time_series!(sys2, Hour(NT), Hour(NT))
#transform_single_time_series!(sys, Hour(1), Hour(1))
#transform_single_time_series!(sys2, Hour(1), Hour(1))
#enforced_region = ["1", "2", "3"] # first one is MRTO, second one is M2M coordinator, third one in RT is loopflow under M2M
enforced_region=["1","3","2"]
#enforced_region=["2","3","1"]
ptdf=PTDF(sys)

buildsubsystem(sys2, enforced_region, union(selected_line,monitoredlined_line))

exchange_1_2 = AreaInterchange(;
name="1_2", available=true, active_power_flow=0.0, from_area=get_component(Area, sys, "1"), to_area=get_component(Area, sys, "2"),
flow_limits=(from_to=99999, to_from=99999),
)
add_component!(sys, exchange_1_2)

exchange_1_3 = AreaInterchange(;
name="1_3", available=true, active_power_flow=0.0, from_area=get_component(Area, sys, "1"), to_area=get_component(Area, sys, "3"),
flow_limits=(from_to=99999, to_from=99999),
)
add_component!(sys, exchange_1_3)

exchange_2_3 = AreaInterchange(;
name="2_3", available=true, active_power_flow=0.0, from_area=get_component(Area, sys, "2"), to_area=get_component(Area, sys, "3"),
flow_limits=(from_to=99999, to_from=99999),
)
add_component!(sys, exchange_2_3)

exchange_1_2 = AreaInterchange(;
name="1_2", available=true, active_power_flow=0.0, from_area=get_component(Area, sys2, "1"), to_area=get_component(Area, sys2, "2"),
flow_limits=(from_to=99999, to_from=99999),
)
add_component!(sys2, exchange_1_2)

add_component_to_subsystem!(sys2, "a", exchange_1_2)
add_component_to_subsystem!(sys2, "b", exchange_1_2)

exchange_1_3 = AreaInterchange(;
name="1_3", available=true, active_power_flow=0.0, from_area=get_component(Area, sys2, "1"), to_area=get_component(Area, sys2, "3"),
flow_limits=(from_to=99999, to_from=99999),
)
add_component!(sys2, exchange_1_3)
add_component_to_subsystem!(sys2, "a", exchange_1_3)
add_component_to_subsystem!(sys2, "b", exchange_1_3)

exchange_2_3 = AreaInterchange(;
name="2_3", available=true, active_power_flow=0.0, from_area=get_component(Area, sys2, "2"), to_area=get_component(Area, sys2, "3"),
flow_limits=(from_to=99999, to_from=99999),
)
add_component!(sys2, exchange_2_3)
add_component_to_subsystem!(sys2, "a", exchange_2_3)
add_component_to_subsystem!(sys2, "b", exchange_2_3)

template_uc = ProblemTemplate(NetworkModel(AreaPTDFPowerModel; use_slacks=true))
set_device_model!(template_uc, ThermalStandard, ThermalBasicUnitCommitment)
set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
set_device_model!(template_uc, AreaInterchange, StaticBranch)

if use_monitoredline==0
set_device_model!(template_uc,
DeviceModel(Line, StaticBranchUnbounded; attributes=Dict("filter_function" => x -> get_name(x) in union(selected_line,monitoredlined_line)),
))
else
set_device_model!(template_uc,
DeviceModel(Line, StaticBranchUnbounded; attributes=Dict("filter_function" => x -> get_name(x) in union(selected_line,monitoredlined_line)),))

set_device_model!(template_uc, DeviceModel(MonitoredLine, StaticBranchUnbounded, use_slacks = true))

for b in monitoredlined_line
line = PSY.get_component(Line, sys, b)
PSY.convert_component!(sys, line, MonitoredLine)
end
end

template_uc2 = MultiProblemTemplate(NetworkModel(SplitAreaPTDFPowerModel; use_slacks=true), ["a", "b"])
set_device_model!(template_uc2, AreaInterchange, StaticBranch)
set_device_model!(template_uc2, ThermalStandard, ThermalBasicUnitCommitment)
set_device_model!(template_uc2, PowerLoad, StaticPowerLoad)

if use_monitoredline==0
set_device_model!(template_uc2,
DeviceModel(Line, StaticBranch; use_slacks=true, attributes=Dict("filter_function" => x -> get_name(x) in union(selected_line,monitoredlined_line)),))
l = PSY.get_component(ACBranch, sys2, "A28")
set_rating!(l,limit)
else
set_device_model!(template_uc2,
DeviceModel(Line, StaticBranchUnbounded; use_slacks=true, attributes=Dict("filter_function" => x -> get_name(x) in union(selected_line,monitoredlined_line)),))

#set_device_model!(template_uc2, MonitoredLine, StaticBranch)
set_device_model!(template_uc2, DeviceModel(MonitoredLine, StaticBranch, use_slacks = true))
for b in monitoredlined_line
line = PSY.get_component(Line, sys2, b)
PSY.convert_component!(sys2, line, MonitoredLine)
end
l = PSY.get_component(MonitoredLine, sys2, "A28")
set_flow_limits!(l,(from_to=limit,to_from=limit))
end
models = SimulationModels(;
decision_models=[
DecisionModel(
template_uc,
sys;
name="UC0",
optimizer=optimizer_with_attributes(
Xpress.Optimizer,
"MIPRELSTOP" => 0.00, # Set the relative mip gap tolerance
"MAXMEMORYSOFT" => 600000, # Set the maximum amount of memory the solver can use (in MB)
),
system_to_file=false,
optimizer_solve_log_print=false,
direct_mode_optimizer=true,
store_variable_names=true,
calculate_conflict=true,
),
DecisionModel(
MultiRegionProblem,
template_uc2,
sys2;
name="UC_Subsystem",
optimizer=optimizer_with_attributes(
Xpress.Optimizer,
"MIPRELSTOP" => 0.00, # Set the relative mip gap tolerance
#"MAXMEMORYSOFT" => 600000, # Set the maximum amount of memory the solver can use (in MB)
),
system_to_file=false,
initialize_model=true,
optimizer_solve_log_print=true,
direct_mode_optimizer=true,
rebuild_model=false,
store_variable_names=true,
calculate_conflict=true,
),
],
)

uc_simulation_ff = Vector{PowerSimulations.AbstractAffectFeedforward}()

FVFF_area_interchange = FixValueFeedforward(;component_type=AreaInterchange,source=FlowActivePowerVariable,affected_values=[FlowActivePowerVariable],)
push!(uc_simulation_ff, FVFF_area_interchange)

#FVFF_moniterd_line = FixValueFeedforward(;
# component_type=MonitoredLine, source=FlowActivePowerVariable,affected_values=[FlowActivePowerVariable],
#)
#push!(uc_simulation_ff, FVFF_moniterd_line)

#FVFF_line = FixValueFeedforward(;
# component_type=Line, source=FlowActivePowerVariable,affected_values=[FlowActivePowerVariable],
#)
#push!(uc_simulation_ff, FVFF_line)

sequence = SimulationSequence(;
models=models,
feedforwards=Dict(
"UC_Subsystem" => uc_simulation_ff,
),
ini_cond_chronology=InterProblemChronology(),
);

sim = Simulation(;
name="sim",
steps=1,
models=models,
sequence=sequence,
initial_time=DateTime("2020-01-01T00:00:00"),
simulation_folder=mktempdir(),
);

build_out = build!(sim; console_level=Logging.Info, serialize=false)
execute_status = execute!(sim; enable_progress_bar=true);

uc0=models.decision_models[1].internal.container
uc2b=models.decision_models[2].internal.container.subproblems["b"]
uc2a=models.decision_models[2].internal.container.subproblems["a"]

println("obj uc0,uc2a,uc2b,",objective_value(uc0.JuMPmodel),",",objective_value(uc2a.JuMPmodel),",",objective_value(uc2b.JuMPmodel),",diff,",objective_value(uc0.JuMPmodel)-objective_value(uc2a.JuMPmodel)-objective_value(uc2b.JuMPmodel))

results = SimulationResults(sim)
results_uc = get_decision_problem_results(results, "UC0")
results_rt = get_decision_problem_results(results, "UC_Subsystem")

bus_df_uc0=read_bus_df(results_uc,"A28",0)
bus_df_uc2=read_bus_df(results_rt,"A28",1)

bus_df_uc0[!,"buscheck"]=bus_df_uc0[!,"ActivePowerBalance__ACBus"]-bus_df_uc0[!,"ThermalStandard"]/100-bus_df_uc0[!,"PowerLoad"]/100
println("UC0 bus check ActivePowerBalance__ACBus<>ThermalStandard/100-PowerLoad/100 ,",filter([:t,:buscheck] => (t,buscheck)-> (t>=1) && (abs(buscheck)>0.000001), bus_df_uc0))
bus_df_uc2[!,"buscheck"]=bus_df_uc2[!,"ActivePowerBalance__ACBus"]-bus_df_uc2[!,"ThermalStandard"]/100-bus_df_uc2[!,"PowerLoad"]/100
println("UC2 bus check ActivePowerBalance__ACBus<>ThermalStandard/100-PowerLoad/100 ,",filter([:t,:buscheck] => (t,buscheck)-> (t>=1) && (abs(buscheck)>0.000001), bus_df_uc2))
println("UC0 ActivePowerBalance__ACBus - UC2 StateEstimationInjections__ACBus",sum(bus_df_uc0[!,"ActivePowerBalance__ACBus"]-bus_df_uc2[!,"StateEstimationInjections__ACBus"]))
#println(filter([:t,:bus] => (t,bus)-> (t>=1) && (bus==203), bus_df_uc2))

####### area flow check #########
areaflow_df0=combine(groupby(bus_df_uc0, [:t, :area]), [:flowcontribution, :loopflowcontribution] .=> sum)
areaflow1_df0=filter([:t,:area] => (t,area)-> (t>=1) && (area=="1"), areaflow_df0)
areaflow2_df0=filter([:t,:area] => (t,area)-> (t>=1) && (area=="2"), areaflow_df0)
areaflow3_df0=filter([:t,:area] => (t,area)-> (t>=1) && (area=="3"), areaflow_df0)

areaflow_df2=combine(groupby(bus_df_uc2, [:t, :area]), [:flowcontribution, :loopflowcontribution] .=> sum)
areaflow1_df2=filter([:t,:area] => (t,area)-> (t>=1) && (area=="1"), areaflow_df2)
areaflow2_df2=filter([:t,:area] => (t,area)-> (t>=1) && (area=="2"), areaflow_df2)
areaflow3_df2=filter([:t,:area] => (t,area)-> (t>=1) && (area=="3"), areaflow_df2)

sum(areaflow2_df0[!,"flowcontribution_sum"]-areaflow2_df2[!,"loopflowcontribution_sum"])
sum(areaflow3_df0[!,"flowcontribution_sum"]-areaflow3_df2[!,"loopflowcontribution_sum"])

println("areaflow1 check sum(uc0 flowcontribution_sum - uc2 loopflowcontribution_sum,",sum(areaflow1_df0[!,"flowcontribution_sum"]-areaflow1_df2[!,"loopflowcontribution_sum"]))
println("areaflow2 check sum(uc0 flowcontribution_sum - uc2 loopflowcontribution_sum,",sum(areaflow2_df0[!,"flowcontribution_sum"]-areaflow2_df2[!,"loopflowcontribution_sum"]))
println("areaflow3 check sum(uc0 flowcontribution_sum - uc2 loopflowcontribution_sum,",sum(areaflow3_df0[!,"flowcontribution_sum"]-areaflow3_df2[!,"loopflowcontribution_sum"]))

println("varflow_uc0,",value.(uc0.variables[InfrastructureSystems.Optimization.VariableKey{FlowActivePowerVariable, MonitoredLine}("")]["A28",:]) )
println("varflow_uc2a,",value.(uc2a.variables[InfrastructureSystems.Optimization.VariableKey{FlowActivePowerVariable, MonitoredLine}("")]["A28",:]) )
println("var_flow_uc2b,",value.(uc2b.variables[InfrastructureSystems.Optimization.VariableKey{FlowActivePowerVariable, MonitoredLine}("")]["A28",:]) )

println(areaflow132_df[:,[:flow_subsystem_a ,:flow_subsystem_b, :seflow]])

areaflow13_df=innerjoin(areaflow1_df2, areaflow3_df2, on=:t,renamecols = "_1" => "_3")
areaflow132_df=innerjoin(areaflow13_df, areaflow2_df2, on=:t,renamecols = "" => "_2")

areaflow132_df[!,"flow_subsystem_a"]=areaflow132_df[!,"flowcontribution_sum_1"]+areaflow132_df[!,"flowcontribution_sum_3"]+areaflow132_df[!,"loopflowcontribution_sum_2"]
areaflow132_df[!,"flow_subsystem_b"]=areaflow132_df[!,"loopflowcontribution_sum_1"]+areaflow132_df[!,"loopflowcontribution_sum_3"]+areaflow132_df[!,"flowcontribution_sum_2"]
areaflow132_df[!,"seflow"]=areaflow132_df[!,"flowcontribution_sum_1"]+areaflow132_df[!,"flowcontribution_sum_3"]+areaflow132_df[!,"flowcontribution_sum_2"]

#StatsPlots.@df(areaflow132_df,Plots.plot(:t,[:flow_subsystem_a :flow_subsystem_b :seflow]))

areaflow13_df0=innerjoin(areaflow1_df0, areaflow3_df0, on=:t,renamecols = "_1" => "_3")
areaflow132_df0=innerjoin(areaflow13_df0, areaflow2_df0, on=:t,renamecols = "" => "_2")
areaflow132_df0[!,"flow_uc0"]=areaflow132_df0[!,"flowcontribution_sum_1"]+areaflow132_df0[!,"flowcontribution_sum_3"]+areaflow132_df0[!,"flowcontribution_sum_2"]
Loading

0 comments on commit 2b67507

Please sign in to comment.