diff --git a/Multi-Stage/YC_areaptdf_DA_2stage.jl b/Multi-Stage/YC_areaptdf_DA_2stage.jl index 26fc130..5e3d134 100644 --- a/Multi-Stage/YC_areaptdf_DA_2stage.jl +++ b/Multi-Stage/YC_areaptdf_DA_2stage.jl @@ -24,6 +24,8 @@ using Dates using DataFrames using HiGHS + +include("YC_test_function.jl") #include("coordination_EnergyOnly.jl") #include("GlobalM2M.jl") #import StatsPlots @@ -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"), @@ -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 @@ -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))