diff --git a/CHANGELOG.md b/CHANGELOG.md
index 299463000..c879bca0d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -23,6 +23,20 @@ Classify the change according to the following categories:
     ### Deprecated
     ### Removed
 
+## Develop LDES
+### Added
+- Added inputs **om_cost_per_kw** and **om_cost_per_kwh** to `ElectricStorage` for modeling capacity-based O&M 
+- Added outputs **lifecycle_om_cost_after_tax** and **year_one_om_cost_before_tax** to `ElectricStorage` 
+- Added new input **self_discharge_fraction_per_timestep** to `ElectricStorage` for modeling battery self-discharge
+- Added input option **can_export_to_grid** (defaults to _false_) to `ElectricStorage` and decision variable **dvStorageToGrid**
+- Added input **allow_simultaneous_charge_discharge** (defaults to _true_) to `ElectricStorage` to give option to disallow battery from simultaneously charging and discharging, which adds binary decision variables often unnecessarily
+- Added input **fixed_duration** to `ElectricStorage`, which fixes the ratio between **size_kw** and **size_kwh** in the optimized solution if provided
+- Added input option **optimize_soc_init_fraction** (defaults to _false_) to `ElectricStorage`, which makes the optimization choose the inital SOC (equal to final SOC) instead of using **soc_init_fraction**. The initial SOC is also constrained to equal the final SOC, which eliminates the "free energy" issue. We currently do not fix SOC when **soc_init_fraction** is used because this has caused infeasibility. 
+- Added output **initial_capital_cost** to all techs
+### Fixed
+- Added missing outputs **lifecycle_export_benefit_before_tax** and **year_one_export_benefit_after_tax** to `ElectricTariff`
+- Add missing output **year_one_om_cost_before_tax** to `PV`
+
 ## Develop 08-09-2024
 ### Changed
 - Improve the full test suite reporting with a verbose summary table, and update the structure to reflect long-term open-source solver usage
diff --git a/src/constraints/electric_utility_constraints.jl b/src/constraints/electric_utility_constraints.jl
index 201317e8f..3c3321ff7 100644
--- a/src/constraints/electric_utility_constraints.jl
+++ b/src/constraints/electric_utility_constraints.jl
@@ -129,8 +129,14 @@ function add_export_constraints(m, p; _n="")
         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( 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("dvStorageToGrid"*_n)][b, ts] for b in p.s.storage.types.elec) + 
+                        sum(m[Symbol("dvProductionToGrid"*_n)][t, :WHL, ts] for t in p.techs_by_exportbin[:WHL])
+                    ) for ts in p.time_steps
+                )
             )
         else
             binWHL = @variable(m, binary = true)
@@ -273,13 +279,13 @@ function add_simultaneous_export_import_constraint(m, p; _n="")
             }
         )
     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_power = 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)+p.s.storage.attr["ElectricStorage"].max_kw
         @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_power*(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]) <= bigM_hourly_load_plus_battery_power * m[Symbol("binNoGridPurchases"*_n)][ts]
         )
     end
 end
diff --git a/src/constraints/storage_constraints.jl b/src/constraints/storage_constraints.jl
index 96089f94a..050953b0c 100644
--- a/src/constraints/storage_constraints.jl
+++ b/src/constraints/storage_constraints.jl
@@ -25,10 +25,21 @@ end
 
 
 function add_general_storage_dispatch_constraints(m, p, b; _n="")
-    # Constraint (4a): initial state of charge
-	@constraint(m,
-        m[Symbol("dvStoredEnergy"*_n)][b, 0] == p.s.storage.attr[b].soc_init_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
-    )
+    # Constraint (4a): initial and final state of charge
+    if hasproperty(p.s.storage.attr[b], :optimize_soc_init_fraction) && p.s.storage.attr[b].optimize_soc_init_fraction
+        # @constraint(m, m[:dvStoredEnergy]["ElectricStorage",maximum(p.time_steps)] == p.s.storage.attr[b].soc_init_fraction * m[Symbol("dvStorageEnergy"*_n)][b] )
+        print("\nOptimizing "*b*" inital SOC and constraining initial SOC = final SOC\n")
+        @constraint(m,
+            m[Symbol("dvStoredEnergy"*_n)][b, 0] == m[:dvStoredEnergy]["ElectricStorage", maximum(p.time_steps)]
+        )
+    else
+        @constraint(m,
+            m[Symbol("dvStoredEnergy"*_n)][b, 0] == p.s.storage.attr[b].soc_init_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+        )
+        # @constraint(m,
+        #     m[Symbol("dvStoredEnergy"*_n)][b, maximum(p.time_steps)] == p.s.storage.attr[b].soc_init_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
+        # )
+    end
 
     #Constraint (4n): State of charge upper bound is storage system size
     @constraint(m, [ts in p.time_steps],
@@ -52,16 +63,18 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
 				
 	# Constraint (4g): 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 * (  
+        m[Symbol("dvStoredEnergy"*_n)][b, ts] == (1-p.s.storage.attr[b].self_discharge_fraction_per_timestep) * 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] 
-            - m[Symbol("dvDischargeFromStorage"*_n)][b,ts] / p.s.storage.attr[b].discharge_efficiency
+            - ((m[Symbol("dvDischargeFromStorage"*_n)][b,ts]+m[Symbol("dvStorageToGrid"*_n)][b, ts])  / p.s.storage.attr[b].discharge_efficiency)
         )
 	)
 
 	# Constraint (4h): 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 * (  
