Merge pull request #586 from NatLabRockies/develop

v0.58.0 CHP load following

Author: Bill-Becker

CHANGELOG.md

@@ -25,6 +25,22 @@ Classify the change according to the following categories:
     ### Deprecated
     ### Removed
 
+## v0.58.0
+### Added
+- New optional attributes for **CHP** object **CHP.serve_absorption_chiller_only**, **CHP.months_serving_absorption_chiller_only**, and **CHP.follow_electrical_load**, which impose constraints on CHP operations if selected.  The default is set to `false` for both attributes.
+- New result **thermal_to_absorption_chiller_series_mmbtu_per_hour** added to heating technologies and new result **storage_to_absorption_chiller_series_mmbtu_per_hour** for hot thermal storage technologies.  This result is included as a part of the thermal site loads served, i.e., the adding this result does not change the existing results.
+- Added results fields to **HighTempThermalStorage** to match those of **HotThermalStorage**.
+
+### Changed
+- Updated heating dispatch results by separating heat flows to absorption chiller from heating load served (formerly, these were aggregated).
+
+### Fixed
+- Fixed a bug in which the CHP system requires a **DomesticHotWater** load.
+- Fixed a bug in which the storage to steam turbine flow was included in the thermal heating load served.
+
+### Changed
+- **HotThermalStorage** and **HighTempThermalStorage** output **storage_to_turbine_series_mmbtu_per_hour** to **storage_to_steamturbine_series_mmbtu_per_hour**
+
 ## v0.57.0
 ### Fixed
 - Include boiler emissions in emissions calculations

Project.toml

@@ -1,7 +1,7 @@
 name = "REopt"
 uuid = "d36ad4e8-d74a-4f7a-ace1-eaea049febf6"
 authors = ["Nick Laws", "Hallie Dunham <hallie.dunham@nlr.gov>", "Bill Becker <william.becker@nlr.gov>", "Bhavesh Rathod <bhavesh.rathod@nlr.gov>", "Alex Zolan <alexander.zolan@nlr.gov>", "Amanda Farthing <amanda.farthing@nlr.gov>", "Xiangkun Li <xiangkun.li@nlr.gov>", "An Pham <an.pham@nlr.gov>", "Byron Pullutasig <byron.pullatasig@nlr.gov>"]
-version = "0.57.0"
+version = "0.58.0"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

src/constraints/chp_constraints.jl

@@ -167,6 +167,68 @@ function add_chp_hourly_om_charges(m, p; _n="")
     nothing
 end
 
