Allow BESS to Export

CHANGELOG.md

@@ -25,6 +25,10 @@ Classify the change according to the following categories:
     ### Deprecated
     ### Removed
 
+## bess-export
+### Added 
+- Add decision variable **dvStorageToGrid** and **ElectricStorage** input options for grid export: `can_net_meter`, `can_wholesale`, `can_export_beyond_nem_limit`
+
 ## v0.59.0
 ### Added
 - Added two new size classes for **SteamTurbine** tech with new tech size ranges.

docs/src/developer/inputs.md

@@ -40,7 +40,7 @@ end
 The set maps are best explained with an example. The `techs_by_exportbin` map uses each technology'sattributes (eg. `PV`) to map each technology to which export bins that technology can access. The export bins include:
 1. `:NEM` (Net Energy Metering)
 2. `:WHL` (Wholesale)
-3. `:EXC` (Excess, beyond NEM))
+3. `:EXC` (Excess, beyond NEM)
 The bins that a technology can access are determined by the technologies attributes `can_net_meter`, `can_wholesale`, and `can_export_beyond_nem_limit`. So if `PV.can_net_meter = true`, `Wind.can_net_meter = true` and all the other attributes are `false` then the `techs_by_exportbin` will only have one non-empty key:
 ```julia
 techs_by_exportbin = Dict(

src/constraints/electric_utility_constraints.jl

@@ -22,31 +22,40 @@ function add_export_constraints(m, p; _n="")
     WHL_benefit = 0
     NEM_techs = String[t for t in p.techs.elec if :NEM in p.export_bins_by_tech[t]]
     WHL_techs = String[t for t in p.techs.elec if :WHL in p.export_bins_by_tech[t]]
+    NEM_storage = String[b for b in p.s.storage.types.elec if :NEM in p.export_bins_by_storage[b]]
+    WHL_storage = String[b for b in p.s.storage.types.elec if :WHL in p.export_bins_by_storage[b]]
 
-    if !isempty(NEM_techs)
+    if !isempty(vcat(NEM_techs, NEM_storage))
         # Constraint (9c): Net metering only -- can't sell more than you purchase
         # hours_per_time_step is cancelled on both sides, but used for unit consistency (convert power to energy)
         @constraint(m,
-            p.hours_per_time_step * sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] 
-            for t in NEM_techs, ts in p.time_steps)
+            p.hours_per_time_step * (
+                sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] for t in NEM_techs, ts in p.time_steps) + 
+                sum( m[Symbol("dvStorageToGrid"*_n)][b, :NEM, ts] for b in NEM_storage, ts in p.time_steps)
+            )
             <= p.hours_per_time_step * sum( m[Symbol("dvGridPurchase"*_n)][ts, tier]
                 for ts in p.time_steps, tier in 1:p.s.electric_tariff.n_energy_tiers)
         )
 
-        if p.s.electric_utility.net_metering_limit_kw == p.s.electric_utility.interconnection_limit_kw && isempty(WHL_techs)
+        if p.s.electric_utility.net_metering_limit_kw == p.s.electric_utility.interconnection_limit_kw && isempty(vcat(WHL_techs, WHL_storage))
             # no need for binNEM nor binWHL
             binNEM = 1
+            # Note: BESS can export but is not included in interconnection limit. 
             @constraint(m,
                 sum(m[Symbol("dvSize"*_n)][t] for t in NEM_techs) <= p.s.electric_utility.interconnection_limit_kw
             )
             NEM_benefit = @expression(m, p.pwf_e * p.hours_per_time_step *
-                sum( sum(p.s.electric_tariff.export_rates[:NEM][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] 
-                    for t in p.techs_by_exportbin[:NEM]) for ts in p.time_steps)
+                sum(p.s.electric_tariff.export_rates[:NEM][ts] * (
+                    sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] for t in p.techs_by_exportbin[:NEM]) +
+                    sum( m[Symbol("dvStorageToGrid"*_n)][b, :NEM, ts] for b in p.storage_by_exportbin[:NEM])
+                ) for ts in p.time_steps)
             )
             if :EXC in p.s.electric_tariff.export_bins
                 EXC_benefit = @expression(m, p.pwf_e * p.hours_per_time_step *
-                    sum( sum(p.s.electric_tariff.export_rates[:EXC][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :EXC, ts] 
-                        for t in p.techs_by_exportbin[:EXC]) for ts in p.time_steps)
+                    sum(p.s.electric_tariff.export_rates[:EXC][ts] * ( 
+                        sum( m[Symbol("dvProductionToGrid"*_n)][t, :EXC, ts] for t in p.techs_by_exportbin[:EXC]) +
+                        sum( m[Symbol("dvStorageToGrid"*_n)][b, :EXC, ts] for b in p.storage_by_exportbin[:EXC])
+                    ) for ts in p.time_steps)
                 )
             end
         else
@@ -63,6 +72,7 @@ function add_export_constraints(m, p; _n="")
             NEM_benefit = @variable(m, lower_bound = max_bene)
 
             # If choosing to take advantage of NEM, must have total capacity less than net_metering_limit_kw
+            # Note: BESS can export but is not included in net_metering_limit_kw or interconnection_limit_kw.
             if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name)
                 @constraint(m,
                     binNEM => {sum(m[Symbol("dvSize"*_n)][t] for t in NEM_techs) <= p.s.electric_utility.net_metering_limit_kw}
@@ -93,57 +103,81 @@ function add_export_constraints(m, p; _n="")
             if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name)
                 @constraint(m,
                     binNEM => {NEM_benefit >= p.pwf_e * p.hours_per_time_step *
-                        sum( sum(p.s.electric_tariff.export_rates[:NEM][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] 
-                            for t in p.techs_by_exportbin[:NEM]) for ts in p.time_steps)
+                        sum(p.s.electric_tariff.export_rates[:NEM][ts] * ( 
+                            sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] for t in p.techs_by_exportbin[:NEM]) +
+                            sum( m[Symbol("dvStorageToGrid"*_n)][b, :NEM, ts] for b in p.storage_by_exportbin[:NEM])
+                            ) for ts in p.time_steps
+                        )
                     }
                 )
                 @constraint(m, !binNEM => {NEM_benefit >= 0})
                 @constraint(m,[ts in p.time_steps_with_grid, t in p.techs_by_exportbin[:NEM]], 
                     !binNEM => { m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] == 0
                         }
                 )
+                @constraint(m,[ts in p.time_steps_with_grid, b in p.storage_by_exportbin[:NEM]], 
+                    !binNEM => { m[Symbol("dvStorageToGrid"*_n)][b, :NEM, ts] == 0
+                        }
+                )
             else
                 @constraint(m,
                     NEM_benefit >= p.pwf_e * p.hours_per_time_step *
-                        sum( sum(p.s.electric_tariff.export_rates[:NEM][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] 
-                            for t in p.techs_by_exportbin[:NEM]) for ts in p.time_steps)
+                        sum(p.s.electric_tariff.export_rates[:NEM][ts] * (
+                            sum( m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] for t in p.techs_by_exportbin[:NEM]) +
+                            sum( m[Symbol("dvStorageToGrid"*_n)][b, :NEM, ts] for b in p.storage_by_exportbin[:NEM])
+                            ) for ts in p.time_steps
+                        )
                 )
                 @constraint(m, NEM_benefit >= max_bene * binNEM)
                 @constraint(m,[ts in p.time_steps_with_grid, t in p.techs_by_exportbin[:NEM]], 
                     m[Symbol("dvProductionToGrid"*_n)][t, :NEM, ts] <= binNEM * sum(p.s.electric_load.loads_kw)
                 )
+                @constraint(m,[ts in p.time_steps_with_grid, b in p.storage_by_exportbin[:NEM]], 
+                    m[Symbol("dvStorageToGrid"*_n)][b, :NEM, ts] <= binNEM * sum(p.s.electric_load.loads_kw)
+                )
             end
 
             EXC_benefit = 0
             if :EXC in p.s.electric_tariff.export_bins
-                EXC_benefit = @variable(m, lower_bound = max_bene)
+                EXC_benefit = @variable(m, lower_bound = max_bene) #lower bound because the benefit is treated as a negative cost
                 if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name)
                     @constraint(m,
                         binNEM => {EXC_benefit >= p.pwf_e * p.hours_per_time_step *
-                            sum( sum(p.s.electric_tariff.export_rates[:EXC][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :EXC, ts] 
-                                for t in p.techs_by_exportbin[:EXC]) for ts in p.time_steps)
+                            sum( p.s.electric_tariff.export_rates[:EXC][ts] * ( 
+                                sum( m[Symbol("dvProductionToGrid"*_n)][t, :EXC, ts] for t in p.techs_by_exportbin[:EXC]) +
+                                sum( m[Symbol("dvStorageToGrid"*_n)][b, :EXC, ts] for b in p.storage_by_exportbin[:EXC])
+                                ) for ts in p.time_steps
+                            )
                         }
                     )
+
+
                     @constraint(m, !binNEM => {EXC_benefit >= 0})
                 else
                     @constraint(m,
                         EXC_benefit >= p.pwf_e * p.hours_per_time_step *
-                            sum( sum(p.s.electric_tariff.export_rates[:EXC][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :EXC, ts] 
-                                for t in p.techs_by_exportbin[:EXC]) for ts in p.time_steps)
+                            sum( p.s.electric_tariff.export_rates[:EXC][ts] * ( 
+                                sum( m[Symbol("dvProductionToGrid"*_n)][t, :EXC, ts] for t in p.techs_by_exportbin[:EXC]) +
+                                sum( m[Symbol("dvStorageToGrid"*_n)][b, :EXC, ts] for b in p.storage_by_exportbin[:EXC])
+                                ) for ts in p.time_steps
+                            )
                     )
                     @constraint(m, EXC_benefit >= max_bene * binNEM)
                 end
             end
         end
     end
 
-    if !isempty(WHL_techs)
+    if !isempty(vcat(WHL_techs, WHL_storage))
 
         if typeof(binNEM) <: Real  # no need for wholesale binary
             binWHL = 1
             WHL_benefit = @expression(m, p.pwf_e * p.hours_per_time_step *
-                sum( sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] 
-                        for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
+                sum( p.s.electric_tariff.export_rates[:WHL][ts] * (
+                    sum(m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] for t in p.techs_by_exportbin[:WHL]) +
+                    sum(m[Symbol("dvStorageToGrid"*_n)][b, :WHL, ts] for b in p.storage_by_exportbin[:WHL])
+                    ) for ts in p.time_steps
+                )
             )
         else
             binWHL = @variable(m, binary = true)
@@ -155,16 +189,22 @@ function add_export_constraints(m, p; _n="")
             if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name)
                 @constraint(m,
                     binWHL => {WHL_benefit >= p.pwf_e * p.hours_per_time_step *
-                        sum( sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] 
-                                for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
+                        sum(p.s.electric_tariff.export_rates[:WHL][ts] * (
+                            sum(m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] for t in p.techs_by_exportbin[:WHL]) + 
+                            sum(m[Symbol("dvStorageToGrid"*_n)][b, :WHL, ts] for b in p.storage_by_exportbin[:WHL])
+                            ) for ts in p.time_steps
+                        )
                     }
                 )
                 @constraint(m, !binWHL => {WHL_benefit >= 0})
             else
                 @constraint(m,
                     WHL_benefit >= p.pwf_e * p.hours_per_time_step *
-                        sum( sum(p.s.electric_tariff.export_rates[:WHL][ts] * m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] 
-                                for t in p.techs_by_exportbin[:WHL]) for ts in p.time_steps)
+                        sum(p.s.electric_tariff.export_rates[:WHL][ts] * (
+                                sum(m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] for t in p.techs_by_exportbin[:WHL]) + 
+                                sum(m[Symbol("dvStorageToGrid"*_n)][b, :WHL, ts] for b in p.storage_by_exportbin[:WHL])                            
+                            ) for ts in p.time_steps
+                        )
                 )
                 @constraint(m, WHL_benefit >= max_bene * binWHL)
             end
@@ -196,7 +236,6 @@ function add_monthly_peak_constraint(m, p; _n="")
             sum(p.production_factor[t, ts] * p.levelization_factor[t] * m[Symbol("dvRatedProduction"*_n)][t, ts] for t in p.techs.chp) - 
             sum(sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) for t in p.techs.chp) -
             sum(sum(m[Symbol("dvProductionToGrid")][t,u,ts] for u in p.export_bins_by_tech[t]) for t in p.techs.chp)
-
         )
     else
         @constraint(m, [mth in p.months, ts in p.s.electric_tariff.time_steps_monthly[mth]],
@@ -287,17 +326,20 @@ function add_simultaneous_export_import_constraint(m, p; _n="")
         )
         @constraint(m, ExportOnlyAfterSiteLoadMetCon[ts in p.time_steps],
             !m[Symbol("binNoGridPurchases"*_n)][ts] => {
-                sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t]) <= 0
+                sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t]) + 
+                sum(m[Symbol("dvStorageToGrid"*_n)][b, u, ts] for b in p.s.storage.types.elec, u in p.export_bins_by_storage[b]) <= 0
             }
         )
     else
-        bigM_hourly_load = maximum(p.s.electric_load.loads_kw)+maximum(p.s.space_heating_load.loads_kw)+maximum(p.s.process_heat_load.loads_kw)+maximum(p.s.dhw_load.loads_kw)+maximum(p.s.cooling_load.loads_kw_thermal)
+        bigM_hourly_load_plus_battery = maximum(p.s.electric_load.loads_kw)+maximum(p.s.space_heating_load.loads_kw)+maximum(p.s.process_heat_load.loads_kw)+maximum(p.s.dhw_load.loads_kw)+maximum(p.s.cooling_load.loads_kw_thermal)+sum(Real[p.s.storage.attr[b].max_kw for b in p.s.storage.types.elec])
         @constraint(m, NoGridPurchasesBinary[ts in p.time_steps],
             sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) +
-            sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec) <= bigM_hourly_load*(1-m[Symbol("binNoGridPurchases"*_n)][ts])
+            sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec) <= bigM_hourly_load_plus_battery*(1-m[Symbol("binNoGridPurchases"*_n)][ts])
         )
         @constraint(m, ExportOnlyAfterSiteLoadMetCon[ts in p.time_steps],
-            sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t]) <= bigM_hourly_load * m[Symbol("binNoGridPurchases"*_n)][ts]
+            sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t]) +
+            sum(m[Symbol("dvStorageToGrid"*_n)][b, u, ts] for b in p.s.storage.types.elec, u in p.export_bins_by_storage[b])
+            <= bigM_hourly_load_plus_battery * m[Symbol("binNoGridPurchases"*_n)][ts]
         )
     end
 end
@@ -391,7 +433,7 @@ end
 function add_elec_utility_expressions(m, p; _n="")
 
     if !isempty(p.s.electric_tariff.export_bins) && !isempty(p.techs.all)
-        # NOTE: levelization_factor is baked into dvProductionToGrid
+        # NOTE: levelization_factor is baked into dvProductionToGrid and dvStorageToGrid 
         m[Symbol("TotalExportBenefit"*_n)] = m[Symbol("NEM_benefit"*_n)] + m[Symbol("WHL_benefit"*_n)] +
                                              m[Symbol("EXC_benefit"*_n)]
     else

src/constraints/emissions_constraints.jl

@@ -120,24 +120,35 @@ function calc_yr1_emissions_offset_from_elec_exports(m, p)
 		return 0.0, 0.0, 0.0, 0.0
 	end
 	yr1_emissions_offset_from_elec_exports_lbs_CO2 = @expression(m, p.hours_per_time_step *
-		sum(m[:dvProductionToGrid][t,u,ts] * (p.s.electric_utility.emissions_factor_series_lb_CO2_per_kwh[ts])
-		for t in p.techs.elec, ts in p.time_steps, u in p.export_bins_by_tech[t])
+		sum( p.s.electric_utility.emissions_factor_series_lb_CO2_per_kwh[ts] * (
+			sum(m[:dvProductionToGrid][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t])
+			+ sum(m[:dvStorageToGrid][b, u, ts] for b in p.s.storage.types.elec, u in p.export_bins_by_storage[b])
+			) for ts in p.time_steps
+		)
 	)
-		# if battery ends up being able to discharge to grid, need to incorporate here- might require complex tracking of what's charging battery
 
 	yr1_emissions_offset_from_elec_exports_lbs_NOx = @expression(m, p.hours_per_time_step *
-		sum(m[:dvProductionToGrid][t,u,ts] * (p.s.electric_utility.emissions_factor_series_lb_NOx_per_kwh[ts])
-		for t in p.techs.elec, ts in p.time_steps, u in p.export_bins_by_tech[t])
+		sum( p.s.electric_utility.emissions_factor_series_lb_NOx_per_kwh[ts] * (
+			sum(m[:dvProductionToGrid][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t])
+			+ sum(m[:dvStorageToGrid][b, u, ts] for b in p.s.storage.types.elec, u in p.export_bins_by_storage[b])
+			) for ts in p.time_steps
+		)
 	)
 
 	yr1_emissions_offset_from_elec_exports_lbs_SO2 = @expression(m, p.hours_per_time_step *
-		sum(m[:dvProductionToGrid][t,u,ts] * (p.s.electric_utility.emissions_factor_series_lb_SO2_per_kwh[ts])
-		for t in p.techs.elec, ts in p.time_steps, u in p.export_bins_by_tech[t])
+		sum( p.s.electric_utility.emissions_factor_series_lb_SO2_per_kwh[ts] * (
+			sum(m[:dvProductionToGrid][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t])
+			+ sum(m[:dvStorageToGrid][b, u, ts] for b in p.s.storage.types.elec, u in p.export_bins_by_storage[b])
+			) for ts in p.time_steps
+		)
 	)
 
 	yr1_emissions_offset_from_elec_exports_lbs_PM25 = @expression(m, p.hours_per_time_step *
-		sum(m[:dvProductionToGrid][t,u,ts] * (p.s.electric_utility.emissions_factor_series_lb_PM25_per_kwh[ts])
-		for t in p.techs.elec, ts in p.time_steps, u in p.export_bins_by_tech[t])
+		sum( p.s.electric_utility.emissions_factor_series_lb_PM25_per_kwh[ts] * (
+			sum(m[:dvProductionToGrid][t,u,ts] for t in p.techs.elec, u in p.export_bins_by_tech[t])
+			+ sum(m[:dvStorageToGrid][b, u, ts] for b in p.s.storage.types.elec, u in p.export_bins_by_storage[b])
+			) for ts in p.time_steps
+		)
 	)
 
 	return yr1_emissions_offset_from_elec_exports_lbs_CO2,