+        m[Symbol("dvStoredEnergy"*_n)][b, ts] == (1-p.s.storage.attr[b].self_discharge_fraction_per_timestep) * 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
         )
@@ -77,21 +90,33 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
 	@constraint(m, [ts in p.time_steps_with_grid],
         m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b, ts] + 
             sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) + m[Symbol("dvGridToStorage"*_n)][b, ts]
+            + m[Symbol("dvStorageToGrid"*_n)][b, ts]
     )
 
+    #Dispatch from electrical storage is no greater than power capacity 
+    @constraint(m, [ts in p.time_steps_without_grid],
+        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] + m[Symbol("dvStorageToGrid"*_n)][b, ts])
+
 	#Constraint (4l)-alt: Dispatch from electrical storage is no greater than power capacity (no grid connection)
 	@constraint(m, [ts in p.time_steps_without_grid],
         m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] + 
             sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec)
     )
-					
-    # Remove grid-to-storage as an option if option to grid charge is turned off
+			
+    #Constraint (4m)-1: Remove grid-to-storage as an option if option to grid charge is turned off
     if !(p.s.storage.attr[b].can_grid_charge)
         for ts in p.time_steps_with_grid
             fix(m[Symbol("dvGridToStorage"*_n)][b, ts], 0.0, force=true)
         end
 	end
 
+    #Constraint (4m)-2: Force storage export to grid to zero if option to grid export is turned off
+    if !p.s.storage.attr[b].can_export_to_grid
+        for ts in p.time_steps
+            fix(m[Symbol("dvStorageToGrid"*_n)][b, ts], 0.0, force=true)
+        end
+    end
+
     if p.s.storage.attr[b].minimum_avg_soc_fraction > 0
         avg_soc = sum(m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps) /
                    (8760. / p.hours_per_time_step)
@@ -99,6 +124,25 @@ function add_elec_storage_dispatch_constraints(m, p, b; _n="")
             sum(m[Symbol("dvStorageEnergy"*_n)][b])
         )
     end
+
+    if p.s.storage.attr[b] isa ElectricStorage && !isnothing(p.s.storage.attr[b].fixed_duration)
+        @constraint(m, m[Symbol("dvStoragePower"*_n)][b] == m[Symbol("dvStorageEnergy"*_n)][b] / p.s.storage.attr[b].fixed_duration)
+    end
+
+    # Prevent charging and discharging of the battery at the same time
+    if !(p.s.storage.attr[b].allow_simultaneous_charge_discharge)
+        #TODO: implement indicator constraint version for solvers that support it
+        @constraint(m, [ts in p.time_steps],
+                    m[Symbol("dvGridToStorage"*_n)][b, ts] + 
+                    sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) <=
+                    p.s.storage.attr[b].max_kw * m[Symbol("binBattCharging")][ts])
+        @constraint(m, [ts in p.time_steps], 
+                    m[Symbol("dvStorageToGrid"*_n)][b, ts] +
+                    m[Symbol("dvDischargeFromStorage"*_n)][b, ts] <= 
+                    p.s.storage.attr[b].max_kw * (1-m[Symbol("binBattCharging")][ts]))
+    end
+
+
 end
 
 function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")
diff --git a/src/core/energy_storage/electric_storage.jl b/src/core/energy_storage/electric_storage.jl
index 626e2993a..63cb028df 100644
--- a/src/core/energy_storage/electric_storage.jl
+++ b/src/core/energy_storage/electric_storage.jl
@@ -147,10 +147,13 @@ end
     soc_min_applies_during_outages::Bool = false
     soc_init_fraction::Float64 = off_grid_flag ? 1.0 : 0.5
     can_grid_charge::Bool = off_grid_flag ? false : true
+    can_export_to_grid::Bool = false
     installed_cost_per_kw::Real = 910.0
     installed_cost_per_kwh::Real = 455.0
     replace_cost_per_kw::Real = 715.0
     replace_cost_per_kwh::Real = 318.0
+    om_cost_per_kw::Real = 0.0 # Capacity-based O&M costs in \$/kW-rated/year. When including battery replacement costs, the user should ensure that O&M costs do not account for augmentation/replacement.
+    om_cost_per_kwh::Real = 0.0 # Capacity-based O&M costs in \$/kWh-rated/year. When including battery replacement costs, the user should ensure that O&M costs do not account for augmentation/replacement.
     inverter_replacement_year::Int = 10
     battery_replacement_year::Int = 10
     macrs_option_years::Int = 7
@@ -165,7 +168,11 @@ end
     model_degradation::Bool = false
     degradation::Dict = Dict()
     minimum_avg_soc_fraction::Float64 = 0.0