+"""
+    add_chp_to_absorption_chiller_only_constraints(m, p; _n="")
+
+Used in the function add_chp_constraints to add constraints that restrict heat output of CHP to only serve absorption chiller  
+load or send to waste in dispatch.
+"""
+function add_chp_to_absorption_chiller_only_constraints(m, p; _n="")
+    monthly_timesteps = get_monthly_time_steps(p.s.electric_load.year; time_steps_per_hour=p.s.settings.time_steps_per_hour)
+    for mth in p.s.chp.months_serving_absorption_chiller_only
+        @constraint(m, [t in p.techs.chp, q in [p.s.absorption_chiller.heating_load_input], ts in monthly_timesteps[mth]], 
+            m[Symbol("dvProductionToWaste"*_n)]["CHP",q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)]["CHP",q,ts] == m[Symbol("dvHeatingProduction"*_n)]["CHP",q,ts]
+        )
+        # During restricted months, CHP cannot serve other heating loads (space heating, DHW)
+        for q in setdiff(p.heating_loads,[p.s.absorption_chiller.heating_load_input])
+            for ts in monthly_timesteps[mth]
+                fix(m[Symbol("dvHeatToAbsorptionChiller"*_n)]["CHP",q,ts], 0.0, force=true)
+                fix(m[Symbol("dvHeatingProduction"*_n)]["CHP",q,ts], 0.0, force=true)
+                fix(m[Symbol("dvProductionToWaste"*_n)]["CHP",q,ts], 0.0, force=true)
+            end
+        end
+    end
+end
+
+
+"""
+    add_chp_electrical_load_following_constraints(m, p; _n="")
+
+Used in function add_chp_constraints to add constraints that restrict output of CHP to serve full electrical load or 
+run at capacity otherwise in dispatch (without regard to heat flows).
+"""
+function add_chp_electrical_load_following_constraints(m, p; _n="")
+    dv = "binCHPSizeExceedsElectricLoad"*_n
+    m[Symbol(dv)] = @variable(m, [p.time_steps], binary=true, base_name=dv)
+    dv = "dvCHPSizeTimesExcess"*_n
+    m[Symbol(dv)] = @variable(m, [p.time_steps], lower_bound=0, base_name=dv)
+    # binary variable enforcement for size >= load
+    max_diff_size_bigM = 2*max(p.max_sizes["CHP"], maximum(p.s.electric_load.loads_kw) #+ sum(p.heating_loads_kw[q][ts] for q in p.heating_loads))  #exclude heating electrification but include elec cooling? 
+    )
+    @constraint(m, [ts in p.time_steps],
+        m[Symbol("binCHPSizeExceedsElectricLoad"*_n)][ts] >= (m[Symbol("dvSize"*_n)]["CHP"] - p.s.electric_load.loads_kw[ts]) / max_diff_size_bigM
+    )
+    @constraint(m, [ts in p.time_steps],
+        m[Symbol("binCHPSizeExceedsElectricLoad"*_n)][ts] <= 1 - (p.s.electric_load.loads_kw[ts] - m[Symbol("dvSize"*_n)]["CHP"]) / max_diff_size_bigM
+    )
+    # set dvCHPSizeTimesExcess = binCHPSizeExceedsElectricLoad * dvSize
+    # big-M is min CF times heat load
+    
+    @constraint(m, [ts in p.time_steps],
+        m[Symbol("dvCHPSizeTimesExcess"*_n)][ts] >= p.production_factor["CHP",ts]*m[Symbol("dvSize"*_n)]["CHP"] - max_diff_size_bigM * (1-m[Symbol("binCHPSizeExceedsElectricLoad"*_n)][ts])  
+    )
+    @constraint(m, [ts in p.time_steps],
+        m[Symbol("dvCHPSizeTimesExcess"*_n)][ts] <= p.production_factor["CHP",ts]*m[Symbol("dvSize"*_n)]["CHP"]
+    )
+    @constraint(m, [ts in p.time_steps],
+        m[Symbol("dvCHPSizeTimesExcess"*_n)][ts] <= max_diff_size_bigM * m[Symbol("binCHPSizeExceedsElectricLoad"*_n)][ts]
+    )
+    #Enforce dispatch: output = system size - (overage)
+    @constraint(m, [ts in p.time_steps],
+        m[Symbol("dvRatedProduction"*_n)]["CHP",ts] >= p.production_factor["CHP",ts]*m[Symbol("dvSize"*_n)]["CHP"] - m[Symbol("dvCHPSizeTimesExcess"*_n)][ts] + p.s.electric_load.loads_kw[ts] * m[Symbol("binCHPSizeExceedsElectricLoad"*_n)][ts]
+    )
+end
+
 
 """
     add_chp_constraints(m, p; _n="")
@@ -208,4 +270,13 @@ function add_chp_constraints(m, p; _n="")
             end
         end
     end
+
+    if !isempty(p.techs.absorption_chiller) && p.s.chp.serve_absorption_chiller_only
+        add_chp_to_absorption_chiller_only_constraints(m, p; _n=_n)
+    end
+
+    if p.s.chp.follow_electrical_load
+        add_chp_electrical_load_following_constraints(m, p; _n=_n)
+    end
 end
+

src/constraints/load_balance.jl

@@ -125,9 +125,10 @@ function add_thermal_load_constraints(m, p; _n="")
                     p.heating_loads_kw[q][ts]
                     + sum(m[Symbol("dvProductionToWaste"*_n)][t,q,ts] for t in union(p.techs.heating, p.techs.chp))
                     + sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot, t in union(p.techs.heating, p.techs.chp))
-                    + sum(m[Symbol("dvCoolingProduction"*_n)][t,ts] / p.thermal_cop[t] for t in p.absorption_chillers_using_heating_load[q])
+                    + sum(m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts] for t in union(p.techs.heating, p.techs.chp))
                     + sum(m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] for t in p.techs.can_supply_steam_turbine)
                     + sum(m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts] for b in p.s.storage.types.hot)
+                    + sum(m[Symbol("dvHeatFromStorageToAbsorptionChiller"*_n)][b,q,ts] for b in p.s.storage.types.hot)
                 )
             else
                 @constraint(m, HeatLoadBalanceCon[q in p.heating_loads, ts in p.time_steps_with_grid],
@@ -137,14 +138,14 @@ function add_thermal_load_constraints(m, p; _n="")
                     p.heating_loads_kw[q][ts]
                     + sum(m[Symbol("dvProductionToWaste"*_n)][t,q,ts] for t in union(p.techs.heating, p.techs.chp))
                     + sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot, t in union(p.techs.heating, p.techs.chp))
-                    + sum(m[Symbol("dvCoolingProduction"*_n)][t,ts] / p.thermal_cop[t] for t in p.absorption_chillers_using_heating_load[q])
+                    + sum(m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts] for t in union(p.techs.heating, p.techs.chp))
+                    + sum(m[Symbol("dvHeatFromStorageToAbsorptionChiller"*_n)][b,q,ts] for b in p.s.storage.types.hot)
                 )
             end
 
         end
 
         if !isempty(p.techs.cooling)
-            
             ##Constraint (5a): Cold thermal loads
             @constraint(m, [ts in p.time_steps_with_grid],
                 sum(m[Symbol("dvCoolingProduction"*_n)][t,ts] for t in p.techs.cooling)
@@ -155,4 +156,39 @@ function add_thermal_load_constraints(m, p; _n="")
             )
         end
     end
+end
+
+function add_absorption_chiller_load_constraints(m,p;_n="")
+    if !isempty(p.techs.absorption_chiller)
+        for q in p.heating_loads
+            if !isempty(p.absorption_chillers_using_heating_load[q])
+                @constraint(m, abschlcon[ts in p.time_steps_with_grid],
+                    sum(m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts] for t in union(p.techs.heating, p.techs.chp))
+                    == 
+                    sum(m[Symbol("dvCoolingProduction"*_n)][t,ts] / p.thermal_cop[t] for t in p.absorption_chillers_using_heating_load[q])
+                )
+            else
+                for t in union(p.techs.heating, p.techs.chp)
+                    for ts in p.time_steps
+                        fix(m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts], 0.0, force=true)
+                    end
+                end
+            end
+        end
+        for t in intersect(p.techs.heating,["GHP","ASHPWaterHeater","ASHPSpaceHeater"])
+            for q in p.heating_loads
+                for ts in p.time_steps
+                    fix(m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts], 0.0, force=true)
+                end
+            end
+        end
+    else
+        for t in union(p.techs.heating, p.techs.chp)
+            for q in p.heating_loads
+                for ts in p.time_steps
+                    fix(m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts], 0.0, force=true)
+                end
+            end
+        end
+    end
 end

src/constraints/steam_turbine_constraints.jl

@@ -5,12 +5,12 @@ function steam_turbine_thermal_input(m, p; _n="")
     # This constraint is already included in storage_constraints.jl if HotThermalStorage and SteamTurbine are considered that also includes dvProductionToStorage["HotThermalStorage"] in LHS
     if isempty(p.s.storage.types.hot)
         @constraint(m, SupplySteamTurbineProductionLimit[t in p.techs.can_supply_steam_turbine, q in p.heating_loads, ts in p.time_steps],
-                    m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] <=
+                    m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts] <=
                     m[Symbol("dvHeatingProduction"*_n)][t,q,ts]
         )
     else
         @constraint(m, SupplySteamTurbineProductionLimit[t in p.techs.can_supply_steam_turbine, q in p.heating_loads, ts in p.time_steps],
-                    m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] + sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot) + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] <=
+                    m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] + sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot) + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts] <=
                     m[Symbol("dvHeatingProduction"*_n)][t,q,ts]
         )
     end
@@ -19,7 +19,7 @@ end
 function steam_turbine_production_constraints(m, p; _n="")
     # Constraint Steam Turbine Thermal Production
     @constraint(m, SteamTurbineThermalProductionCon[t in p.techs.steam_turbine, ts in p.time_steps],
-                sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] for q in p.heating_loads) == p.s.steam_turbine.thermal_produced_to_thermal_consumed_ratio * (
+                sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts] for q in p.heating_loads) == p.s.steam_turbine.thermal_produced_to_thermal_consumed_ratio * (
                     sum(m[Symbol("dvThermalToSteamTurbine"*_n)][tst,q,ts] for tst in p.techs.can_supply_steam_turbine, q in p.heating_loads) + 
                     sum(m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts] for b in p.s.storage.types.hot, q in p.heating_loads)
                     )

src/constraints/storage_constraints.jl

@@ -152,7 +152,7 @@ function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")
 			m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] <=  m[Symbol("dvHeatingProduction"*_n)][t,q,ts]
 			)
         @constraint(m, [q in p.heating_loads, ts in p.time_steps],
-            m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts] <= m[Symbol("dvHeatFromStorage"*_n)][b,q,ts]
+            m[Symbol("dvHeatFromStorageToTurbine"*_n)][b,q,ts] + m[Symbol("dvHeatFromStorageToAbsorptionChiller"*_n)][b,q,ts] <= m[Symbol("dvHeatFromStorage"*_n)][b,q,ts]
         )
         if !p.s.storage.attr[b].can_supply_steam_turbine
             for q in p.heating_loads
@@ -170,6 +170,23 @@ function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")
                 end
             end
         end
+    else
+        if !isempty(p.techs.absorption_chiller)
+            @constraint(m, [q in [p.s.absorption_chiller.heating_load_input], ts in p.time_steps],
+                m[Symbol("dvHeatFromStorageToAbsorptionChiller"*_n)][b,q,ts] <= m[Symbol("dvHeatFromStorage"*_n)][b,q,ts]
+            )
+            for q in setdiff(p.heating_loads,[p.s.absorption_chiller.heating_load_input])
+                for ts in p.time_steps
+                    fix(m[Symbol("dvHeatFromStorageToAbsorptionChiller"*_n)][b,q,ts], 0.0, force=true)
+                end
+            end
+        else
+            for q in p.heating_loads
+                for ts in p.time_steps
+                    fix(m[Symbol("dvHeatFromStorageToAbsorptionChiller"*_n)][b,q,ts], 0.0, force=true)
+                end
+            end
+        end
 	end
 
     # Constraint (4j)-1: Reconcile state-of-charge for (hot) thermal storage

src/constraints/thermal_tech_constraints.jl

@@ -28,18 +28,18 @@ function add_heating_tech_constraints(m, p; _n="")
     if !isempty(p.techs.steam_turbine)
         if !isempty(p.s.storage.types.hot)
             @constraint(m, [t in p.techs.can_supply_steam_turbine, q in p.heating_loads, ts in p.time_steps],
-                sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot) + m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] + m[Symbol("dvProductionToWaste"*_n)][t,q,ts]  <=
+                sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot) + m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts]  <=
                 m[Symbol("dvHeatingProduction"*_n)][t,q,ts]
             )
             if !isempty(setdiff(union(p.techs.heating, p.techs.chp),p.techs.can_supply_steam_turbine))
                 @constraint(m, [t in setdiff(union(p.techs.heating,p.techs.chp),p.techs.can_supply_steam_turbine), q in p.heating_loads, ts in p.time_steps],
-                    sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot) + m[Symbol("dvProductionToWaste"*_n)][t,q,ts]  <=  
+                    sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot) + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts]  <=  
                     m[Symbol("dvHeatingProduction"*_n)][t,q,ts]
                 )
             end
         else
             @constraint(m, [t in p.techs.can_supply_steam_turbine, q in p.heating_loads, ts in p.time_steps],
-                m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] + m[Symbol("dvProductionToWaste"*_n)][t,q,ts]  <=
+                m[Symbol("dvThermalToSteamTurbine"*_n)][t,q,ts] + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts]  <=
                 m[Symbol("dvHeatingProduction"*_n)][t,q,ts]
             )
             if !isempty(setdiff(union(p.techs.heating, p.techs.chp),p.techs.can_supply_steam_turbine))
@@ -51,7 +51,7 @@ function add_heating_tech_constraints(m, p; _n="")
     else
         if !isempty(p.s.storage.types.hot)
             @constraint(m, [t in union(p.techs.heating, p.techs.chp), q in p.heating_loads, ts in p.time_steps],
-                sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot) + m[Symbol("dvProductionToWaste"*_n)][t,q,ts]  <=
+                sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for b in p.s.storage.types.hot) + m[Symbol("dvProductionToWaste"*_n)][t,q,ts] + m[Symbol("dvHeatToAbsorptionChiller"*_n)][t,q,ts]  <=
                 m[Symbol("dvHeatingProduction"*_n)][t,q,ts]
             )
         else