Skip to content

Commit

Permalink
Add missing Maps to the genesysmod function.
Browse files Browse the repository at this point in the history
removed remaining older formulation using model[:   ] inside different files
  • Loading branch information
dqpinel committed Oct 10, 2024
1 parent 606a1a0 commit c1d8b95
Show file tree
Hide file tree
Showing 5 changed files with 70 additions and 55 deletions.
4 changes: 2 additions & 2 deletions src/genesysmod.jl
Original file line number Diff line number Diff line change
Expand Up @@ -191,10 +191,10 @@ function genesysmod(;elmod_daystep, elmod_hourstep, solver, DNLPsolver, year=201
end

elseif termination_status(model) == MOI.OPTIMAL
VarPar = genesysmod_variable_parameter(model, Sets, Params)
VarPar = genesysmod_variable_parameter(model, Sets, Params, Vars)
if switch_processed_results == 1
GENeSYS_MOD.genesysmod_results(model, Sets, Params, VarPar, Vars, Switch,
Settings, elapsed,(switch_dispatch == 1 ? "dispatch" : ""))
Settings, Maps, elapsed,(switch_dispatch == 1 ? "dispatch" : ""))
# GENeSYS_MOD.genesysmod_results_old(model, Sets, Params, VarPar, Vars, Switch,
# Settings, elapsed,"dispatch")
end
Expand Down
51 changes: 35 additions & 16 deletions src/genesysmod_dec.jl
Original file line number Diff line number Diff line change
Expand Up @@ -60,15 +60,24 @@ function genesysmod_dec(model,Sets, Params,Switch, Maps)

############### Activity Variables #############

@variable(model, RateOfActivity[𝓨,𝓛,𝓣,𝓜,𝓡] >= 0)
RateOfActivity = def_daa(𝓨,𝓛,𝓣,𝓜,𝓡)
TotalAnnualTechnologyActivityByMode = def_daa(𝓨,𝓣,𝓜,𝓡)
ProductionByTechnologyAnnual = def_daa(𝓨,𝓣,𝓕,𝓡)
UseByTechnologyAnnual = def_daa(𝓨,𝓣,𝓕,𝓡)
for y 𝓨 for r 𝓡 for t 𝓣
for m Maps.Tech_MO[t]
for l 𝓛
RateOfActivity[y,l,t,m,r] = @variable(model, lower_bound = 0, base_name= "RateOfActivity[$y,$l,$t,$m,$r]")
end
TotalAnnualTechnologyActivityByMode[y,t,m,r] = @variable(model, lower_bound = 0, base_name= "TotalAnnualTechnologyActivityByMode[$y,$t,$m,$r]")
end
for f Maps.Tech_Fuel[t]
ProductionByTechnologyAnnual[y,t,f,r] = @variable(model, lower_bound = 0, base_name= "ProductionByTechnologyAnnual[$y,$t,$f,$r]")
UseByTechnologyAnnual[y,t,f,r] = @variable(model, lower_bound = 0, base_name= "UseByTechnologyAnnual[$y,$t,$f,$r]")
end
end end end
@variable(model, TotalTechnologyAnnualActivity[𝓨,𝓣,𝓡] >= 0)

@variable(model, TotalAnnualTechnologyActivityByMode[𝓨,𝓣,𝓜,𝓡] >= 0)

@variable(model, ProductionByTechnologyAnnual[𝓨,𝓣,𝓕,𝓡] >= 0)

@variable(model, UseByTechnologyAnnual[𝓨,𝓣,𝓕,𝓡] >= 0)


@variable(model, TotalActivityPerYear[𝓡,𝓛,𝓣,𝓨] >= 0)
@variable(model, CurtailedEnergyAnnual[𝓨,𝓕,𝓡] >= 0)
@variable(model, CurtailedCapacity[𝓡,𝓛,𝓣,𝓨] >= 0)
Expand Down Expand Up @@ -158,14 +167,24 @@ function genesysmod_dec(model,Sets, Params,Switch, Maps)


######### Trade #############