-```
+    self_discharge_fraction_per_timestep::Float64 = 0.0 # Battery self-discharge as a fraction per timestep loss based on kWh stored in each timestep
+    fixed_duration::Union{Real, Nothing} = nothing
+    optimize_soc_init_fraction::Bool = false
+    allow_simultaneous_charge_discharge::Bool = true # Simultaneous charge/discharge is typically suboptimal anyway and allowing this avoids additional binary variables (which can slow solve time)
+``` 
 """
 Base.@kwdef struct ElectricStorageDefaults
     off_grid_flag::Bool = false
@@ -180,10 +187,13 @@ Base.@kwdef struct ElectricStorageDefaults
     soc_min_applies_during_outages::Bool = false
     soc_init_fraction::Float64 = off_grid_flag ? 1.0 : 0.5
     can_grid_charge::Bool = off_grid_flag ? false : true
+    can_export_to_grid::Bool = false
     installed_cost_per_kw::Real = 910.0
     installed_cost_per_kwh::Real = 455.0
     replace_cost_per_kw::Real = 715.0
     replace_cost_per_kwh::Real = 318.0
+    om_cost_per_kw::Real = 0.0
+    om_cost_per_kwh::Real = 0.0
     inverter_replacement_year::Int = 10
     battery_replacement_year::Int = 10
     macrs_option_years::Int = 7
@@ -198,6 +208,10 @@ Base.@kwdef struct ElectricStorageDefaults
     model_degradation::Bool = false
     degradation::Dict = Dict()
     minimum_avg_soc_fraction::Float64 = 0.0
+    self_discharge_fraction_per_timestep::Float64 = 0.0
+    fixed_duration::Union{Real, Nothing} = nothing
+    optimize_soc_init_fraction::Bool = false
+    allow_simultaneous_charge_discharge::Bool = true
 end
 
 
@@ -219,10 +233,13 @@ struct ElectricStorage <: AbstractElectricStorage
     soc_min_applies_during_outages::Bool
     soc_init_fraction::Float64
     can_grid_charge::Bool
+    can_export_to_grid::Bool
     installed_cost_per_kw::Real
     installed_cost_per_kwh::Real
     replace_cost_per_kw::Real
     replace_cost_per_kwh::Real
+    om_cost_per_kw::Real
+    om_cost_per_kwh::Real
     inverter_replacement_year::Int
     battery_replacement_year::Int
     macrs_option_years::Int
@@ -239,6 +256,10 @@ struct ElectricStorage <: AbstractElectricStorage
     model_degradation::Bool
     degradation::Degradation
     minimum_avg_soc_fraction::Float64
+    self_discharge_fraction_per_timestep::Float64
+    fixed_duration::Union{Real, Nothing}
+    optimize_soc_init_fraction::Bool
+    allow_simultaneous_charge_discharge::Bool
 
     function ElectricStorage(d::Dict, f::Financial)  
         s = ElectricStorageDefaults(;d...)
@@ -307,10 +328,13 @@ struct ElectricStorage <: AbstractElectricStorage
             s.soc_min_applies_during_outages,
             s.soc_init_fraction,
             s.can_grid_charge,
+            s.can_export_to_grid,
             s.installed_cost_per_kw,
             s.installed_cost_per_kwh,
             replace_cost_per_kw,
             replace_cost_per_kwh,
+            s.om_cost_per_kw,
+            s.om_cost_per_kwh,
             s.inverter_replacement_year,
             s.battery_replacement_year,
             s.macrs_option_years,
@@ -326,7 +350,11 @@ struct ElectricStorage <: AbstractElectricStorage
             net_present_cost_per_kwh,
             s.model_degradation,
             degr,
-            s.minimum_avg_soc_fraction
+            s.minimum_avg_soc_fraction,
+            s.self_discharge_fraction_per_timestep,
+            s.fixed_duration,
+            s.optimize_soc_init_fraction,
+            s.allow_simultaneous_charge_discharge
         )
     end
 end
diff --git a/src/core/reopt.jl b/src/core/reopt.jl
index e9344b703..f3d081c07 100644
--- a/src/core/reopt.jl
+++ b/src/core/reopt.jl
@@ -213,6 +213,7 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 				@constraint(m, [ts in p.time_steps], m[:dvGridToStorage][b, ts] == 0)
 				@constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
 						m[:dvProductionToStorage][b, t, ts] == 0)
+				@constraint(m, [ts in p.time_steps], m[Symbol("dvStorageToGrid"*_n)][b, ts] == 0)  # if there isn't a battery, then the battery can't export power to the grid
 			elseif b in p.s.storage.types.hot
 				@constraint(m, [q in q in setdiff(p.heating_loads, p.heating_loads_served_by_tes[b]), ts in p.time_steps], m[:dvHeatFromStorage][b,q,ts] == 0)
 				if "DomesticHotWater" in p.heating_loads_served_by_tes[b]
@@ -392,9 +393,11 @@ function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
 		sum( p.s.storage.attr[b].net_present_cost_per_kwh * m[:dvStorageEnergy][b] for b in p.s.storage.types.all )
 	))
 	
-	@expression(m, TotalPerUnitSizeOMCosts, p.third_party_factor * p.pwf_om *
-		sum( p.om_cost_per_kw[t] * m[:dvSize][t] for t in p.techs.all )
-	)
+	@expression(m, TotalPerUnitSizeOMCosts, p.third_party_factor * p.pwf_om * (
+		sum(p.om_cost_per_kw[t] * m[:dvSize][t] for t in p.techs.all) + 
+		sum(p.s.storage.attr[b].om_cost_per_kw * m[:dvStoragePower][b] for b in p.s.storage.types.elec) +
+		sum(p.s.storage.attr[b].om_cost_per_kwh * m[:dvStorageEnergy][b] for b in p.s.storage.types.elec)
+	))
 
 	add_elec_utility_expressions(m, p)
 
@@ -587,6 +590,7 @@ function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
 		dvProductionToStorage[p.s.storage.types.all, union(p.techs.ghp,p.techs.all), p.time_steps] >= 0  # Power from technology t used to charge storage system b [kW]
 		dvDischargeFromStorage[p.s.storage.types.all, p.time_steps] >= 0 # Power discharged from storage system b [kW]
 		dvGridToStorage[p.s.storage.types.elec, p.time_steps] >= 0 # Electrical power delivered to storage by the grid [kW]
+		dvStorageToGrid[p.s.storage.types.elec, p.time_steps] >= 0 # TODO, add: "p.StorageSalesTiers" as well? export of energy from storage to the grid
 		dvStoredEnergy[p.s.storage.types.all, 0:p.time_steps[end]] >= 0  # State of charge of storage system b
 		dvStoragePower[p.s.storage.types.all] >= 0   # Power capacity of storage system b [kW]
 		dvStorageEnergy[p.s.storage.types.all] >= 0   # Energy capacity of storage system b [kWh]
@@ -594,6 +598,14 @@ function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
 		dvPeakDemandMonth[p.months, 1:p.s.electric_tariff.n_monthly_demand_tiers] >= 0  # Peak electrical power demand during month m [kW]
 		MinChargeAdder >= 0
         binGHP[p.ghp_options], Bin  # Can be <= 1 if require_ghp_purchase=0, and is ==1 if require_ghp_purchase=1
+				
+	end
+
+	for b in p.s.storage.types.elec
+		if !(p.s.storage.attr[b].allow_simultaneous_charge_discharge)
+			@warn "Adding binary variable to prevent simultaneous battery charge/discharge. Some solvers are very slow with integer variables."
+			@variable(m, binBattCharging[p.time_steps], Bin) # Binary for battery charging (vs discharging)
+		end
 	end
 
 	if !isempty(p.techs.gen)  # Problem becomes a MILP
@@ -612,7 +624,7 @@ function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
     end
 
 	if !(p.s.electric_utility.allow_simultaneous_export_import) & !isempty(p.s.electric_tariff.export_bins)
-		@warn "Adding binary variable to prevent simultaneous grid import/export. Some solvers are very slow with integer variables"
+		@warn "Adding binary variable to prevent simultaneous grid import/export. Some solvers are very slow with integer variables."
 		@variable(m, binNoGridPurchases[p.time_steps], Bin)
 	end
 
@@ -644,7 +656,7 @@ function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
     end
 
 	if !isempty(p.s.electric_utility.outage_durations) # add dvUnserved Load if there is at least one outage
-		@warn "Adding binary variable to model outages. Some solvers are very slow with integer variables"
+		@warn "Adding binary variable to model outages. Some solvers are very slow with integer variables."
 		max_outage_duration = maximum(p.s.electric_utility.outage_durations)
 		outage_time_steps = p.s.electric_utility.outage_time_steps
 		tZeros = p.s.electric_utility.outage_start_time_steps
diff --git a/src/core/reopt_multinode.jl b/src/core/reopt_multinode.jl
index 240956f6e..7693e79b8 100644
--- a/src/core/reopt_multinode.jl
+++ b/src/core/reopt_multinode.jl
@@ -2,7 +2,6 @@
 
 
 function add_variables!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}}) where T <: AbstractScenario
-	
 	dvs_idx_on_techs = String[
 		"dvSize",
 		"dvPurchaseSize",
@@ -16,7 +15,8 @@ function add_variables!(m::JuMP.AbstractModel, ps::AbstractVector{REoptInputs{T}
 		"dvStorageEnergy",
 	]
 	dvs_idx_on_storagetypes_time_steps = String[
-		"dvDischargeFromStorage"
+		"dvDischargeFromStorage",
+		"dvStorageToGrid"
 	]
 	for p in ps
 		_n = string("_", p.s.site.node)
diff --git a/src/mpc/model.jl b/src/mpc/model.jl
index dc8145033..7eabc8e82 100644
--- a/src/mpc/model.jl
+++ b/src/mpc/model.jl
@@ -99,6 +99,7 @@ function build_mpc!(m::JuMP.AbstractModel, p::MPCInputs)
 			@constraint(m, [ts in p.time_steps], m[:dvDischargeFromStorage][b, ts] == 0)
 			if b in p.s.storage.types.elec
 				@constraint(m, [ts in p.time_steps], m[:dvGridToStorage][b, ts] == 0)
+				@constraint(m, [ts in p.time_steps], m[:dvStorageToGrid][b, ts] == 0)
 			end
 		else
 			add_general_storage_dispatch_constraints(m, p, b)
@@ -230,6 +231,7 @@ function add_variables!(m::JuMP.AbstractModel, p::MPCInputs)
 		dvCurtail[p.techs.all, p.time_steps] >= 0  # [kW]
 		dvProductionToStorage[p.s.storage.types.all, p.techs.all, p.time_steps] >= 0  # Power from technology t used to charge storage system b [kW]
 		dvDischargeFromStorage[p.s.storage.types.all, p.time_steps] >= 0 # Power discharged from storage system b [kW]
+		dvStorageToGrid[p.s.storage.types.elec, p.time_steps] >= 0 # TODO, add: "p.StorageSalesTiers" as well? export of energy from storage to the grid
 		dvGridToStorage[p.s.storage.types.elec, p.time_steps] >= 0 # Electrical power delivered to storage by the grid [kW]
 		dvStoredEnergy[p.s.storage.types.all, 0:p.time_steps[end]] >= 0  # State of charge of storage system b
 		dvStoragePower[p.s.storage.types.all] >= 0   # Power capacity of storage system b [kW]
diff --git a/src/mpc/model_multinode.jl b/src/mpc/model_multinode.jl
index ab71d5765..03d1677d1 100644
--- a/src/mpc/model_multinode.jl
+++ b/src/mpc/model_multinode.jl
@@ -152,7 +152,8 @@ function add_variables!(m::JuMP.AbstractModel, ps::AbstractVector{MPCInputs})
 		"dvRatedProduction",
 	]
 	dvs_idx_on_storagetypes_time_steps = String[
-		"dvDischargeFromStorage"
+		"dvDischargeFromStorage",
+		"dvStorageToGrid"
 	]
 	for p in ps
 		_n = string("_", p.s.node)
diff --git a/src/mpc/structs.jl b/src/mpc/structs.jl
index 340bb07c0..ddc8cd665 100644
--- a/src/mpc/structs.jl
+++ b/src/mpc/structs.jl
@@ -226,7 +226,9 @@ Base.@kwdef struct MPCElectricStorage < AbstractElectricStorage
     soc_min_fraction::Float64 = 0.2
     soc_init_fraction::Float64 = 0.5
     can_grid_charge::Bool = true
+    can_export_to_grid::Bool = false
     grid_charge_efficiency::Float64 = 0.96 * 0.975^2
+    allow_simultaneous_charge_discharge::Bool = true
 end
 ```
 """
