danielolsen commented 2 years ago

:beetle:

[x] I have checked that this issue has not already been reported.

Bug summary

If all buses within a zone have a Pd of zero, model building will crash within the _add_constraint_power_balance! function, due to a NaN value within bus_demand. The origin of this is NaN value is the following line: bus_share = bus_df[:, :load] ./ bus_df_with_zone_load[:, :load_sum] which divides by zero for at least one zone since the sum of all Pd values for the zone is zero.

https://github.com/Breakthrough-Energy/REISE.jl/blob/420dc059a785dc3db3781734ec0e264f31f04a28/src/model.jl#L55

Code for reproduction

Can be reproduced using the current HIFLD grid model for the Eastern interconnect, which has all zeros for Pd values for buses within the "Colorado Eastern" zone (zone_id = 6). To reproduce, you can call REISE.run_scenario_gurobi or REISE.run_scenario. Checking the data with python we see the root of the data problem.

>>> from powersimdata import Grid
>>> eastern = Grid("Eastern", "hifld")
>>> eastern.bus.groupby("zone_id")["Pd"].sum()
zone_id
1     10011.178932
2      6151.994940
6         0.000000
7      7185.898740
8      1924.068480
9     42205.325171
10    21388.176215
11     6307.999616
13    25755.139020
14    13525.860755
15     5916.248707
16     9188.076334
17     9426.187824
18    13712.366730
19    12097.884480
20     2684.338920
21    20058.991980
22    11412.078189
23    12401.671927
24     6132.018522
25      107.257620
27    20757.000844
28     1584.745305
29     3870.868983
31     2709.729240
32    17958.609315
33      555.734850
36    39349.735830
37    23618.239177
38     7985.510489
40    25689.231120
41     2125.034310
42    10309.792483
43     1495.255011
45    13807.662019
46     5488.202250
50    17033.393106
51     1254.869130
53    11788.683081
54     3645.733980
Name: Pd, dtype: float64

Actual outcome

The output produced by the above code, which may be a screenshot, console output, etc.

Expression contains an invalid NaN constant. This could be produced by `Inf - Inf`.