Import = @variable(model, Import[𝓨,𝓛,𝓕,𝓡,𝓡] >= 0, container=JuMP.Containers.DenseAxisArray)
Export = @variable(model, Export[𝓨,𝓛,𝓕,𝓡,𝓡] >= 0, container=JuMP.Containers.DenseAxisArray)

NewTradeCapacity = @variable(model, NewTradeCapacity[𝓨, 𝓕, 𝓡, 𝓡] >= 0, container=JuMP.Containers.DenseAxisArray)
TotalTradeCapacity = @variable(model, TotalTradeCapacity[𝓨, 𝓕, 𝓡, 𝓡] >= 0, container=JuMP.Containers.DenseAxisArray)
NewTradeCapacityCosts = @variable(model, NewTradeCapacityCosts[𝓨, 𝓕, 𝓡, 𝓡] >= 0, container=JuMP.Containers.DenseAxisArray)
DiscountedNewTradeCapacityCosts = @variable(model, DiscountedNewTradeCapacityCosts[𝓨, 𝓕, 𝓡, 𝓡] >= 0, container=JuMP.Containers.DenseAxisArray)
Import = def_daa(𝓨,𝓛,𝓕,𝓡,𝓡)
Export = def_daa(𝓨,𝓛,𝓕,𝓡,𝓡)
NewTradeCapacity = def_daa(𝓨,𝓕,𝓡,𝓡)
TotalTradeCapacity = def_daa(𝓨,𝓕,𝓡,𝓡)
NewTradeCapacityCosts = def_daa(𝓨,𝓕,𝓡,𝓡)
DiscountedNewTradeCapacityCosts = def_daa(𝓨,𝓕,𝓡,𝓡)
for y 𝓨 for f 𝓕 for r1 𝓡 for r2 𝓡
if Params.TradeRoute[r1,r2,f,y] != 0
for l 𝓛
Import[y,l,f,r1,r2] = @variable(model, lower_bound= 0, base_name="Import[$y,$l,$f,$r1,$r2]")
Export[y,l,f,r1,r2] = @variable(model, lower_bound= 0, base_name="Export[$y,$l,$f,$r1,$r2]")
end
NewTradeCapacity[y,f,r1,r2] = @variable(model, lower_bound= 0, base_name="NewTradeCapacity[$y,$f,$r1,$r2]")
TotalTradeCapacity[y,f,r1,r2] = @variable(model, lower_bound= 0, base_name="TotalTradeCapacity[$y,$f,$r1,$r2]")
NewTradeCapacityCosts[y,f,r1,r2] = @variable(model, lower_bound= 0, base_name="NewTradeCapacityCosts[$y,$f,$r1,$r2]")
DiscountedNewTradeCapacityCosts[y,f,r1,r2] = @variable(model, lower_bound= 0, base_name="DiscountedNewTradeCapacityCosts[$y,$f,$r1,$r2]")
end
end end end end

NetTrade = @variable(model, NetTrade[𝓨,𝓛,𝓕,𝓡], container=JuMP.Containers.DenseAxisArray)
NetTradeAnnual = @variable(model, NetTradeAnnual[𝓨,𝓕,𝓡], container=JuMP.Containers.DenseAxisArray)
Expand Down
16 changes: 8 additions & 8 deletions src/genesysmod_equ.jl
Original file line number Diff line number Diff line change
Expand Up @@ -315,7 +315,7 @@ function genesysmod_equ(model,Sets,Params, Vars,Emp_Sets,Settings,Switch, Maps)
end

for y 𝓨 for t 𝓣 for r 𝓡 for l 𝓛
@constraint(model, model[:TotalCapacityAnnual][y,t,r] >= model[:CurtailedCapacity][r,l,t,y], base_name="CA3c_CurtailedCapacity|$(r)|$(l)|$(t)|$(y)")
@constraint(model, Vars.TotalCapacityAnnual[y,t,r] >= Vars.CurtailedCapacity[r,l,t,y], base_name="CA3c_CurtailedCapacity|$(r)|$(l)|$(t)|$(y)")
end end end end
print("Cstr: Cap Adequacy A3 : ",Dates.now()-start,"\n")

