diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7e8a5f1d8..223ad0202 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## Guidelines
-- When making a Pull Request into `develop` start a new double-hash header for "Develop - YYYY-MM-DD"
+- When working in feature branch, start a new double-hash header with the name of the branch and record changes under that
+- When merging `develop` into a feature branch, keep the feature branch section and the "Develop" section separate to simplify merge conflicts
+- When making a Pull Request into `develop`, merge the feature branch section into the "Develop" section (if it exists), else rename the feature branch header to "Develop"
 - When making a Pull Request into `master` change "Develop" to the next version number
 
 ### Formatting
@@ -32,6 +34,24 @@ Classify the change according to the following categories:
 ### Removed
 - 80% scaling of battery maintenance costs when using augmentation strategy
 
+## Develop
+### Changed
+- Replace all `1/p.s.settings.time_steps_per_hour` with `p.hours_per_time_step` for simplicity/consistency
+- Rename function `add_storage_sum_constraints` to `add_storage_sum_grid_constraints` for clarity
+### Added
+- Constraints to prevent simultaneous charge/discharge of storage
+- Specify in docstrings that **PV** **max_kw** and **size_kw** are kW-DC
+- Add the Logging package to `test/Project.toml` because it is used in `runtests.jl`
+### Fixed
+- Force **ElectricLoad** **critical_load_kw** to be _nothing_ when **off_grid_flag** is _true_ (**critical_load_fraction** was already being forced to 1, but the user was still able to get around this by providing **critical_load_kw**)
+- Removed looping over storage name in functions `add_hot_thermal_storage_dispatch_constraints` and `add_cold_thermal_storage_dispatch_constraints` because this loop is already done when calling these functions and storage name is passed in as argument `b`
+- Remove extraneous line of code in `results/wind.jl`
+- Change type of **value_of_lost_load** in **FinancialInputs** struct to fix convert error when user provides an _Int_
+- Change international location in "Solar Dataset" test set from Cameroon to Oulu because the locations in the NSRDB have been expanded significantly so there is now an NSRDB point at Cameroon
+- Handle edge case where the values of **outage_start_time_steps** and **outage_durations** makes an outage extend beyond the end of the year. The outage will now wrap around to the beginning of the year.
+- Enforce minimum allowable sizes for ASHP technologies by introducing improved big-M values for segmented size constraints.
+- Removed default values from ASHP functions that calculate minimum allowable size and performance.
+
 ## v0.48.0
 ### Added
 - Added new file `src/core/ASHP.jl` with new technology **ASHP**, which uses electricity as input and provides heating and/or cooling as output; load balancing and technology-specific constraints have been updated and added accordingly
diff --git a/src/constraints/outage_constraints.jl b/src/constraints/outage_constraints.jl
index 49b9244e1..2fdc6e92a 100644
--- a/src/constraints/outage_constraints.jl
+++ b/src/constraints/outage_constraints.jl
@@ -2,8 +2,8 @@
 function add_dv_UnservedLoad_constraints(m,p)
     # Effective load balance
     @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
-        m[:dvUnservedLoad][s, tz, ts] == p.s.electric_load.critical_loads_kw[tz+ts-1]
-        - sum(  m[:dvMGRatedProduction][t, s, tz, ts] * (p.production_factor[t, tz+ts-1] + p.unavailability[t][tz+ts-1]) * p.levelization_factor[t]
+        m[:dvUnservedLoad][s, tz, ts] == p.s.electric_load.critical_loads_kw[time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)]
+        - sum(  m[:dvMGRatedProduction][t, s, tz, ts] * (p.production_factor[t, time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)] + p.unavailability[t][time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)]) * p.levelization_factor[t]
               - m[:dvMGProductionToStorage][t, s, tz, ts] - m[:dvMGCurtail][t, s, tz, ts]
             for t in p.techs.elec
         )
@@ -22,7 +22,7 @@ end
 
 function add_outage_cost_constraints(m,p)
     @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps],
-        m[:dvMaxOutageCost][s] >= p.pwf_e * sum(p.value_of_lost_load_per_kwh[tz+ts-1] * m[:dvUnservedLoad][s, tz, ts] for ts in 1:p.s.electric_utility.outage_durations[s])
+        m[:dvMaxOutageCost][s] >= p.pwf_e * sum(p.value_of_lost_load_per_kwh[time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)] * m[:dvUnservedLoad][s, tz, ts] for ts in 1:p.s.electric_utility.outage_durations[s])
     )
 
     @expression(m, ExpectedOutageCost,
@@ -116,7 +116,7 @@ function add_MG_production_constraints(m,p)
 	# Electrical production sent to storage or export must be less than technology's rated production
 	@constraint(m, [t in p.techs.elec, s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
 		m[:dvMGProductionToStorage][t, s, tz, ts] + m[:dvMGCurtail][t, s, tz, ts] <=
-		(p.production_factor[t, tz+ts-1] + p.unavailability[t][tz+ts-1]) * p.levelization_factor[t] * m[:dvMGRatedProduction][t, s, tz, ts]
+		(p.production_factor[t, time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)] + p.unavailability[t][time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)]) * p.levelization_factor[t] * m[:dvMGRatedProduction][t, s, tz, ts]
     )
 
     @constraint(m, [t in p.techs.elec, s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps], 
@@ -138,7 +138,7 @@ function add_MG_Gen_fuel_burn_constraints(m,p)
     # Define dvMGFuelUsed by summing over outage time_steps.
     @constraint(m, [t in p.techs.gen, s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps],
         m[:dvMGFuelUsed][t, s, tz] == fuel_slope_gal_per_kwhe * p.hours_per_time_step * p.levelization_factor[t] *
-        sum( (p.production_factor[t, tz+ts-1] + p.unavailability[t][tz+ts-1]) * m[:dvMGRatedProduction][t, s, tz, ts] for ts in 1:p.s.electric_utility.outage_durations[s])
+        sum( (p.production_factor[t, time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)] + p.unavailability[t][time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)]) * m[:dvMGRatedProduction][t, s, tz, ts] for ts in 1:p.s.electric_utility.outage_durations[s])
         + fuel_intercept_gal_per_hr * p.hours_per_time_step * 
         sum( m[:binMGGenIsOnInTS][s, tz, ts] for ts in 1:p.s.electric_utility.outage_durations[s])
     )