@@ -238,10 +240,13 @@ Base.@kwdef struct MPCElectricStorage <: AbstractElectricStorage
     soc_min_fraction::Float64 = 0.2
     soc_init_fraction::Float64 = 0.5
     can_grid_charge::Bool = true
+    can_export_to_grid::Bool = false
     grid_charge_efficiency::Float64 = 0.96 * 0.975^2
     max_kw::Float64 = size_kw
     max_kwh::Float64 = size_kwh
     minimum_avg_soc_fraction::Float64 = 0.0
+    self_discharge_fraction_per_timestep::Float64 = 0.0
+    allow_simultaneous_charge_discharge::Bool = true
 end
 
 
diff --git a/src/results/absorption_chiller.jl b/src/results/absorption_chiller.jl
index 86f424b99..2e974261f 100644
--- a/src/results/absorption_chiller.jl
+++ b/src/results/absorption_chiller.jl
@@ -25,6 +25,7 @@ function add_absorption_chiller_results(m::JuMP.AbstractModel, p::REoptInputs, d
 
 	r["size_kw"] = value(sum(m[:dvSize][t] for t in p.techs.absorption_chiller))
 	r["size_ton"] = r["size_kw"] / KWH_THERMAL_PER_TONHOUR
+	r["initial_capital_cost"] = round(value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.absorption_chiller)) * p.s.absorption_chiller.installed_cost_per_kw, digits=3)
 	@expression(m, ABSORPCHLtoTESKW[ts in p.time_steps],
 		sum(m[:dvProductionToStorage][b,t,ts] for b in p.s.storage.types.cold, t in p.techs.absorption_chiller))
 	r["thermal_to_storage_series_ton"] = round.(value.(ABSORPCHLtoTESKW) ./ KWH_THERMAL_PER_TONHOUR, digits=5)