Expand All @@ -342,11 +342,11 @@ function genesysmod_equ(model,Sets,Params, Vars,Emp_Sets,Settings,Switch, Maps)
for l 𝓛
@constraint(model, Vars.Import[y,l,f,r,rr] == Vars.Export[y,l,f,rr,r], base_name="EB1_TradeBalanceEachTS|$(y)|$(l)|$(f)|$(r)|$(rr)")
end
else
#= else
for l ∈ 𝓛
JuMP.fix(Vars.Import[y,l,f,r,rr], 0; force=true)
JuMP.fix(Vars.Export[y,l,f,rr,r], 0; force=true)
end
end =#
end
end

Expand Down Expand Up @@ -405,7 +405,7 @@ function genesysmod_equ(model,Sets,Params, Vars,Emp_Sets,Settings,Switch, Maps)
for i eachindex(𝓨) for r 𝓡 for rr 𝓡
if Params.TradeRoute[r,rr,"Power",𝓨[i]] > 0
for l 𝓛
@constraint(model, (model[:Import][𝓨[i],l,"Power",r,rr]) <= model[:TotalTradeCapacity][𝓨[i],"Power",rr,r]*Params.YearSplit[l,𝓨[i]]*31.536 , base_name="TrC1_TradeCapacityPowerLinesImport|$(𝓨[i])|$(l)_Power|$(r)|$(rr)")
@constraint(model, (Vars.Import[𝓨[i],l,"Power",r,rr]) <= Vars.TotalTradeCapacity[𝓨[i],"Power",rr,r]*Params.YearSplit[l,𝓨[i]]*31.536 , base_name="TrC1_TradeCapacityPowerLinesImport|$(𝓨[i])|$(l)_Power|$(r)|$(rr)")
end
for f 𝓕
if Params.TradeCapacityGrowthCosts[r,rr,f] != 0
Expand All @@ -415,7 +415,7 @@ function genesysmod_equ(model,Sets,Params, Vars,Emp_Sets,Settings,Switch, Maps)
end
end
for f 𝓕
if Params.TradeRoute[r,rr,f,𝓨[i]] == 0 || Params.TradeCapacityGrowthCosts[r,rr,f] == 0
if Params.TradeRoute[r,rr,f,𝓨[i]] != 0 && Params.TradeCapacityGrowthCosts[r,rr,f] == 0
JuMP.fix(Vars.DiscountedNewTradeCapacityCosts[𝓨[i],f,r,rr],0; force=true)
end
end
Expand Down Expand Up @@ -449,9 +449,9 @@ function genesysmod_equ(model,Sets,Params, Vars,Emp_Sets,Settings,Switch, Maps)
base_name="TrC5a_NewTradeCapacityLimitH2|$(𝓨[i])|H2|$(r)|$(rr)")
end
for f 𝓕
if Params.TradeRoute[r,rr,f,𝓨[i]] == 0
#= if Params.TradeRoute[r,rr,f,𝓨[i]] == 0
JuMP.fix(Vars.NewTradeCapacity[𝓨[i],f,r,rr],0; force=true)
end
end =#
if Params.TradeCapacityGrowthCosts[r,rr,f] > 0 && f != "Power"
@constraint(model, sum(Vars.Import[𝓨[i],l,f,rr,r] for l 𝓛) <= Vars.TotalTradeCapacity[𝓨[i],f,r,rr],
base_name="TrC7_TradeCapacityLimitNonPower$(𝓨[i])|$(f)|$(r)|$(rr)")
Expand All @@ -463,7 +463,7 @@ function genesysmod_equ(model,Sets,Params, Vars,Emp_Sets,Settings,Switch, Maps)
base_name="TrC6_SymmetricalTransmissionExpansion|$(𝓨[i])|$(r)|$(rr)")
end

if Params.TradeRoute[r,rr,"Power",𝓨[i]] == 0 || Params.GrowthRateTradeCapacity[r,rr,"Power",𝓨[i]] == 0
if Params.TradeRoute[r,rr,"Power",𝓨[i]] != 0 && Params.GrowthRateTradeCapacity[r,rr,"Power",𝓨[i]] == 0
JuMP.fix(Vars.NewTradeCapacity[𝓨[i],"Power",r,rr],0; force=true)
end