@@ -287,6 +287,14 @@ function add_MG_storage_dispatch_constraints(m,p)
         )
     )
 
+	# Prevent simultaneous charge and discharge by limitting charging alone to not make the SOC exceed 100%
+    @constraint(m, [ts in p.time_steps_without_grid],
+        m[:dvStorageEnergy]["ElectricStorage"] >= m[:dvMGStoredEnergy][s, tz, ts-1] + p.hours_per_time_step * (  
+            p.s.storage.attr["ElectricStorage"].charge_efficiency * sum(m[:dvMGProductionToStorage][t, s, tz, ts] for t in p.techs.elec) 
+        )
+    )
+
+    # Min SOC
     if p.s.storage.attr["ElectricStorage"].soc_min_applies_during_outages
         # Minimum state of charge
         @constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
diff --git a/src/constraints/storage_constraints.jl b/src/constraints/storage_constraints.jl
index 96089f94a..cbf9219b7 100644
--- a/src/constraints/storage_constraints.jl
+++ b/src/constraints/storage_constraints.jl
@@ -50,7 +50,7 @@ end
 
 function add_elec_storage_dispatch_constraints(m, p, b; _n="")
 				
-	# Constraint (4g): state-of-charge for electrical storage - with grid
+	# Constraint (4g)-1: state-of-charge for electrical storage - with grid
 	@constraint(m, [ts in p.time_steps_with_grid],
         m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (  
             sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) 
@@ -58,13 +58,27 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
             - m[Symbol("dvDischargeFromStorage"*_n)][b,ts] / p.s.storage.attr[b].discharge_efficiency
         )
 	)
-
-	# Constraint (4h): state-of-charge for electrical storage - no grid
+	# Constraint (4g)-2: state-of-charge for electrical storage - no grid
 	@constraint(m, [ts in p.time_steps_without_grid],
         m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (  
             sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.elec) 
             - m[Symbol("dvDischargeFromStorage"*_n)][b, ts] / p.s.storage.attr[b].discharge_efficiency
         )
+    )	
+
+	# Constraint (4h): prevent simultaneous charge and discharge by limitting charging alone to not make the SOC exceed 100%
+    # (4h)-1: with grid
+	@constraint(m, [ts in p.time_steps_with_grid],
+        m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (  
+            sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) 
+            + p.s.storage.attr[b].grid_charge_efficiency * m[Symbol("dvGridToStorage"*_n)][b, ts] 
+        )
+	)
+	# (4h)-2: no grid
+	@constraint(m, [ts in p.time_steps_without_grid],
+        m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (  
+            sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.elec) 
+        )
     )
 
 	# Constraint (4i)-1: Dispatch to electrical storage is no greater than power capacity
@@ -104,16 +118,24 @@ end
 function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")
 
     # Constraint (4j)-1: Reconcile state-of-charge for (hot) thermal storage
-	@constraint(m, [b in p.s.storage.types.hot, 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 -
-        p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+	@constraint(m, [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]
         )
     )
     
+    # Prevent simultaneous charge and discharge by limitting charging alone to not make the SOC exceed 100%
+	@constraint(m, [ts in p.time_steps],
+        m[Symbol("dvStorageEnergy"*_n)][b] >= 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)
+            - p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+        )
+    )
+
     #Constraint (4n)-1: Dispatch to and from thermal storage is no greater than power capacity
