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 2b67507 commit 18035e6
Showing 1 changed file with 6 additions and 115 deletions.
121 changes: 6 additions & 115 deletions Multi-Stage/YC_areaptdf_DA_2stage.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ using Dates
using DataFrames

using HiGHS

include("YC_test_function.jl")
#include("coordination_EnergyOnly.jl")
#include("GlobalM2M.jl")
#import StatsPlots
Expand All @@ -50,118 +52,7 @@ 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))

function read_bus_df(results_rt,se=0)
ac = read_realized_variable(results_rt, "ActivePowerVariable__ThermalStandard")
ab = read_realized_variable(results_rt, "ActivePowerBalance__ACBus")
if se==1 al = read_realized_variable(results_rt, "StateEstimationInjections__ACBus") end
ald = read_realized_variable(results_rt, "ActivePowerTimeSeriesParameter__PowerLoad")
am = read_realized_variable(results_rt, "FlowActivePowerVariable__MonitoredLine")

ac0 = read_realized_variable(results_uc, "ActivePowerVariable__ThermalStandard")
ab0 = read_realized_variable(results_uc, "ActivePowerBalance__ACBus")
ald0 = read_realized_variable(results_uc, "ActivePowerTimeSeriesParameter__PowerLoad")
am0 = read_realized_variable(results_uc, "FlowActivePowerVariable__MonitoredLine")

bus_ActivePowerVariable__ThermalStandard=Dict()
bus_ActivePowerTimeSeriesParameter__PowerLoad=Dict()
bus_ActivePowerBalance__ACBus=Dict()
bus_StateEstimationInjections__ACBus=Dict()

#101,[-37.327065809684214, -37.353785256, -37.9549728]
bus_df=DataFrames.DataFrame(;
t =Int64[], bus=[], area = [], ThermalStandard=[], PowerLoad=[], ActivePowerBalance__ACBus=[],
StateEstimationInjections__ACBus=[], ptdf=[], flowcontribution=[], loopflowcontribution=[]
)
for t in 1:NT
for b in PSY.get_components(PSY.Bus,sys)
bn=get_number(b); bus_ActivePowerVariable__ThermalStandard[t,bn]=0
bus_ActivePowerTimeSeriesParameter__PowerLoad[t,bn]=0; bus_ActivePowerBalance__ACBus[t,bn]=0
bus_StateEstimationInjections__ACBus[t,bn]=0
end
end

for i in names(ab)
try b=parse(Int,i)
for t=1:NT bus_ActivePowerBalance__ACBus[t,b]=ab[t,i] end
catch e println("error= ",e,",i=",i) end
end
if se==1
for i in names(al)
try b=parse(Int,i)
for t=1:NT bus_StateEstimationInjections__ACBus[t,b]=al[t,i] end
catch e println("error= ",e,",i=",i) end
end
end

for i in names(ald)
try b=get_number(get_bus(get_component(PowerLoad, sys, i)))
for t=1:NT bus_ActivePowerTimeSeriesParameter__PowerLoad[t,b]=bus_ActivePowerTimeSeriesParameter__PowerLoad[t,b]+ald[t,i] end
catch e println("error= ",e,",i=",i) end
end

for i in names(ac)
try b=get_number(get_bus(get_component(ThermalStandard, sys, i)))
for t=1:NT bus_ActivePowerVariable__ThermalStandard[t,b]=bus_ActivePowerVariable__ThermalStandard[t,b]+ac[t,i] end
catch e println("error ",e,",i=",i) end
end

for t in 1:NT
for b in PSY.get_components(PSY.Bus,sys)
bn=get_number(b); area=get_name(get_area(b));
push!(bus_df,(t,bn,area, bus_ActivePowerVariable__ThermalStandard[t,bn], bus_ActivePowerTimeSeriesParameter__PowerLoad[t,bn],
bus_ActivePowerBalance__ACBus[t,bn], bus_StateEstimationInjections__ACBus[t,bn], ptdf[line,bn],
bus_ActivePowerBalance__ACBus[t,bn]*ptdf[line,bn],bus_StateEstimationInjections__ACBus[t,bn]*ptdf[line,bn]))
end
end

bus_df[!,"buscheck"]=bus_df[!,"ActivePowerBalance__ACBus"]-bus_df[!,"ThermalStandard"]/100-bus_df[!,"PowerLoad"]/100
println("bus check ActivePowerBalance__ACBus<>ThermalStandard/100-PowerLoad/100 ,",filter([:t,:buscheck] => (t,buscheck)-> (t>=1) && (abs(buscheck)>0.000001), bus_df))
#println(filter([:t,:bus] => (t,bus)-> (t>=1) && (bus==203), bus_df))
return(bus_df)
end

function check_models(uc0,uc2a)
sol0=Dict()
v0=uc0.variables[PowerSimulations.VariableKey{ActivePowerVariable, ThermalStandard}("")]
sol0=Dict(zip(name.(v0),value.(v0)))
for v in v0
fix(v,value.(v),force=true)
end

v01=uc0.variables[InfrastructureSystems.Optimization.VariableKey{FlowActivePowerVariable, AreaInterchange}("")]
for v in v01
fix(v,value.(v),force=true)
end

v2a=uc2a.variables[PowerSimulations.VariableKey{ActivePowerVariable, ThermalStandard}("")]
sol2a=Dict(zip(name.(v2a),value.(v2a)))

for (k,v) in sol2a
if abs(sol2a[k]-sol0[k])>0.00001
println("k,",k,",",sol2a[k],",",sol0[k])
var0=variable_by_name(uc0.JuMPmodel,k)
fix(var0,sol2a[k],force=true)
end
end
solve!(models.decision_models[1])
end

function objterms(uc0,ff)
objd0=Dict(); sol0=Dict(); objv0=0
obj0=objective_function(uc0.JuMPmodel); objd0=Dict(zip(name.(keys(obj0.terms)),values(obj0.terms)));
av0 = JuMP.all_variables(uc0.JuMPmodel); sol0=Dict(zip(name.(av0), value.(av0)))

open(ff,"w") do io
redirect_stdout(io) do
for (k,v) in objd0
objv0=objv0+sol0[k]*objd0[k]
println("name;value;objcoef;obj_contribution;",k,";",sol0[k],";",objd0[k],";", sol0[k]*objd0[k])
end
end
end
return(objd0,sol0, objv0)
end
ptdf=PTDF(sys)

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"),
Expand Down Expand Up @@ -301,7 +192,7 @@ execute_status = execute!(sim; enable_progress_bar=true);
uc0=models.decision_models[1].internal.container
uc2=models.decision_models[2].internal.container

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

######### Model check ############
open("mod_uc0.txt","w") do io
Expand Down Expand Up @@ -335,8 +226,8 @@ results = SimulationResults(sim)
results_uc0 = get_decision_problem_results(results, "UC0")
results_uc1 = get_decision_problem_results(results, "UC_Subsystem")

bus_df_uc0=read_bus_df(results_uc0,0)
bus_df_uc2=read_bus_df(results_uc1,0)
bus_df_uc0=read_bus_df(results_uc0,"A28",0)
bus_df_uc2=read_bus_df(results_uc1,"A28",0)

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))
Expand Down

0 comments on commit 18035e6

Please sign in to comment.