add argument b to function add_hot_thermal_storage_dispatch_constraints

src/constraints/storage_constraints.jl

@@ -139,49 +139,45 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
     end
 end
 
-function add_hot_thermal_storage_dispatch_constraints(m, p; _n="")
+function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")
 
     # Constraint (4f)-1b: SteamTurbineTechs
 	if !isempty(p.techs.steam_turbine)
-		@constraint(m, SteamTurbineTechProductionFlowCon[b in p.s.storage.types.hot, t in p.techs.steam_turbine, q in p.heating_loads, ts in p.time_steps],
+		@constraint(m, SteamTurbineTechProductionFlowCon[t in p.techs.steam_turbine, q in p.heating_loads, ts in p.time_steps],
 			m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] <=  m[Symbol("dvHeatingProduction"*_n)][t,q,ts]
 			)
-        @constraint(m, StorageToTurbineProductionFlowCon[b in p.s.storage.types.hot, q in p.heating_loads, ts in p.time_steps],
+        @constraint(m, StorageToTurbineProductionFlowCon[q in p.heating_loads, ts in p.time_steps],
             m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts] <= m[Symbol("dvHeatFromStorage"*_n)][b,q,ts]
         )
-        for b in p.s.storage.types.hot
-            if !p.s.storage.attr[b].can_supply_steam_turbine
-                for q in p.heating_loads
-                    for ts in p.time_steps
-                        fix(m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts], 0.0, force=true)
-                    end
+        if !p.s.storage.attr[b].can_supply_steam_turbine
+            for q in p.heating_loads
+                for ts in p.time_steps
+                    fix(m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts], 0.0, force=true)
                 end
-            elseif p.s.storage.attr[b].supply_turbine_only
-                @constraint(m, StorageOnlyToSteamTurbineCon[q in p.heating_loads, ts in p.time_steps],
-                    m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] ==  m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts]
-                    )
-                for t in p.techs.steam_turbine
-                    for ts in p.time_steps, q in p.heating_loads
-                        fix(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts], 0.0, force=true)
-                    end
+            end
+        elseif p.s.storage.attr[b].supply_turbine_only
+            @constraint(m, StorageOnlyToSteamTurbineCon[q in p.heating_loads, ts in p.time_steps],
+                m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] ==  m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts]
+                )
+            for t in p.techs.steam_turbine
+                for ts in p.time_steps, q in p.heating_loads
+                    fix(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts], 0.0, force=true)
                 end
             end
         end
 	end
 
     # Constraint (4j)-1: Reconcile state-of-charge for (hot) thermal storage
-    if !isempty(setdiff(p.s.storage.types.hot, ["HotSensibleTes"]))
+    if b != "HotSensibleTes"
         @constraint(m, [b in setdiff(p.s.storage.types.hot, ["HotSensibleTes"]), ts in p.time_steps],
             m[Symbol("dvStoredEnergy"*_n)][b,ts] == m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + p.hours_per_time_step * (
                 p.s.storage.attr[b].charge_efficiency * sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for t in union(p.techs.heating, p.techs.chp), q in p.heating_loads) -
                 sum(m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] for q in p.heating_loads) / p.s.storage.attr[b].discharge_efficiency -
                 p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
             )
         )
-    end
-
-    if "HotSensibleTes" in p.s.storage.types.hot
-        @constraint(m, [b in intersect(p.s.storage.types.hot, ["HotSensibleTes"]), ts in p.time_steps],
+    else
+        @constraint(m, [ts in p.time_steps],
         m[Symbol("dvStoredEnergy"*_n)][b,ts] == m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + (1/p.s.settings.time_steps_per_hour) * (
             p.s.storage.attr[b].charge_efficiency * sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for t in union(p.techs.heating, p.techs.chp), q in p.heating_loads) -
             sum(m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] for q in p.heating_loads) / p.s.storage.attr[b].discharge_efficiency -
@@ -206,10 +202,10 @@ function add_hot_thermal_storage_dispatch_constraints(m, p; _n="")
         for q in p.heating_loads)
     )
 