diff --git a/src/results/boiler.jl b/src/results/boiler.jl
index cacdfb695..b4cf002b9 100644
--- a/src/results/boiler.jl
+++ b/src/results/boiler.jl
@@ -22,6 +22,7 @@
 function add_boiler_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
     r = Dict{String, Any}()
     r["size_mmbtu_per_hour"] = round(value(m[Symbol("dvSize"*_n)]["Boiler"]) / KWH_PER_MMBTU, digits=3)
+	r["initial_capital_cost"] = round(value(m[Symbol("dvSize"*_n)]["Boiler"]) * p.s.boiler.installed_cost_per_kw, digits=3)
 	r["fuel_consumption_series_mmbtu_per_hour"] = 
         round.(value.(m[:dvFuelUsage]["Boiler", ts] for ts in p.time_steps) / KWH_PER_MMBTU, digits=3)
     r["annual_fuel_consumption_mmbtu"] = round(sum(r["fuel_consumption_series_mmbtu_per_hour"]), digits=3)
diff --git a/src/results/chp.jl b/src/results/chp.jl
index ccd5aa5f0..21be26d47 100644
--- a/src/results/chp.jl
+++ b/src/results/chp.jl
@@ -30,7 +30,8 @@ function add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 	# Note: the node number is an empty string if evaluating a single `Site`.
     r = Dict{String, Any}()
 	r["size_kw"] = value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.chp))