Expand Down
30 changes: 15 additions & 15 deletions src/genesysmod_results.jl
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
append!(output_energy_balance_annual, combine(groupby(df_tmp, [:Region, :Sector, :Technology, :Fuel, :Type, :Unit, :PathwayScenario, :Year]), :Value=> sum; renamecols=false))
end

df_imp= convert_jump_container_to_df(value.(model[:Import]);dim_names=[:Year, :Timeslice, :Fuel, :Region, :Region2])
df_imp= convert_jump_container_to_df(value.(Vars.Import);dim_names=[:Year, :Timeslice, :Fuel, :Region, :Region2])
df_tmp = combine(groupby(df_imp, [:Year, :Timeslice, :Fuel, :Region]), :Value => sum; renamecols=false)

dict_col_value = Dict(:Sector=>"Trade", :Type=>"Import", :Unit=>"PJ",
Expand All @@ -181,7 +181,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
append!(output_energy_balance_annual, combine(groupby(df_tmp, [:Region, :Sector, :Technology, :Fuel, :Type, :Unit, :PathwayScenario, :Year]), :Value=> sum; renamecols=false))
end

df_exp= convert_jump_container_to_df(-(1) * value.(model[:Export]);dim_names=[:Year, :Timeslice, :Fuel, :Region, :Region2])
df_exp= convert_jump_container_to_df(-(1) * value.(Vars.Export);dim_names=[:Year, :Timeslice, :Fuel, :Region, :Region2])
df_tmp = combine(groupby(df_exp, [:Year, :Timeslice, :Fuel, :Region]), :Value => sum; renamecols=false)
setindex!(dict_col_value, "Export", :Type)
merge_df(df_tmp, dict_col_value, output_energy_balance, colnames)
Expand Down Expand Up @@ -211,15 +211,15 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
df_peak_capacity[!,:Type] .= "PeakCapacity"
df_peak_capacity[!,:PathwayScenario] .= "$(Switch.emissionPathway)_$(Switch.emissionScenario)"

df_new_capacity = convert_jump_container_to_df((value.(model[:NewCapacity]));dim_names=[:Year, :Technology, :Region])
df_new_capacity = convert_jump_container_to_df((value.(Vars.NewCapacity));dim_names=[:Year, :Technology, :Region])
df_new_capacity[!,:Type] .= "NewCapacity"
df_new_capacity[!,:PathwayScenario] .= "$(Switch.emissionPathway)_$(Switch.emissionScenario)"

df_residual_capacity = convert_jump_container_to_df(Params.ResidualCapacity;dim_names=[:Region, :Technology, :Year])
df_residual_capacity[!,:Type] .= "ResidualCapacity"
df_residual_capacity[!,:PathwayScenario] .= "$(Switch.emissionPathway)_$(Switch.emissionScenario)"

df_total_capacity = convert_jump_container_to_df(value.(model[:TotalCapacityAnnual][:,tmp_techs,:]);dim_names=[:Year, :Technology, :Region])
df_total_capacity = convert_jump_container_to_df(value.(Vars.TotalCapacityAnnual[:,tmp_techs,:]);dim_names=[:Year, :Technology, :Region])
df_total_capacity[!,:Type] .= "TotalCapacity"
df_total_capacity[!,:PathwayScenario] .= "$(Switch.emissionPathway)_$(Switch.emissionScenario)"

Expand Down Expand Up @@ -249,7 +249,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
colnames = [:Region, :Sector, :Emission, :Technology, :Type, :PathwayScenario, :Year, :Value]
output_emissions = DataFrame([name => [] for name in colnames])

df_technology_emission = convert_jump_container_to_df(value.(model[:AnnualTechnologyEmission]);dim_names=[:Year, :Technology, :Emission, :Region])
df_technology_emission = convert_jump_container_to_df(value.(Vars.AnnualTechnologyEmission);dim_names=[:Year, :Technology, :Emission, :Region])
df_technology_emission[!,:Type] .= "Emissions"
df_technology_emission[!,:PathwayScenario] .= "$(Switch.emissionPathway)_$(Switch.emissionScenario)"

Expand Down Expand Up @@ -478,7 +478,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
colnames = [:Region, :Region2, :Type, :Year, :Value]
output_trade_capacity = DataFrame([name => [] for name in colnames])

df_tmp = convert_jump_container_to_df(value.(model[:TotalTradeCapacity][:,["Power"],:,:]);dim_names=[:Year, :Fuel, :Region, :Region2])
df_tmp = convert_jump_container_to_df(value.(Vars.TotalTradeCapacity[:,["Power"],:,:]);dim_names=[:Year, :Fuel, :Region, :Region2])
df_tmp[!,:Type] .= "Power Transmissions Capacity"
if !isempty(df_tmp)
select!(df_tmp,colnames)
Expand Down Expand Up @@ -513,7 +513,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
end
for se Sets.Sector
if sum(( Params.TagDemandFuelToSector[f,se]>0 ? Params.TagDemandFuelToSector[f,se]*VarPar.ProductionAnnual[y,f,r] : 0) for f Sets.Fuel for r Sets.Region_full ) > 0
ElectrificationRate[se,y] = sum(Params.TagDemandFuelToSector[f,se]*Params.TagElectricTechnology[t]*value(model[:ProductionByTechnologyAnnual][y,t,f,r]) for f Sets.Fuel for t Sets.Technology for r Sets.Region_full if value(model[:ProductionByTechnologyAnnual][y,t,f,r]) > 0)/sum(Params.TagDemandFuelToSector[f,se]*VarPar.ProductionAnnual[y,f,r] for f Sets.Fuel for r Sets.Region_full if Params.TagDemandFuelToSector[f,se]>0)
ElectrificationRate[se,y] = sum(Params.TagDemandFuelToSector[f,se]*Params.TagElectricTechnology[t]*value(Vars.ProductionByTechnologyAnnual[y,t,f,r]) for f Sets.Fuel for t Sets.Technology for r Sets.Region_full if value(Vars.ProductionByTechnologyAnnual[y,t,f,r]) > 0)/sum(Params.TagDemandFuelToSector[f,se]*VarPar.ProductionAnnual[y,f,r] for f Sets.Fuel for r Sets.Region_full if Params.TagDemandFuelToSector[f,se]>0)
end
end
end
Expand Down Expand Up @@ -566,7 +566,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
append!(output_other, df_tmp)
end

df_tmp = convert_jump_container_to_df(value.(model[:UseByTechnologyAnnual][:,:,:,:])/3.6;dim_names=[:Year, :Dim1, :Dim2, :Region])
df_tmp = convert_jump_container_to_df(value.(Vars.UseByTechnologyAnnual[:,:,:,:])/3.6;dim_names=[:Year, :Dim1, :Dim2, :Region])
df_tmp[!,:Type] .= "FinalEnergyConsumption"
if !isempty(df_tmp)
select!(df_tmp,colnames)
Expand Down Expand Up @@ -622,7 +622,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
fed= JuMP.Containers.DenseAxisArray(zeros(length(Sets.Year),length(Sets.Fuel),length(Sets.Region_full),length(Sets.Sector)), Sets.Year, Sets.Fuel, Sets.Region_full, Sets.Sector)
for se FinalDemandSector for y Sets.Year for f Sets.Fuel for r Sets.Region_full
if f ["Area_Rooftop_Residential","Area_Rooftop_Commercial","Heat_District"]
fed[y,f,r,se] = sum(value(model[:UseByTechnologyAnnual][y,t,f,r]) for t Sets.Technology if Params.TagTechnologyToSector[t,se] != 0)/3.6
fed[y,f,r,se] = sum(value(Vars.UseByTechnologyAnnual[y,t,f,r]) for t Sets.Technology if Params.TagTechnologyToSector[t,se] != 0)/3.6
end
end end end end

Expand Down Expand Up @@ -717,7 +717,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
pe_p= JuMP.Containers.DenseAxisArray(zeros(length(Sets.Year),length(Sets.Fuel),length(Sets.Region_full)), Sets.Year, Sets.Fuel, Sets.Region_full)
for y Sets.Year for f Sets.Fuel for r Sets.Region_full
if f ["Area_Rooftop_Residential","Area_Rooftop_Commercial"]
pe[y,f,r] = sum(value(model[:ProductionByTechnologyAnnual][y,t,f,r]) for t Sets.Technology if sum(Params.InputActivityRatio[r,t,:,:,y]) == 0)/3.6
pe[y,f,r] = sum(value(Vars.ProductionByTechnologyAnnual[y,t,f,r]) for t Sets.Technology if sum(Params.InputActivityRatio[r,t,:,:,y]) == 0)/3.6
pe_p[y,f,r] = pe[y,f,r] / sum(pe[y,:,r])
end
end end end
Expand Down Expand Up @@ -780,7 +780,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
eg_o_p= JuMP.Containers.DenseAxisArray(zeros(length(Sets.Year),length(Sets.Region_full)), Sets.Year, Sets.Region_full)

for y Sets.Year for r Sets.Region_full
div= sum(value(model[:ProductionByTechnologyAnnual][y,t,"Power",r]) for t Sets.Technology if Params.TagTechnologyToSector[t,"Storages"]==0)/3.6
div= sum(value(Vars.ProductionByTechnologyAnnual[y,t,"Power",r]) for t Sets.Technology if Params.TagTechnologyToSector[t,"Storages"]==0)/3.6

eg_temp = Dict{Tuple, Float64}()
for t in Sets.Technology
Expand All @@ -803,9 +803,9 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
eg_p[key...] = val / div
end

eg_s[y,r] = sum(value(model[:ProductionByTechnologyAnnual][y,t,"Power",r]) for t Params.TagTechnologyToSubsets["Solar"]) /3.6
eg_w[y,r] = sum(value(model[:ProductionByTechnologyAnnual][y,t,"Power",r]) for t Params.TagTechnologyToSubsets["Wind"]) /3.6
eg_h[y,r] = sum(value(model[:ProductionByTechnologyAnnual][y,t,"Power",r]) for t Params.TagTechnologyToSubsets["Hydro"]) /3.6
eg_s[y,r] = sum(value(Vars.ProductionByTechnologyAnnual[y,t,"Power",r]) for t Params.TagTechnologyToSubsets["Solar"]) /3.6
eg_w[y,r] = sum(value(Vars.ProductionByTechnologyAnnual[y,t,"Power",r]) for t Params.TagTechnologyToSubsets["Wind"]) /3.6
eg_h[y,r] = sum(value(Vars.ProductionByTechnologyAnnual[y,t,"Power",r]) for t Params.TagTechnologyToSubsets["Hydro"]) /3.6
eg_s_p[y,r] = eg_s[y,r]/div
eg_w_p[y,r] = eg_w[y,r]/div
eg_h_p[y,r] = eg_h[y,r]/div
Expand Down Expand Up @@ -894,7 +894,7 @@ function genesysmod_results(model,Sets, Params, VarPar, Vars, Switch, Settings,
for y Sets.Year for r Sets.Region_full
for f Sets.Fuel
if sum(Params.OutputActivityRatio[r,t,f,m,y] for t Params.TagTechnologyToSubsets["ImportTechnology"] for m Sets.Mode_of_operation) != 0
ispe[y,f,r] = (sum(value(model[:ProductionByTechnologyAnnual][y,t,f,r]) for t Params.TagTechnologyToSubsets["ImportTechnology"])/3.6)/sum(pe[y,:,r])
ispe[y,f,r] = (sum(value(Vars.ProductionByTechnologyAnnual[y,t,f,r]) for t Params.TagTechnologyToSubsets["ImportTechnology"])/3.6)/sum(pe[y,:,r])
end
end
end end
Expand Down
Loading

0 comments on commit c1d8b95

Please sign in to comment.