-    if "HotSensibleTes" in p.s.storage.types.hot
+    if b == "HotSensibleTes"
         @constraint(m, [ts in p.time_steps],
-            m[Symbol("dvStorageEnergy"*_n)]["HotSensibleTes"] / p.s.storage.attr["HotSensibleTes"].num_charge_hours >= 
-            sum(m[Symbol("dvHeatToStorage"*_n)]["HotSensibleTes",t,q,ts] for t in union(p.techs.heating, p.techs.chp), q in p.heating_loads)
+            m[Symbol("dvStorageEnergy"*_n)][b] / p.s.storage.attr[b].num_charge_hours >= 
+            sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for t in union(p.techs.heating, p.techs.chp), q in p.heating_loads)
         )
         @constraint(m, [ts in p.time_steps],
         m[Symbol("dvStorageEnergy"*_n)]["HotSensibleTes"] / p.s.storage.attr["HotSensibleTes"].num_discharge_hours  >= 
@@ -232,26 +228,26 @@ function add_hot_thermal_storage_dispatch_constraints(m, p; _n="")
         end
     end
 
-    if "HotSensibleTes" in p.s.storage.types.hot && p.s.storage.attr["HotSensibleTes"].one_direction_flow
+    if b =="HotSensibleTes" && p.s.storage.attr[b].one_direction_flow
         dv = "binStorageCharge"*_n
-        m[Symbol(dv)] = @variable(m, [p.s.storage.types.hot, p.time_steps], base_name=dv, binary=true)
+        m[Symbol(dv)] = @variable(m, [["HotSensibleTes"], p.time_steps], base_name=dv, binary=true)
         dv = "binStorageDischarge"*_n
-        m[Symbol(dv)] = @variable(m, [p.s.storage.types.hot, p.time_steps], base_name=dv, binary=true)
+        m[Symbol(dv)] = @variable(m, [["HotSensibleTes"], p.time_steps], base_name=dv, binary=true)
 
-        max_storage_power = min(p.s.storage.attr["HotSensibleTes"].max_kw, 
+        max_storage_power = min(p.s.storage.attr[b].max_kw, 
             100 * maximum(sum(p.heating_loads_kw[q][ts] for q in p.heating_loads) for ts in p.time_steps)    
         )
 
         @constraint(m, SensibleTesChargeMax[ts in p.time_steps],
-            sum(m[Symbol("dvHeatToStorage"*_n)]["HotSensibleTes",t,q,ts] for t in union(p.techs.heating, p.techs.chp), q in p.heating_loads) <=
-            max_storage_power * m[Symbol("binStorageCharge"*_n)]["HotSensibleTes",ts]
+            sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for t in union(p.techs.heating, p.techs.chp), q in p.heating_loads) <=
+            max_storage_power * m[Symbol("binStorageCharge"*_n)][b,ts]
         )
         @constraint(m, SensibleTesDischargeMax[ts in p.time_steps],
-            sum(m[Symbol("dvHeatFromStorage"*_n)]["HotSensibleTes",q,ts] for q in p.heating_loads) <=
-            max_storage_power * m[Symbol("binStorageDischarge"*_n)]["HotSensibleTes",ts]
+            sum(m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] for q in p.heating_loads) <=
+            max_storage_power * m[Symbol("binStorageDischarge"*_n)][b,ts]
         )
         @constraint(m, SensibleTesFlowDirection[ts in p.time_steps],
-            m[Symbol("binStorageDischarge"*_n)]["HotSensibleTes",ts] + m[Symbol("binStorageCharge"*_n)]["HotSensibleTes",ts] <= 1
+            m[Symbol("binStorageDischarge"*_n)][b,ts] + m[Symbol("binStorageCharge"*_n)][b,ts] <= 1
         )
     end
 

src/core/reopt.jl

@@ -242,15 +242,13 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 				add_elec_storage_dispatch_constraints(m, p, b)
 			elseif b in p.s.storage.types.cold
 				add_cold_thermal_storage_dispatch_constraints(m, p, b)
-			elseif !(b in p.s.storage.types.hot)
-				println(b)
+			elseif b in p.s.storage.types.hot
+				add_hot_thermal_storage_dispatch_constraints(m, p, b)
+			else
 				throw(@error("Invalid storage does not fall in a thermal or electrical set"))
 			end
 		end
 	end
-	if !isempty(p.s.storage.types.hot)
-		add_hot_thermal_storage_dispatch_constraints(m, p)
-	end
 
 	if any(max_kw->max_kw > 0, (p.s.storage.attr[b].max_kw for b in p.s.storage.types.elec))
 		add_storage_sum_grid_constraints(m, p)

src/core/reopt_multinode.jl

@@ -127,14 +127,14 @@ function build_reopt!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}})
 					add_elec_storage_dispatch_constraints(m, p, b; _n=_n)
 				elseif b in p.s.storage.types.cold
 					add_cold_thermal_storage_dispatch_constraints(m, p, b; _n=_n)
+				elseif b in p.s.storage.types.hot
+					add_hot_thermal_storage_dispatch_constraints(m, p, b; _n=_n)
+				else
+					throw(@error("Invalid storage does not fall in a thermal or electrical set"))
 				end
             end
         end
 
-		if !isempty(p.s.storage.types.hot)
-			add_hot_thermal_storage_dispatch_constraints(m, p; _n=_n)
-		end
-
         if any(max_kw->max_kw > 0, (p.s.storage.attr[b].max_kw for b in p.s.storage.types.elec))
             add_storage_sum_grid_constraints(m, p; _n=_n)
         end

src/mpc/model.jl

@@ -106,14 +106,13 @@ function build_mpc!(m::JuMP.AbstractModel, p::MPCInputs)
 				add_elec_storage_dispatch_constraints(m, p, b)
 			elseif b in p.s.storage.types.cold
 				add_cold_thermal_storage_dispatch_constraints(m, p, b)