-    r["size_supplemental_firing_kw"] = value(sum(m[Symbol("dvSupplementaryFiringSize"*_n)][t] for t in p.techs.chp))
+	#TODO: add initial_capital_cost to results (need to handle installed_cost_per_kw being a vector when cost curve is used)
+	r["size_supplemental_firing_kw"] = value(sum(m[Symbol("dvSupplementaryFiringSize"*_n)][t] for t in p.techs.chp))
 	@expression(m, CHPFuelUsedKWH, sum(m[Symbol("dvFuelUsage"*_n)][t, ts] for t in p.techs.chp, ts in p.time_steps))
 	r["annual_fuel_consumption_mmbtu"] = round(value(CHPFuelUsedKWH) / KWH_PER_MMBTU, digits=3)
 	@expression(m, Year1CHPElecProd,
diff --git a/src/results/electric_heater.jl b/src/results/electric_heater.jl
index af6e4bad2..4bac8aeb5 100644
--- a/src/results/electric_heater.jl
+++ b/src/results/electric_heater.jl
@@ -20,6 +20,7 @@
 function add_electric_heater_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
     r = Dict{String, Any}()
     r["size_mmbtu_per_hour"] = round(value(m[Symbol("dvSize"*_n)]["ElectricHeater"]) / KWH_PER_MMBTU, digits=3)
+    r["initial_capital_cost"] = round(value(m[Symbol("dvSize"*_n)]["ElectricHeater"]) * p.s.electric_heater.installed_cost_per_kw, digits=3)
     @expression(m, ElectricHeaterElectricConsumptionSeries[ts in p.time_steps],
         p.hours_per_time_step * sum(m[:dvHeatingProduction][t,q,ts] / p.heating_cop[t] 
         for q in p.heating_loads, t in p.techs.electric_heater))
diff --git a/src/results/electric_storage.jl b/src/results/electric_storage.jl
index 4c2b1483c..95a8ef783 100644
--- a/src/results/electric_storage.jl
+++ b/src/results/electric_storage.jl
@@ -34,6 +34,14 @@ function add_electric_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::
         r["initial_capital_cost"] = r["size_kwh"] * p.s.storage.attr[b].installed_cost_per_kwh +
             r["size_kw"] * p.s.storage.attr[b].installed_cost_per_kw
 
+        StoragePerUnitOMCosts = p.third_party_factor * p.pwf_om * (p.s.storage.attr[b].om_cost_per_kw * m[Symbol("dvStoragePower"*_n)][b] +
+                                                                 p.s.storage.attr[b].om_cost_per_kwh * m[Symbol("dvStorageEnergy"*_n)][b])
+
+        r["lifecycle_om_cost_after_tax"] = round(value(StoragePerUnitOMCosts) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
+        r["lifecycle_om_cost_before_tax"] = round(value(StoragePerUnitOMCosts), digits=0)
+        r["year_one_om_cost_before_tax"] = round(value(StoragePerUnitOMCosts) / (p.pwf_om * p.third_party_factor), digits=0)
+        r["year_one_om_cost_after_tax"] = round(value(StoragePerUnitOMCosts) * (1 - p.s.financial.owner_tax_rate_fraction) / (p.pwf_om * p.third_party_factor), digits=0)
+            
         if p.s.storage.attr[b].model_degradation
             r["state_of_health"] = value.(m[:SOH]).data / value.(m[:dvStorageEnergy])["ElectricStorage"];
             r["maintenance_cost"] = value(m[:degr_cost])
@@ -43,9 +51,13 @@ function add_electric_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::
                 ))
             end
          end
+         # report the exported electricity from the battery:
+         r["storage_to_grid_series_kw"] = round.(value.(m[Symbol("dvStorageToGrid"*_n)][b, :]), digits = 3)
+
     else
         r["soc_series_fraction"] = []
         r["storage_to_load_series_kw"] = []
+        r["storage_to_grid_series_kw"] = []
     end
 
     d[b] = r