Stacktrace:
 [1] error(::String) at .\error.jl:33
 [2] _assert_isfinite(::JuMP.GenericAffExpr{Float64,JuMP.VariableRef}) at C:\Users\DanielOlsen\.julia\packages\JuMP\Xrr7O\src\aff_expr.jl:476
 [3] MathOptInterface.ScalarAffineFunction(::JuMP.GenericAffExpr{Float64,JuMP.VariableRef}) at C:\Users\DanielOlsen\.julia\packages\JuMP\Xrr7O\src\aff_expr.jl:507
 [4] moi_function at C:\Users\DanielOlsen\.julia\packages\JuMP\Xrr7O\src\aff_expr.jl:550 [inlined]
 [5] moi_function at C:\Users\DanielOlsen\.julia\packages\JuMP\Xrr7O\src\constraints.jl:396 [inlined]
 [6] add_constraint(::JuMP.Model, ::JuMP.ScalarConstraint{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},MathOptInterface.EqualTo{Float64}}, ::String) at C:\Users\DanielOlsen\.julia\packages\JuMP\Xrr7O\src\constraints.jl:547
 [7] add_constraint at C:\Users\DanielOlsen\.julia\packages\JuMP\Xrr7O\src\constraints.jl:546 [inlined]
 [8] _broadcast_getindex_evalf at .\broadcast.jl:648 [inlined]
 [9] _broadcast_getindex at .\broadcast.jl:621 [inlined]
 [10] getindex at .\broadcast.jl:575 [inlined]
 [11] copyto_nonleaf!(::Array{JuMP.ConstraintRef{JuMP.Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}},JuMP.ScalarShape},2}, ::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Tuple{Base.OneTo{Int64},Base.OneTo{Int64}},typeof(JuMP.add_constraint),Tuple{Base.RefValue{JuMP.Model},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(JuMP.build_constraint),Tuple{Base.RefValue{JuMP.var"#_error#92"{Tuple{Symbol,Symbol,Expr},Symbol,LineNumberNode}},Base.Broadcast.Extruded{Array{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},2},Tuple{Bool,Bool},Tuple{Int64,Int64}},Base.RefValue{MathOptInterface.EqualTo{Float64}}}}}}, ::CartesianIndices{2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}, ::CartesianIndex{2}, ::Int64) at .\broadcast.jl:1026
 [12] copy at .\broadcast.jl:880 [inlined]
 [13] materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(JuMP.add_constraint),Tuple{Base.RefValue{JuMP.Model},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(JuMP.build_constraint),Tuple{Base.RefValue{JuMP.var"#_error#92"{Tuple{Symbol,Symbol,Expr},Symbol,LineNumberNode}},Array{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},2},Base.RefValue{MathOptInterface.EqualTo{Float64}}}}}}) at .\broadcast.jl:837
 [14] macro expansion at C:\Users\DanielOlsen\.julia\packages\JuMP\Xrr7O\src\macros.jl:677 [inlined]
 [15] _add_constraint_power_balance!(::JuMP.Model, ::REISE.Case, ::REISE.Sets, ::REISE.Storage, ::REISE.DemandFlexibility, ::Array{Float64,2}, ::Bool) at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\model.jl:278
 [16] _build_model(::JuMP.Model; case::REISE.Case, storage::REISE.Storage, demand_flexibility::REISE.DemandFlexibility, start_index::Int64, interval_length::Int64, demand_scaling::Float64, load_shed_enabled::Bool, load_shed_penalty::Int64, trans_viol_enabled::Bool, trans_viol_penalty::Int64, initial_ramp_enabled::Bool, initial_ramp_g0::Array{Float64,1}, storage_e0::Array{Float64,1}, init_shifted_demand::Array{Float64,1}) at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\model.jl:728
 [17] interval_loop(::Gurobi.Env, ::Dict{String,Any}, ::Dict{String,Int64}, ::Int64, ::Int64, ::Int64, ::String, ::String) at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\loop.jl:73
 [18] (::REISE.var"#112#113"{Int64,Int64,Int64,String,Gurobi.Env})() at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\REISE.jl:91
 [19] #100 at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\save.jl:120 [inlined]
 [20] redirect_stderr(::REISE.var"#100#104"{REISE.var"#112#113"{Int64,Int64,Int64,String,Gurobi.Env}}, ::IOStream) at .\stream.jl:1150
 [21] #99 at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\save.jl:119 [inlined]
 [22] redirect_stdout(::REISE.var"#99#103"{IOStream,REISE.var"#112#113"{Int64,Int64,Int64,String,Gurobi.Env}}, ::IOStream) at .\stream.jl:1150
 [23] #98 at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\save.jl:118 [inlined]
 [24] open(::REISE.var"#98#102"{IOStream,REISE.var"#112#113"{Int64,Int64,Int64,String,Gurobi.Env}}, ::String, ::Vararg{String,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at .\io.jl:325
 [25] open at .\io.jl:323 [inlined]
 [26] #97 at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\save.jl:117 [inlined]
 [27] open(::REISE.var"#97#101"{REISE.var"#112#113"{Int64,Int64,Int64,String,Gurobi.Env},String}, ::String, ::Vararg{String,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at .\io.jl:325
 [28] open at .\io.jl:323 [inlined]
 [29] redirect_stdout_stderr(::REISE.var"#112#113"{Int64,Int64,Int64,String,Gurobi.Env}, ::String, ::String) at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\save.jl:116
 [30] run_scenario(; num_segments::Int64, interval::Int64, n_interval::Int64, start_index::Int64, inputfolder::String, outputfolder::String, threads::Nothing, optimizer_factory::Gurobi.Env, solver_kwargs::Dict{String,Int64}, model_kwargs::Dict{String,Integer}) at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\REISE.jl:89
 [31] run_scenario_gurobi(; solver_kwargs::Nothing, kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{6,Symbol},NamedTuple{(:interval, :n_interval, :start_index, :inputfolder, :outputfolder, :model_kwargs),Tuple{Int64,Int64,Int64,String,String,Dict{String,Integer}}}}) at C:\Users\DanielOlsen\repos\bes\REISE.jl\src\solver_specific\gurobi.jl:13

Expected outcome

In the case I was running into, I had set the demand to zero for that zone anyway, so I would have hoped that bus_demand would be all zeros for that zone. In this situation, I think we can safely replace any NaN values with zeros, and potentially warn the user/log of what's going on. In a situation for which there is non-zero demand for a zone with all-zero Pd values, I think we should detect this when loading the input files and halt before trying to build the model.

Environment

This bug should be environment agnostic.

Please specify your platform and versions of the relevant libraries you are using:

Operating system: Windows 10
REISE.jl revision (run git rev-parse origin/HEAD): e8639dd8a2fe42299e58e0c7ed49bc84527f5d07
Julia version: 1.5.4

lanesmith commented 2 years ago

In this situation, I think we can safely replace any NaN values with zeros, and potentially warn the user/log of what's going on. In a situation for which there is non-zero demand for a zone with all-zero Pd values, I think we should detect this when loading the input files and halt before trying to build the model.

I agree with both of these suggestions.