-			elseif !(b in p.s.storage.types.hot)
+			elseif b in p.s.storage.types.hot
+				add_hot_thermal_storage_dispatch_constraints(m, p, b)
+			else
 				throw(@error("Invalid storage does not fall in a thermal or electrical set"))
 			end
 		end
 	end
-	if !isempty(p.s.storage.types.hot)
-		add_hot_thermal_storage_dispatch_constraints(m, p)
-	end
 
 	if any(size_kw->size_kw > 0, (p.s.storage.attr[b].size_kw for b in p.s.storage.types.all))
 		add_storage_sum_grid_constraints(m, p)

src/mpc/model_multinode.jl

@@ -73,14 +73,13 @@ function build_mpc!(m::JuMP.AbstractModel, ps::AbstractVector{MPCInputs})
 					add_elec_storage_dispatch_constraints(m, p, b; _n=_n)
 				elseif b in p.s.storage.types.cold
 					add_cold_thermal_storage_dispatch_constraints(m, p, b; _n=_n)
-				elseif !(b in p.s.storage.types.hot)
-					@error("Invalid storage does not fall in a thermal or electrical set")
+				elseif b in p.s.storage.types.hot
+					add_hot_thermal_storage_dispatch_constraints(m, p, b; _n=_n)
+				else
+					throw(@error("Invalid storage does not fall in a thermal or electrical set"))
 				end
 			end
 		end
-		if !isempty(p.s.storage.types.hot)
-			add_hot_thermal_storage_dispatch_constraints(m, p; _n=_n)
-		end
 
 		if any(size_kw->size_kw > 0, (p.s.storage.attr[b].size_kw for b in p.s.storage.types.elec))
 			add_storage_sum_grid_constraints(m, p; _n=_n)

test/testing1.jl

@@ -0,0 +1,48 @@
+using REopt
+using JuMP
+# using Cbc
+# using HiGHS
+using Xpress
+using JSON
+# using Plots
+using Test
+
+ENV["NREL_DEVELOPER_API_KEY"]="ogQAO0gClijQdYn7WOKeIS02zTUYLbwYJJczH9St"
+
+@testset "Electric Heater" begin
+    d = JSON.parsefile("./scenarios/electric_heater.json")
+    d["SpaceHeatingLoad"]["annual_mmbtu"] = 0.4 * 8760
+    d["DomesticHotWaterLoad"]["annual_mmbtu"] = 0.4 * 8760
+    d["ProcessHeatLoad"]["annual_mmbtu"] = 0.2 * 8760
+    s = Scenario(d)
+    p = REoptInputs(s)
+    m = Model(optimizer_with_attributes(Xpress.Optimizer))
+    results = run_reopt(m, p)
+    println(results["Messages"])
+
+    #first run: Boiler produces the required heat instead of the electric heater - electric heater should not be purchased
+    @test results["ElectricHeater"]["size_mmbtu_per_hour"] ≈ 0.0 atol=0.1
+    @test results["ElectricHeater"]["annual_thermal_production_mmbtu"] ≈ 0.0 atol=0.1
+    @test results["ElectricHeater"]["annual_electric_consumption_kwh"] ≈ 0.0 atol=0.1
+    @test results["ElectricUtility"]["annual_energy_supplied_kwh"] ≈ 87600.0 atol=0.1
+
+    d["ExistingBoiler"]["fuel_cost_per_mmbtu"] = 100
+    d["ElectricHeater"]["installed_cost_per_mmbtu_per_hour"] = 1.0
+    d["ElectricTariff"]["monthly_energy_rates"] = [0,0,0,0,0,0,0,0,0,0,0,0]
+    s = Scenario(d)
+    p = REoptInputs(s)
+    m = Model(optimizer_with_attributes(Xpress.Optimizer))
+    results = run_reopt(m, p)
+    println(results["Messages"])
+
+    annual_thermal_prod = 0.8 * 8760  #80% efficient boiler --> 0.8 MMBTU of heat load per hour
+    annual_electric_heater_consumption = annual_thermal_prod * REopt.KWH_PER_MMBTU  #1.0 COP
+    annual_energy_supplied = 87600 + annual_electric_heater_consumption
+
+    #Second run: ElectricHeater produces the required heat with free electricity
+    @test results["ElectricHeater"]["size_mmbtu_per_hour"] ≈ 0.8 atol=0.1
+    @test results["ElectricHeater"]["annual_thermal_production_mmbtu"] ≈ annual_thermal_prod rtol=1e-4
+    @test results["ElectricHeater"]["annual_electric_consumption_kwh"] ≈ annual_electric_heater_consumption rtol=1e-4
+    @test results["ElectricUtility"]["annual_energy_supplied_kwh"] ≈ annual_energy_supplied rtol=1e-4
+
+end