diff --git a/src/results/electric_tariff.jl b/src/results/electric_tariff.jl
index 8bc2c05a9..331c07a8e 100644
--- a/src/results/electric_tariff.jl
+++ b/src/results/electric_tariff.jl
@@ -36,7 +36,9 @@ function add_electric_tariff_results(m::JuMP.AbstractModel, p::REoptInputs, d::D
     r["year_one_min_charge_adder_before_tax"] = round(value(m[Symbol("MinChargeAdder"*_n)]) / p.pwf_e, digits=2)
                                 
     r["lifecycle_export_benefit_after_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=2)
+    r["lifecycle_export_benefit_before_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]), digits=2)
     r["year_one_export_benefit_before_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) / p.pwf_e, digits=0)
+    r["year_one_export_benefit_after_tax"] = -1 * round(value(m[Symbol("TotalExportBenefit"*_n)]) * (1 - p.s.financial.offtaker_tax_rate_fraction) / p.pwf_e, digits=0)
 
     r["lifecycle_coincident_peak_cost_after_tax"] = round(value(m[Symbol("TotalCPCharges"*_n)]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=2)
     r["year_one_coincident_peak_cost_before_tax"] = round(value(m[Symbol("TotalCPCharges"*_n)]) / p.pwf_e, digits=2)
diff --git a/src/results/generator.jl b/src/results/generator.jl
index ec9c2efb7..22e6ab975 100644
--- a/src/results/generator.jl
+++ b/src/results/generator.jl
@@ -32,6 +32,7 @@ function add_generator_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _
 			for t in p.techs.gen, ts in p.time_steps)
 	)
 	r["size_kw"] = round(value(sum(m[:dvSize][t] for t in p.techs.gen)), digits=2)
+	r["initial_capital_cost"] = round(value(sum(m[:dvSize][t] for t in p.techs.gen)) * p.s.generator.installed_cost_per_kw, digits=3)
 	r["lifecycle_fixed_om_cost_after_tax"] = round(value(GenPerUnitSizeOMCosts) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
 	r["lifecycle_variable_om_cost_after_tax"] = round(value(m[:TotalPerUnitProdOMCosts]) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
 	r["lifecycle_fuel_cost_after_tax"] = round(value(m[:TotalGenFuelCosts]) * (1 - p.s.financial.offtaker_tax_rate_fraction), digits=2)
diff --git a/src/results/pv.jl b/src/results/pv.jl
index 439f5e102..514ab0934 100644
--- a/src/results/pv.jl
+++ b/src/results/pv.jl
@@ -66,6 +66,8 @@ function add_pv_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 		PVPerUnitSizeOMCosts = p.third_party_factor * p.om_cost_per_kw[t] * p.pwf_om * m[Symbol("dvSize"*_n)][t]
 		r["lifecycle_om_cost_after_tax"] = round(value(PVPerUnitSizeOMCosts) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
         r["lcoe_per_kwh"] = calculate_lcoe(p, r, get_pv_by_name(t, p.s.pvs))
+        r["year_one_om_cost_before_tax"] = round(value(PVPerUnitSizeOMCosts) / p.pwf_om, digits=0)
+        r["initial_capital_cost"] = round(value(m[Symbol("dvSize"*_n)][t]) * get_pv_by_name(t, p.s.pvs).installed_cost_per_kw, digits=3)
         d[t] = r
 	end
     nothing
diff --git a/src/results/steam_turbine.jl b/src/results/steam_turbine.jl
index e25aa37a2..d0b6d283d 100644
--- a/src/results/steam_turbine.jl
+++ b/src/results/steam_turbine.jl
@@ -25,6 +25,7 @@ function add_steam_turbine_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dic
     r = Dict{String, Any}()
 
 	r["size_kw"] = round(value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.steam_turbine)), digits=3)
+	r["initial_capital_cost"] = round(value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.steam_turbine)) * p.s.steam_turbine.installed_cost_per_kw, digits=3)
     @expression(m, Year1SteamTurbineThermalConsumptionKWH,
 		p.hours_per_time_step * sum(m[Symbol("dvThermalToSteamTurbine"*_n)][tst,q,ts] for tst in p.techs.can_supply_steam_turbine, q in p.heating_loads, ts in p.time_steps))
     r["annual_thermal_consumption_mmbtu"] = round(value(Year1SteamTurbineThermalConsumptionKWH) / KWH_PER_MMBTU, digits=5)
@@ -99,7 +100,6 @@ function add_steam_turbine_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dic
         @expression(m, SteamTurbineToProcessHeatKW[ts in p.time_steps], 0.0)
     end
     r["thermal_to_process_heat_load_series_mmbtu_per_hour"] = round.(value.(SteamTurbineToProcessHeatKW ./ KWH_PER_MMBTU), digits=5)
-
 	
 	d["SteamTurbine"] = r
 	nothing
diff --git a/src/results/thermal_storage.jl b/src/results/thermal_storage.jl
index 48e9f607c..602269cb4 100644
--- a/src/results/thermal_storage.jl
+++ b/src/results/thermal_storage.jl
@@ -20,6 +20,7 @@ function add_hot_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict,
     r = Dict{String, Any}()
     size_kwh = round(value(m[Symbol("dvStorageEnergy"*_n)][b]), digits=3)
     r["size_gal"] = round(size_kwh / kwh_per_gal, digits=0)
+	r["initial_capital_cost"] = round(value(m[Symbol("dvStorageEnergy"*_n)][b]) * p.s.storage.attr[b].installed_cost_per_kwh, digits=3)
 
     if size_kwh != 0
     	soc = (m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps)
@@ -93,12 +94,13 @@ function add_cold_storage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict
     Note: the node number is an empty string if evaluating a single `Site`.
     =#
 
-    kwh_per_gal = get_kwh_per_gal(p.s.storage.attr["ColdThermalStorage"].hot_water_temp_degF,
-                                    p.s.storage.attr["ColdThermalStorage"].cool_water_temp_degF)
+    kwh_per_gal = get_kwh_per_gal(p.s.storage.attr[b].hot_water_temp_degF,
+                                    p.s.storage.attr[b].cool_water_temp_degF)
     
     r = Dict{String, Any}()
     size_kwh = round(value(m[Symbol("dvStorageEnergy"*_n)][b]), digits=3)
     r["size_gal"] = round(size_kwh / kwh_per_gal, digits=0)
+	r["initial_capital_cost"] = round(value(m[Symbol("dvStorageEnergy"*_n)][b]) * p.s.storage.attr[b].installed_cost_per_kwh, digits=3)
 
     if size_kwh != 0
     	soc = (m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps)
diff --git a/src/results/wind.jl b/src/results/wind.jl
index fe16dc065..39a713f52 100644
--- a/src/results/wind.jl
+++ b/src/results/wind.jl
@@ -26,7 +26,7 @@ function add_wind_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 	r["production_factor_series"] = Vector(p.production_factor[t, :])
 	per_unit_size_om = @expression(m, p.third_party_factor * p.pwf_om * m[:dvSize][t] * p.om_cost_per_kw[t])
 
-	r["size_kw"] = round(value(m[:dvSize][t]), digits=2)
+	r["size_kw"] = round(value(m[Symbol("dvSize"*_n)][t]), digits=2)
 	r["lifecycle_om_cost_after_tax"] = round(value(per_unit_size_om) * (1 - p.s.financial.owner_tax_rate_fraction), digits=0)
 	r["year_one_om_cost_before_tax"] = round(value(per_unit_size_om) / (p.pwf_om * p.third_party_factor), digits=0)
 
@@ -66,6 +66,7 @@ function add_wind_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
 	r["annual_energy_produced_kwh"] = round(value(AvgWindProd), digits=0)
 
     r["lcoe_per_kwh"] = calculate_lcoe(p, r, p.s.wind)
+	r["initial_capital_cost"] = round(value(m[Symbol("dvSize"*_n)][t]) * p.s.wind.installed_cost_per_kw, digits=3)
 	d[t] = r
     nothing
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 4aee1d809..bed58fa6d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -606,6 +606,38 @@ else  # run HiGHS tests
             end
         end
 
+        @testset "Electric Storage O&M and Self-Discharge" begin
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            data = JSON.parsefile("./scenarios/storage_om.json")
+
+            data["ElectricStorage"]["om_cost_per_kw"] = 10
+            data["ElectricStorage"]["om_cost_per_kwh"] = 0
+
+            inputs = REoptInputs(data)
+            results1 = run_reopt(model, inputs)
+
+            @test results1["ElectricStorage"]["year_one_om_cost_before_tax"] ≈ 400 atol=1
+            @test results1["ElectricStorage"]["lifecycle_om_cost_after_tax"] ≈ 6272 atol=1 
+
+            data["ElectricStorage"]["om_cost_per_kw"] = 10
+            data["ElectricStorage"]["om_cost_per_kwh"] = 5
+
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            inputs = REoptInputs(data)
+            results2 = run_reopt(model, inputs)
+
+            @test results2["ElectricStorage"]["year_one_om_cost_before_tax"] ≈ 800 atol=1
+            @test results2["ElectricStorage"]["lifecycle_om_cost_after_tax"] ≈ 12543 atol=1 
+            @test results2["Financial"]["lcc"] ≈ results1["Financial"]["lcc"] + 12543 - 6272 atol=1
+            
+            data["ElectricStorage"]["self_discharge_fraction_per_timestep"] = 0.0025/24
+            model = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
+            inputs = REoptInputs(data)
+            results3 = run_reopt(model, inputs)
+            @test sum(results2["ElectricStorage"]["storage_to_load_series_kw"]) - 
+                sum(results3["ElectricStorage"]["storage_to_load_series_kw"]) ≈ 15.382 atol=0.01 
+        end
+        
         @testset "Heating loads and addressable load fraction" begin
             # Default LargeOffice CRB with SpaceHeatingLoad and DomesticHotWaterLoad are served by ExistingBoiler
             m = Model(optimizer_with_attributes(HiGHS.Optimizer, "output_flag" => false, "log_to_console" => false))
diff --git a/test/scenarios/storage_om.json b/test/scenarios/storage_om.json
new file mode 100644
index 000000000..c188b51b8
--- /dev/null
+++ b/test/scenarios/storage_om.json
@@ -0,0 +1,46 @@
+{
+    "Site": {
+        "longitude": -118.1164613,
+        "latitude": 34.5794343,
+        "roof_squarefeet": 5000.0,
+        "land_acres": 1.0
+    },
+    "ElectricLoad": {
+        "doe_reference_name": "RetailStore",
+        "annual_kwh": 10000000.0,
+        "year": 2017
+    },
+    "ElectricStorage": {
+        "min_kw": 40,
+        "max_kw": 40,
+        "min_kwh": 80,
+        "max_kwh": 80,   
+        "can_grid_charge": false,
+        "replace_cost_per_kw": 460.0,
+        "replace_cost_per_kwh": 230.0,
+        "installed_cost_per_kw": 1000.0,
+        "installed_cost_per_kwh": 500.0,
+        "total_itc_fraction": 0.0,
+        "macrs_option_years": 0,
+        "macrs_bonus_fraction": 0,
+        "soc_init_fraction": 1.0,
+        "soc_min_fraction": 0.0,
+        "can_export_to_grid": true
+    },
+    "ElectricTariff": {
+        "urdb_label": "5ed6c1a15457a3367add15ae"
+    },
+    "PV": {
+        "min_kw": 0.0,
+        "max_kw": 0.0
+    },
+    "Financial": {
+        "elec_cost_escalation_rate_fraction": 0.026,
+        "offtaker_discount_rate_fraction": 0.05,
+        "owner_discount_rate_fraction": 0.05,
+        "analysis_years": 20,
+        "offtaker_tax_rate_fraction": 0,
+        "owner_tax_rate_fraction": 0,
+        "om_cost_escalation_rate_fraction": 0.025
+    }
+}