-	@constraint(m, [b in p.s.storage.types.hot, ts in p.time_steps],
+	@constraint(m, [ts in p.time_steps],
         m[Symbol("dvStoragePower"*_n)][b] >= 
         sum(m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] + 
         sum(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts] for t in union(p.techs.heating, p.techs.chp))
@@ -122,14 +144,14 @@ function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")
     # TODO missing thermal storage constraints from API ???
 
     # Constraint (4o): Discharge from storage is equal to sum of heat from storage for all qualities
-    @constraint(m, HeatDischargeReconciliation[b in p.s.storage.types.hot, ts in p.time_steps],
+    @constraint(m, HeatDischargeReconciliation[ts in p.time_steps],
         m[Symbol("dvDischargeFromStorage"*_n)][b,ts] == 
         sum(m[Symbol("dvHeatFromStorage"*_n)][b,q,ts] for q in p.heating_loads)
     )
 
     #Do not allow GHP to charge storage
     if !isempty(p.techs.ghp)
-        for b in p.s.storage.types.hot, t in p.techs.ghp, q in p.heating_loads, ts in p.time_steps
+        for t in p.techs.ghp, q in p.heating_loads, ts in p.time_steps
             fix(m[Symbol("dvHeatToStorage"*_n)][b,t,q,ts], 0.0, force=true)
         end
     end
@@ -140,36 +162,44 @@ function add_cold_thermal_storage_dispatch_constraints(m, p, b; _n="")
 
     # Constraint (4f)-2: (Cold) Thermal production sent to storage or grid must be less than technology's rated production
 	if !isempty(p.techs.cooling)
-		@constraint(m, CoolingTechProductionFlowCon[b in p.s.storage.types.cold, t in p.techs.cooling, ts in p.time_steps],
+		@constraint(m, CoolingTechProductionFlowCon[t in p.techs.cooling, ts in p.time_steps],
     	        m[Symbol("dvProductionToStorage"*_n)][b,t,ts]  <=
 				m[Symbol("dvCoolingProduction"*_n)][t,ts]
 				)
 	end
 
     # Constraint (4j)-2: Reconcile state-of-charge for (cold) thermal storage
-	@constraint(m, ColdTESInventoryCon[b in p.s.storage.types.cold, 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) * (
-        sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.cooling) -
-        m[Symbol("dvDischargeFromStorage"*_n)][b,ts]/p.s.storage.attr[b].discharge_efficiency -
-        p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+	@constraint(m, ColdTESInventoryCon[ts in p.time_steps],
+        m[Symbol("dvStoredEnergy"*_n)][b,ts] == m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + p.hours_per_time_step * (
+            sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.cooling) -
+            m[Symbol("dvDischargeFromStorage"*_n)][b,ts]/p.s.storage.attr[b].discharge_efficiency -
+            p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+        )
+    )
+
+    # Prevent simultaneous charge and discharge by limitting charging alone to not make the SOC exceed 100%
+	@constraint(m, [ts in p.time_steps],
+        m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + p.hours_per_time_step * (  
+            p.s.storage.attr[b].charge_efficiency * sum(m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.cooling)
+            - p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
         )
     )
 
     #Constraint (4n)-2: Dispatch to and from thermal storage is no greater than power capacity
-    @constraint(m, [b in p.s.storage.types.cold, ts in p.time_steps],
+    @constraint(m, [ts in p.time_steps],
         m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] + 
         sum(m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.cooling)
     )
 
     #Do not allow GHP to charge storage
     if !isempty(p.techs.ghp)
-        for b in p.s.storage.types.cold, t in p.techs.ghp, ts in p.time_steps
+        for t in p.techs.ghp, ts in p.time_steps
                     fix(m[Symbol("dvProductionToStorage"*_n)][b,t,ts], 0.0, force=true)
         end
     end
 end
 
-function add_storage_sum_constraints(m, p; _n="")
+function add_storage_sum_grid_constraints(m, p; _n="")
 
 	##Constraint (8c): Grid-to-storage no greater than grid purchases 
 	@constraint(m, [ts in p.time_steps_with_grid],
diff --git a/src/constraints/thermal_tech_constraints.jl b/src/constraints/thermal_tech_constraints.jl
index ef3fafc36..3daccdf89 100644
--- a/src/constraints/thermal_tech_constraints.jl
+++ b/src/constraints/thermal_tech_constraints.jl
@@ -15,7 +15,7 @@ function add_boiler_tech_constraints(m, p; _n="")
     )
     if "Boiler" in p.techs.boiler  # ExistingBoiler does not have om_cost_per_kwh
         m[:TotalBoilerPerUnitProdOMCosts] = @expression(m, p.third_party_factor * p.pwf_om *
-            sum(p.s.boiler.om_cost_per_kwh / p.s.settings.time_steps_per_hour *
+            sum(p.s.boiler.om_cost_per_kwh * p.hours_per_time_step *
             m[Symbol("dvHeatingProduction"*_n)]["Boiler",q,ts] for q in p.heating_loads, ts in p.time_steps)
         )
     else
diff --git a/src/core/ashp.jl b/src/core/ashp.jl
index f2f2bfcc0..b29a4fe7f 100644
--- a/src/core/ashp.jl
+++ b/src/core/ashp.jl
@@ -444,14 +444,14 @@ function get_ashp_performance(cop_reference,
                 cf_reference,
                 reference_temps,
                 ambient_temp_degF,
-                back_up_temp_threshold_degF = 10.0
+                back_up_temp_threshold_degF
                 )
 """
 function get_ashp_performance(cop_reference,
     cf_reference,
     reference_temps,
     ambient_temp_degF,
-    back_up_temp_threshold_degF = 10.0
+    back_up_temp_threshold_degF
     )
     num_timesteps = length(ambient_temp_degF)
     cop = zeros(num_timesteps)
@@ -493,9 +493,9 @@ Obtains the default minimum allowable size for ASHP system.  This is calculated
 """
 function get_ashp_default_min_allowable_size(heating_load::Array{<:Real,1}, 
     heating_cf::Array{<:Real,1}, 
-    cooling_load::Array{<:Real,1} = Real[], 
-    cooling_cf::Array{<:Real,1} = Real[],
-    peak_load_thermal_factor::Real = 0.5
+    cooling_load::Array{<:Real,1}, 
+    cooling_cf::Array{<:Real,1},
+    peak_load_thermal_factor::Real
     )
 
     if isempty(cooling_cf)
diff --git a/src/core/bau_inputs.jl b/src/core/bau_inputs.jl
index c8b7661b5..7a01ecf7a 100644
--- a/src/core/bau_inputs.jl
+++ b/src/core/bau_inputs.jl
@@ -287,7 +287,7 @@ function setup_bau_emissions_inputs(p::REoptInputs, s_bau::BAUScenario, generato
         bau_grid_to_load = [max(i,0) for i in bau_grid_to_load]
     end
 
-    bau_grid_emissions_lb_CO2_per_year = sum(p.s.electric_utility.emissions_factor_series_lb_CO2_per_kwh .* bau_grid_to_load) / p.s.settings.time_steps_per_hour
+    bau_grid_emissions_lb_CO2_per_year = sum(p.s.electric_utility.emissions_factor_series_lb_CO2_per_kwh .* bau_grid_to_load) * p.hours_per_time_step
     bau_emissions_lb_CO2_per_year += bau_grid_emissions_lb_CO2_per_year
 
     ## Generator emissions (during outages)
diff --git a/src/core/electric_load.jl b/src/core/electric_load.jl
index 84fd080c7..4f6e3ba29 100644
--- a/src/core/electric_load.jl
+++ b/src/core/electric_load.jl
@@ -105,9 +105,15 @@ mutable struct ElectricLoad  # mutable to adjust (critical_)loads_kw based off o
         min_load_met_annual_fraction::Real = off_grid_flag ? 0.99999 : 1.0 # if off grid, 99.999%, else must be 100%. Applied to each time_step as a % of electric load.
         )
         
-        if off_grid_flag  && !(critical_load_fraction == 1.0)
-            @warn "ElectricLoad critical_load_fraction must be 1.0 (100%) for off-grid scenarios. Any other value will be overriden when `off_grid_flag` is true. If you wish to alter the load profile or load met, adjust the loads_kw or min_load_met_annual_fraction."
-            critical_load_fraction = 1.0
+        if off_grid_flag
+            if !isnothing(critical_loads_kw)
+                @warn "ElectricLoad critical_loads_kw will be ignored because `off_grid_flag` is true. If you wish to alter the load profile or load met, adjust the loads_kw or min_load_met_annual_fraction."
+                critical_loads_kw = nothing
+            end
+            if critical_load_fraction != 1.0
+                @warn "ElectricLoad critical_load_fraction must be 1.0 (100%) for off-grid scenarios. Any other value will be overriden when `off_grid_flag` is true. If you wish to alter the load profile or load met, adjust the loads_kw or min_load_met_annual_fraction."
+                critical_load_fraction = 1.0
+            end
         end
 
         if !(off_grid_flag)
diff --git a/src/core/financial.jl b/src/core/financial.jl
index 234fdc662..105ebe7e2 100644
--- a/src/core/financial.jl
+++ b/src/core/financial.jl
@@ -56,7 +56,7 @@ struct Financial
     owner_tax_rate_fraction::Float64
     owner_discount_rate_fraction::Float64
     analysis_years::Int
-    value_of_lost_load_per_kwh::Union{Array{Float64,1}, Float64}
+    value_of_lost_load_per_kwh::Union{Array{<:Real,1}, Real}
     microgrid_upgrade_cost_fraction::Float64
     macrs_five_year::Array{Float64,1}
     macrs_seven_year::Array{Float64,1}
diff --git a/src/core/pv.jl b/src/core/pv.jl
index 1f90132f1..8ae36b5dc 100644
--- a/src/core/pv.jl
+++ b/src/core/pv.jl
@@ -13,7 +13,7 @@
     location::String="both", # one of ["roof", "ground", "both"]
     existing_kw::Real=0,
     min_kw::Real=0,
-    max_kw::Real=1.0e9, # max new capacity (beyond existing_kw)
+    max_kw::Real=1.0e9, # max new DC capacity (beyond existing_kw)
     installed_cost_per_kw::Real=1790.0,
     om_cost_per_kw::Real=18.0,
     degradation_fraction::Real=0.005,
diff --git a/src/core/reopt.jl b/src/core/reopt.jl
index a901e87c0..662d249cd 100644
--- a/src/core/reopt.jl
+++ b/src/core/reopt.jl
@@ -68,7 +68,7 @@ Solve the model using a `Scenario` or `BAUScenario`.
 function run_reopt(m::JuMP.AbstractModel, s::AbstractScenario)
 	
 	try
-		if s.site.CO2_emissions_reduction_min_fraction > 0.0 || s.site.CO2_emissions_reduction_max_fraction < 1.0
+		if (!isnothing(s.site.CO2_emissions_reduction_min_fraction) && s.site.CO2_emissions_reduction_min_fraction > 0.0) || (!isnothing(s.site.CO2_emissions_reduction_max_fraction) && s.site.CO2_emissions_reduction_max_fraction < 1.0)
 			throw(@error("To constrain CO2 emissions reduction min or max percentages, the optimal and business as usual scenarios must be run in parallel. Use a version of run_reopt() that takes an array of two models."))
 		end
 		run_reopt(m, REoptInputs(s))
@@ -247,7 +247,7 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 	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_constraints(m, p)
+		add_storage_sum_grid_constraints(m, p)
 	end
 
 	add_production_constraints(m, p)
diff --git a/src/core/reopt_inputs.jl b/src/core/reopt_inputs.jl
index 409458e6c..bb1d1bf94 100644
--- a/src/core/reopt_inputs.jl
+++ b/src/core/reopt_inputs.jl
@@ -953,7 +953,14 @@ function setup_ASHPSpaceHeater_inputs(s, max_sizes, min_sizes, cap_cost_slope, o
     if s.ashp.min_allowable_kw > 0.0
         cap_cost_slope["ASHPSpaceHeater"] = s.ashp.installed_cost_per_kw
         push!(segmented_techs, "ASHPSpaceHeater")
-        seg_max_size["ASHPSpaceHeater"] = Dict{Int,Float64}(1 => s.ashp.max_kw)
+        #as a reasonable big-M, assume 10 times the size to serve 100% of peak load.
+        ashp_sh_max_size = get_ashp_default_min_allowable_size(s.space_heating_load.loads_kw, 
+            s.ashp.heating_cf, 
+            s.cooling_load.loads_kw_thermal, 
+            s.ashp.cooling_cf,
+            10.0
+        )
+        seg_max_size["ASHPSpaceHeater"] = Dict{Int,Float64}(1 => min(s.ashp.max_kw, ashp_sh_max_size))
         seg_min_size["ASHPSpaceHeater"] = Dict{Int,Float64}(1 => s.ashp.min_allowable_kw)
         n_segs_by_tech["ASHPSpaceHeater"] = 1
         seg_yint["ASHPSpaceHeater"] = Dict{Int,Float64}(1 => 0.0)
@@ -990,7 +997,14 @@ function setup_ASHPWaterHeater_inputs(s, max_sizes, min_sizes, cap_cost_slope, o
     if s.ashp_wh.min_allowable_kw > 0.0
         cap_cost_slope["ASHPWaterHeater"] = s.ashp_wh.installed_cost_per_kw
         push!(segmented_techs, "ASHPWaterHeater")
-        seg_max_size["ASHPWaterHeater"] = Dict{Int,Float64}(1 => s.ashp_wh.max_kw)
+        #as a reasonable big-M, assume 10 times the size to serve 100% of peak load.
+        ashp_wh_max_size = get_ashp_default_min_allowable_size(s.dhw_load.loads_kw, 
+            s.ashp_wh.heating_cf, 
+            Real[], 
+            Real[],
+            10.0
+        )
+        seg_max_size["ASHPWaterHeater"] = Dict{Int,Float64}(1 => min(s.ashp_wh.max_kw, ashp_wh_max_size))
         seg_min_size["ASHPWaterHeater"] = Dict{Int,Float64}(1 => s.ashp_wh.min_allowable_kw)
         n_segs_by_tech["ASHPWaterHeater"] = 1
         seg_yint["ASHPWaterHeater"] = Dict{Int,Float64}(1 => 0.0)
diff --git a/src/core/reopt_multinode.jl b/src/core/reopt_multinode.jl
index 240956f6e..a41a80a35 100644
--- a/src/core/reopt_multinode.jl
+++ b/src/core/reopt_multinode.jl
@@ -134,7 +134,7 @@ function build_reopt!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}})
         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_constraints(m, p; _n=_n)
+            add_storage_sum_grid_constraints(m, p; _n=_n)
         end
     
         add_production_constraints(m, p; _n=_n)
diff --git a/src/core/scenario.jl b/src/core/scenario.jl
index 5e5c83c42..621781f53 100644
--- a/src/core/scenario.jl
+++ b/src/core/scenario.jl
@@ -40,11 +40,15 @@ A Scenario struct can contain the following keys:
 - [PV](@ref) (optional, can be Array)
 - [Wind](@ref) (optional)
 - [ElectricStorage](@ref) (optional)
+- [HotThermalStorage](@ref) (optional)
+- [ColdThermalStorage](@ref) (optional)
+- [ElectricStorage](@ref) (optional)
 - [ElectricUtility](@ref) (optional)
 - [Generator](@ref) (optional)
 - [DomesticHotWaterLoad](@ref) (optional)
 - [SpaceHeatingLoad](@ref) (optional)
 - [ProcessHeatLoad](@ref) (optional)
+- [CoolingLoad](@ref) (optional)
 - [ExistingBoiler](@ref) (optional)
 - [Boiler](@ref) (optional)
 - [CHP](@ref) (optional)
diff --git a/src/core/utils.jl b/src/core/utils.jl
index 2ba9f4d4b..8f92c9f6f 100644
--- a/src/core/utils.jl
+++ b/src/core/utils.jl
@@ -1,4 +1,9 @@
 # REopt®, Copyright (c) Alliance for Sustainable Energy, LLC. See also https://github.com/NREL/REopt.jl/blob/master/LICENSE.
+function time_step_wrap_around(time_step::Int; time_steps_per_hour::Int=1)::Int
+    time_steps_per_year = 8760 * time_steps_per_hour
+    ((time_step - 1) % time_steps_per_year) + 1
+end
+
 function solver_is_compatible_with_indicator_constraints(solver_name::String)::Bool
     return any(lowercase.(INDICATOR_COMPATIBLE_SOLVERS) .== lowercase(solver_name))
 end
diff --git a/src/mpc/model.jl b/src/mpc/model.jl
index dc8145033..d40a80da3 100644
--- a/src/mpc/model.jl
+++ b/src/mpc/model.jl
@@ -115,7 +115,7 @@ function build_mpc!(m::JuMP.AbstractModel, p::MPCInputs)
 	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_constraints(m, p)
+		add_storage_sum_grid_constraints(m, p)
 	end
 
 	add_production_constraints(m, p)
diff --git a/src/mpc/model_multinode.jl b/src/mpc/model_multinode.jl
index ab71d5765..68a3ef34a 100644
--- a/src/mpc/model_multinode.jl
+++ b/src/mpc/model_multinode.jl
@@ -82,7 +82,7 @@ function build_mpc!(m::JuMP.AbstractModel, ps::AbstractVector{MPCInputs})
 		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_constraints(m, p; _n=_n)
+			add_storage_sum_grid_constraints(m, p; _n=_n)
 		end
 
 		add_production_constraints(m, p; _n=_n)
diff --git a/src/results/electric_load.jl b/src/results/electric_load.jl
index 3024f6943..761e20658 100644
--- a/src/results/electric_load.jl
+++ b/src/results/electric_load.jl
@@ -23,7 +23,7 @@ function add_electric_load_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dic
     r["load_series_kw"] = p.s.electric_load.loads_kw
     r["critical_load_series_kw"] = p.s.electric_load.critical_loads_kw
     r["annual_calculated_kwh"] = round(
-        sum(r["load_series_kw"]) / p.s.settings.time_steps_per_hour, digits=2
+        sum(r["load_series_kw"]) * p.hours_per_time_step, digits=2
     )
     
     if p.s.settings.off_grid_flag
diff --git a/src/results/heating_cooling_load.jl b/src/results/heating_cooling_load.jl
index 1a1d9fe19..4c9a270be 100644
--- a/src/results/heating_cooling_load.jl
+++ b/src/results/heating_cooling_load.jl
@@ -23,11 +23,11 @@ function add_cooling_load_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict
     r["electric_chiller_base_load_series_kw"] = load_series_kw ./ p.s.cooling_load.existing_chiller_cop
 
     r["annual_calculated_tonhour"] = round(
-        sum(r["load_series_ton"]) / p.s.settings.time_steps_per_hour, digits=2
+        sum(r["load_series_ton"]) * p.hours_per_time_step, digits=2
     )
     
     r["annual_electric_chiller_base_load_kwh"] = round(
-        sum(r["electric_chiller_base_load_series_kw"]) / p.s.settings.time_steps_per_hour, digits=2
+        sum(r["electric_chiller_base_load_series_kw"]) * p.hours_per_time_step, digits=2
     )
 
     d["CoolingLoad"] = r
@@ -76,13 +76,13 @@ function add_heating_load_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict
     r["total_heating_boiler_fuel_load_series_mmbtu_per_hour"] =   r["dhw_boiler_fuel_load_series_mmbtu_per_hour"] .+ r["space_heating_boiler_fuel_load_series_mmbtu_per_hour"] .+ r["process_heat_boiler_fuel_load_series_mmbtu_per_hour"]
 
     r["annual_calculated_dhw_thermal_load_mmbtu"] = round(
-        sum(r["dhw_thermal_load_series_mmbtu_per_hour"]) / p.s.settings.time_steps_per_hour, digits =   2
+        sum(r["dhw_thermal_load_series_mmbtu_per_hour"]) * p.hours_per_time_step, digits =   2
     )
     r["annual_calculated_space_heating_thermal_load_mmbtu"] = round(
-        sum(r["space_heating_thermal_load_series_mmbtu_per_hour"]) / p.s.settings.time_steps_per_hour, digits =   2
+        sum(r["space_heating_thermal_load_series_mmbtu_per_hour"]) * p.hours_per_time_step, digits =   2
     )
     r["annual_calculated_process_heat_thermal_load_mmbtu"] = round(
-        sum(r["process_heat_thermal_load_series_mmbtu_per_hour"]) / p.s.settings.time_steps_per_hour, digits =   2
+        sum(r["process_heat_thermal_load_series_mmbtu_per_hour"]) * p.hours_per_time_step, digits =   2
     )
     r["annual_calculated_total_heating_thermal_load_mmbtu"] =   r["annual_calculated_dhw_thermal_load_mmbtu"] + r["annual_calculated_space_heating_thermal_load_mmbtu"] + r["annual_calculated_process_heat_thermal_load_mmbtu"]
     
diff --git a/src/results/outages.jl b/src/results/outages.jl
index 318ce766b..be2bfd1a2 100644
--- a/src/results/outages.jl
+++ b/src/results/outages.jl
@@ -37,7 +37,7 @@
 	The output keys for "Outages" are subject to change.
 
 !!! note 
-	`Outage` results only added to results when multiple outages are modeled via the `ElectricUtility.outage_durations` input.
+	`Outage` results only added to results when multiple outages are modeled via the `ElectricUtility.outage_start_time_steps` and `ElectricUtility.outage_durations` inputs.
 
 !!! note
 	When modeling PV the name of the PV system is used for the output keys to allow for modeling multiple PV systems. The default PV name is `PV`.
@@ -88,7 +88,7 @@ function add_outage_results(m, p, d::Dict)
     for ts in p.s.electric_utility.outage_time_steps
         for (t, tz) in enumerate(p.s.electric_utility.outage_start_time_steps)
             for s in p.s.electric_utility.scenarios
-                r["critical_loads_per_outage_series_kw"][s,t,ts] = p.s.electric_load.critical_loads_kw[tz+ts-1]
+                r["critical_loads_per_outage_series_kw"][s,t,ts] = p.s.electric_load.critical_loads_kw[time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)]
             end
         end
     end
@@ -146,7 +146,7 @@ function add_outage_results(m, p, d::Dict)
 				sum(
 					(
 						value.(
-							m[:dvMGRatedProduction][t, s, tz, ts] * (p.production_factor[t, tz+ts-1] + p.unavailability[t][tz+ts-1]) * p.levelization_factor[t]
+							m[:dvMGRatedProduction][t, s, tz, ts] * (p.production_factor[t, time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)] + p.unavailability[t][time_step_wrap_around(tz+ts-1, time_steps_per_hour=p.s.settings.time_steps_per_hour)]) * p.levelization_factor[t]
 							- m[:dvMGCurtail][t, s, tz, ts]
 							- m[:dvMGProductionToStorage][t, s, tz, ts]
 							for s in p.s.electric_utility.scenarios,
diff --git a/src/results/pv.jl b/src/results/pv.jl
index 439f5e102..4c5774faa 100644
--- a/src/results/pv.jl
+++ b/src/results/pv.jl
@@ -1,7 +1,7 @@
 # REopt®, Copyright (c) Alliance for Sustainable Energy, LLC. See also https://github.com/NREL/REopt.jl/blob/master/LICENSE.
 """
 `PV` results keys:
-- `size_kw` Optimal PV capacity
+- `size_kw` Optimal PV DC capacity
 - `lifecycle_om_cost_after_tax` Lifecycle operations and maintenance cost in present value, after tax
 - `year_one_energy_produced_kwh` Energy produced over the first year
 - `annual_energy_produced_kwh` Average annual energy produced when accounting for degradation
diff --git a/src/results/wind.jl b/src/results/wind.jl
index fe16dc065..0e89e19dc 100644
--- a/src/results/wind.jl
+++ b/src/results/wind.jl
@@ -32,8 +32,6 @@ function add_wind_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 
 	if !isempty(p.s.storage.types.elec)
 		WindToStorage = (sum(m[:dvProductionToStorage][b, t, ts] for b in p.s.storage.types.elec) for ts in p.time_steps)
-		PVtoBatt = (sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for b in p.s.storage.types.elec) for ts in p.time_steps)
-
 	else
 		WindToStorage = zeros(length(p.time_steps))
 	end
diff --git a/test/Project.toml b/test/Project.toml
index dad017e0e..066554eb4 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,6 +4,7 @@ DotEnv = "4dc1fcf4-5e3b-5448-94ab-0c38ec0385c1"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Xpress = "9e70acf3-d6c9-5be6-b5bd-4e2c73e3e054"
diff --git a/test/runtests.jl b/test/runtests.jl
index 6976abac5..a998c75b2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,9 +21,29 @@ elseif "CPLEX" in ARGS
     @testset "test_with_cplex" begin
         include("test_with_cplex.jl")
     end
-
 else  # run HiGHS tests
     @testset verbose=true "REopt test set using HiGHS solver" begin
+        @testset "Prevent simultaneous charge and discharge" begin
+            logger = SimpleLogger()
+            results = nothing
+            with_logger(logger) do
+                model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+                results = run_reopt(model, "./scenarios/simultaneous_charge_discharge.json")
+            end
+            @test any(.&(
+                    results["ElectricStorage"]["storage_to_load_series_kw"] .!= 0.0,
+                    (
+                        results["ElectricUtility"]["electric_to_storage_series_kw"] .+ 
+                        results["PV"]["electric_to_storage_series_kw"]
+                    ) .!= 0.0
+                )
+                ) ≈ false
+            @test any(.&(
+                    results["Outages"]["storage_discharge_series_kw"] .!= 0.0,
+                    results["Outages"]["pv_to_storage_series_kw"] .!= 0.0
+                )
+                ) ≈ false
+        end
         @testset "hybrid profile" begin
             electric_load = REopt.ElectricLoad(; 
                 blended_doe_reference_percents = [0.2, 0.2, 0.2, 0.2, 0.2],
@@ -50,12 +70,11 @@ else  # run HiGHS tests
             dataset, distance, datasource = REopt.call_solar_dataset_api(latitude, longitude, radius)
             @test dataset == "nsrdb"
 
-            # 3. Younde, Cameroon
-            latitude, longitude = 3.8603988398663125, 11.528880303663136
+            # 3. Oulu, Findland
+            latitude, longitude = 65.0102196310875, 25.465387094897675
             radius = 0
             dataset, distance, datasource = REopt.call_solar_dataset_api(latitude, longitude, radius)
-            @test dataset ≈ "nsrdb"
-            @test dataset == "nsrdb"
+            @test dataset ≈ "intl"
 
             # 4. Fairbanks, AK 
             site = "Fairbanks"
@@ -99,8 +118,8 @@ else  # run HiGHS tests
             # min allowable size
             heating_load = Array{Real}([10.0,10.0,10.0,10.0])
             cooling_load = Array{Real}([10.0,10.0,10.0,10.0])
-            space_heating_min_allowable_size = REopt.get_ashp_default_min_allowable_size(heating_load, heating_cf, cooling_load, cooling_cf)
-            wh_min_allowable_size = REopt.get_ashp_default_min_allowable_size(heating_load, heating_cf)
+            space_heating_min_allowable_size = REopt.get_ashp_default_min_allowable_size(heating_load, heating_cf, cooling_load, cooling_cf, 0.5)
+            wh_min_allowable_size = REopt.get_ashp_default_min_allowable_size(heating_load, heating_cf, Real[], Real[], 0.5)
             @test space_heating_min_allowable_size ≈ 9.166666666666666 atol=1e-8
             @test wh_min_allowable_size ≈ 5.0 atol=1e-8
         end
@@ -141,17 +160,6 @@ else  # run HiGHS tests
             @test r["ElectricStorage"]["size_kwh"] ≈ 83.3 atol=0.1
         end
 
-        @testset "Outage with Generator" begin
-            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
-            results = run_reopt(model, "./scenarios/generator.json")
-            @test results["Generator"]["size_kw"] ≈ 9.55 atol=0.01
-            @test (sum(results["Generator"]["electric_to_load_series_kw"][i] for i in 1:9) + 
-                sum(results["Generator"]["electric_to_load_series_kw"][i] for i in 13:8760)) == 0
-            p = REoptInputs("./scenarios/generator.json")
-            simresults = simulate_outages(results, p)
-            @test simresults["resilience_hours_max"] == 11
-        end
-
         # TODO test MPC with outages
         @testset "MPC" begin
             model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
@@ -159,6 +167,10 @@ else  # run HiGHS tests
             @test maximum(r["ElectricUtility"]["to_load_series_kw"][1:15]) <= 98.0 
             @test maximum(r["ElectricUtility"]["to_load_series_kw"][16:24]) <= 97.0
             @test sum(r["PV"]["to_grid_series_kw"]) ≈ 0
+            grid_draw = r["ElectricUtility"]["to_load_series_kw"] .+ r["ElectricUtility"]["to_battery_series_kw"]
+            # the grid draw limit in the 10th time step is set to 90
+            # without the 90 limit the grid draw is 98 in the 10th time step
+            @test grid_draw[10] <= 90
         end
 
         @testset "MPC Multi-node" begin
@@ -260,10 +272,7 @@ else  # run HiGHS tests
             post["PV"]["tilt"] = 17
             scen = Scenario(post)
             @test scen.pvs[1].tilt ≈ 17
-        
-        
-        end    
-        
+        end
 
         @testset "AlternativeFlatLoads" begin
             input_data = JSON.parsefile("./scenarios/flatloads.json")
@@ -1211,30 +1220,6 @@ else  # run HiGHS tests
             @test results[3]["Financial"]["lcc"] + results[10]["Financial"]["lcc"] ≈ 1.2830872235e7 rtol=1e-5
         end
 
-        @testset "MPC" begin
-            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
-            r = run_mpc(model, "./scenarios/mpc.json")
-            @test maximum(r["ElectricUtility"]["to_load_series_kw"][1:15]) <= 98.0 
-            @test maximum(r["ElectricUtility"]["to_load_series_kw"][16:24]) <= 97.0
-            @test sum(r["PV"]["to_grid_series_kw"]) ≈ 0
-            grid_draw = r["ElectricUtility"]["to_load_series_kw"] .+ r["ElectricUtility"]["to_battery_series_kw"]
-            # the grid draw limit in the 10th time step is set to 90
-            # without the 90 limit the grid draw is 98 in the 10th time step
-            @test grid_draw[10] <= 90
-        end
-
-        @testset "Complex Incentives" begin
-            """
-            This test was compared against the API test:
-                reo.tests.test_reopt_url.EntryResourceTest.test_complex_incentives
-            when using the hardcoded levelization_factor in this package's REoptInputs function.
-            The two LCC's matched within 0.00005%. (The Julia pkg LCC is  1.0971991e7)
-            """
-            m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
-            results = run_reopt(m, "./scenarios/incentives.json")
-            @test results["Financial"]["lcc"] ≈ 1.094596365e7 atol=5e4  
-        end
-
         @testset verbose=true "Rate Structures" begin
 
             @testset "Tiered Energy" begin
@@ -1848,8 +1833,8 @@ else  # run HiGHS tests
             post_name = "off_grid.json" 
             post = JSON.parsefile("./scenarios/$post_name")
             m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
-            r = run_reopt(m, post)
             scen = Scenario(post)
+            r = run_reopt(m, scen)
             
             # Test default values 
             @test scen.electric_utility.outage_start_time_step ≈ 1
diff --git a/test/scenarios/nogridcost_multiscenario.json b/test/scenarios/nogridcost_multiscenario.json
index 29a0eb83c..14eb6be34 100644
--- a/test/scenarios/nogridcost_multiscenario.json
+++ b/test/scenarios/nogridcost_multiscenario.json
@@ -5,7 +5,7 @@
       5
     ],
     "outage_start_time_steps": [
-      10
+      8758
     ],
     "outage_probabilities": [
       0.5,
diff --git a/test/scenarios/simultaneous_charge_discharge.json b/test/scenarios/simultaneous_charge_discharge.json
new file mode 100644
index 000000000..ec82ffcfc
--- /dev/null
+++ b/test/scenarios/simultaneous_charge_discharge.json
@@ -0,0 +1,32 @@
+{
+    "ElectricLoad": {
+        "doe_reference_name": "LargeHotel"
+    },
+    "Site": {
+        "latitude": 33.068653, 
+        "longitude": -116.546593,
+        "land_acres": 1000000,
+        "min_resil_time_steps": 0
+    },
+    "ElectricTariff": {
+        "blended_annual_demand_rate": 10.0,
+        "blended_annual_energy_rate": 0.2
+    },
+    "PV": {
+        "min_kw": 2000,
+        "max_kw": 2000
+    },
+    "ElectricStorage": {
+        "min_kw": 300,
+        "min_kwh": 600,
+        "max_kw": 300,
+        "max_kwh": 600
+    },
+    "ElectricUtility": {
+        "outage_start_time_steps": [12, 3612],
+        "outage_durations": [24]
+    },
+    "Financial": {
+        "value_of_lost_load_per_kwh": 10000
+    }
+}
\ No newline at end of file