degradation model cleanup

Moving old degradation branch changes to this "cleaned up" branch

src/constraints/battery_degradation.jl

@@ -7,6 +7,7 @@ function add_degradation_variables(m, p)
     @variable(m, Eplus_sum[days] >= 0)
     @variable(m, Eminus_sum[days] >= 0)
     @variable(m, EFC[days] >= 0)
+    @variable(m, SOH[days])
 end
 
 
@@ -47,33 +48,35 @@ NOTE the average SOC and EFC variables are in absolute units. For example, the S
     at the battery capacity in kWh.
 """
 function add_degradation(m, p; b="ElectricStorage")
+
+    # Indices
     days = 1:365*p.s.financial.analysis_years
+    months = 1:p.s.financial.analysis_years*12
+
     strategy = p.s.storage.attr[b].degradation.maintenance_strategy
 
     if isempty(p.s.storage.attr[b].degradation.maintenance_cost_per_kwh)
-        function pwf(day::Int)
+        # Correctly account for discount rate and install cost declination rate for days over analysis period
+        function pwf_bess_replacements(day::Int)
             (1-p.s.storage.attr[b].degradation.installed_cost_per_kwh_declination_rate)^(day/365) / 
             (1+p.s.financial.owner_discount_rate_fraction)^(day/365)
         end
-        # for the augmentation strategy the maintenance cost curve (function of time) starts at 
-        # 80% of the installed cost since we are not replacing the entire battery
-        f = strategy == "augmentation" ? 0.8 : 1.0
-        p.s.storage.attr[b].degradation.maintenance_cost_per_kwh = [ f * 
-            p.s.storage.attr[b].installed_cost_per_kwh * pwf(d) for d in days[1:end-1]
+        p.s.storage.attr[b].degradation.maintenance_cost_per_kwh = [ 
+            p.s.storage.attr[b].installed_cost_per_kwh * pwf_bess_replacements(d) for d in days[1:end-1]
         ]
     end
 
-    @assert(length(p.s.storage.attr[b].degradation.maintenance_cost_per_kwh) == length(days) - 1,
-        "The degradation maintenance_cost_per_kwh must have a length of $(length(days)-1)."
-    )
-
-    @variable(m, SOH[days])
+    # Under augmentation scenario, each day's battery augmentation cost is calculated using day-1 value from maintenance_cost_per_kwh vector
+    # Therefore, on last day, day-1's maintenance cost is utilized.
+    if length(p.s.storage.attr[b].degradation.maintenance_cost_per_kwh) != length(days) - 1
+        throw(@error("The degradation maintenance_cost_per_kwh must have a length of $(length(days)-1)."))
+    end
 
     add_degradation_variables(m, p)
     constrain_degradation_variables(m, p, b=b)
 
     @constraint(m, [d in 2:days[end]],
-        SOH[d] == SOH[d-1] - p.hours_per_time_step * (
+        m[:SOH][d] == m[:SOH][d-1] - p.hours_per_time_step * (
             p.s.storage.attr[b].degradation.calendar_fade_coefficient * 
             p.s.storage.attr[b].degradation.time_exponent * 
             m[:Eavg][d-1] * d^(p.s.storage.attr[b].degradation.time_exponent-1) + 
@@ -82,7 +85,7 @@ function add_degradation(m, p; b="ElectricStorage")
     )
     # NOTE SOH can be negative
 
-    @constraint(m, SOH[1] == m[:dvStorageEnergy][b])
+    @constraint(m, m[:SOH][1] == m[:dvStorageEnergy][b])
     # NOTE SOH is _not_ normalized, and has units of kWh
 
     if strategy == "replacement"
@@ -100,29 +103,9 @@ function add_degradation(m, p; b="ElectricStorage")
         The first month that the battery is replaced is determined by d_0p8, which is the integer 
         number of days that the SOH is at least 80% of the purchased capacity.
         We define a binary for each month and only allow one month to be chosen.
-        =#
-
-        # define d_0p8
-        @warn "Adding binary and indicator constraints for 
-         ElectricStorage.degradation.maintenance_strategy = \"replacement\". 
-         Not all solvers support indicators and some are slow with integers."
-        # TODO import the latest battery degradation model in the degradation branch
-        @variable(m, soh_indicator[days], Bin)
-        @constraint(m, [d in days],
-            soh_indicator[d] => {SOH[d] >= 0.8*m[:dvStorageEnergy][b]}
-        )
-        @expression(m, d_0p8, sum(soh_indicator[d] for d in days))
 
-        # define binaries for the finding the month that battery must be replaced
-        months = 1:p.s.financial.analysis_years*12
-        @variable(m, bmth[months], Bin)
-        # can only pick one month (or no month if SOH is >= 80% in last day)
-        @constraint(m, sum(bmth[mth] for mth in months) == 1-soh_indicator[length(days)])
-        # the month picked is at most the month in which the SOH hits 80%
-        @constraint(m, sum(mth*bmth[mth] for mth in months) <= d_0p8 / 30.42)
-        # 30.42 is the average number of days in a month
-
-        #=
+            # maintenance_cost_per_kwh must have length == length(days) - 1, i.e. starts on day 2
+        
         number of replacments as function of d_0p8
          ^
          |
@@ -139,41 +122,77 @@ function add_degradation(m, p; b="ElectricStorage")
         
         The above curve is multiplied by the maintenance_cost_per_kwh to create the cost coefficients
         =#
-        c = zeros(length(months))  # initialize cost coefficients
-        N = 365*p.s.financial.analysis_years
+
+        @warn "Adding binary decision variables for 
+        ElectricStorage.degradation.maintenance_strategy = \"replacement\". 
+        Some solvers are slow with integers."
+
+        @variable(m, binSOHIndicator[months], Bin) # track SOH levels, should be 1 if SOH >= 80%, 0 otherwise
+        @variable(m, binSOHIndicatorChange[months], Bin) # track which month SOH indicator drops to < 80%
+        @variable(m, 0 <= dvSOHChangeTimesEnergy[months]) # track the kwh to be replaced in a replacement month
+
+        # the big M
+        if p.s.storage.attr[b].max_kwh == 1.0e6 || p.s.storage.attr[b].max_kwh == 0
+            # Under default max_kwh (i.e. not modeling large batteries) or max_kwh = 0
+            bigM_StorageEnergy = 24*maximum(p.s.electric_load.loads_kw)
+        else
+            # Select the larger value of maximum electric load or provided max_kwh size.
+            bigM_StorageEnergy = max(24*maximum(p.s.electric_load.loads_kw), p.s.storage.attr[b].max_kwh)
+        end
+
+        # HEALTHY: if binSOHIndicator is 1, then SOH >= 80%. If binSOHIndicator is 0 and SOH >= very negative number
+        @constraint(m, [mth in months], m[:SOH][Int(round(30.4167*mth))] >= 0.8*m[:dvStorageEnergy][b] - bigM_StorageEnergy * (1-binSOHIndicator[mth]))
+
+        # UNHEALTHY: if binSOHIndicator is 1, then SOH <= large number. If binSOHIndicator is 0 and SOH <= 80%
+        @constraint(m, [mth in months], m[:SOH][Int(round(30.4167*mth))] <= 0.8*m[:dvStorageEnergy][b] + bigM_StorageEnergy * (binSOHIndicator[mth]))
+
+        # binSOHIndicatorChange[mth] = binSOHIndicator[mth-1] - binSOHIndicator[mth].
+        # If replacement month is x, then binSOHIndicatorChange[x] = 1. All other binSOHIndicatorChange values will be 0s (either 1-1 or 0-0)
+        @constraint(m, m[:binSOHIndicatorChange][1] == 1 - m[:binSOHIndicator][1])
+        @constraint(m, [mth in 2:months[end]], m[:binSOHIndicatorChange][mth] == m[:binSOHIndicator][mth-1] - m[:binSOHIndicator][mth])
+
+        @expression(m, months_to_first_replacement, sum(m[:binSOHIndicator][mth] for mth in months))
+
+        # -> linearize the product of binSOHIndicatorChange & m[:dvStorageEnergy][b]
+        @constraint(m, [mth in months], m[:dvSOHChangeTimesEnergy][mth] >= m[:dvStorageEnergy][b] - bigM_StorageEnergy * (1 - m[:binSOHIndicatorChange][mth]))
+        @constraint(m, [mth in months], m[:dvSOHChangeTimesEnergy][mth] <= m[:dvStorageEnergy][b] + bigM_StorageEnergy * (1 - m[:binSOHIndicatorChange][mth]))
+        @constraint(m, [mth in months], m[:dvSOHChangeTimesEnergy][mth] <= bigM_StorageEnergy * m[:binSOHIndicatorChange][mth])
+
+        replacement_costs = zeros(length(months))  # initialize cost coefficients
+        residual_values = zeros(length(months))  # initialize cost coefficients for residual_value
+        N = 365*p.s.financial.analysis_years # number of days
+
         for mth in months
-            day = Int(round((mth-1)*30.42 + 15, digits=0))
-            c[mth] = p.s.storage.attr[b].degradation.maintenance_cost_per_kwh[day] *
-                ceil(N/day - 1)
+            day = Int(round((mth-1)*30.4167 + 15, digits=0))
+            batt_replace_count = Int(ceil(N/day - 1)) # number of battery replacements in analysis period if they periodically happened on "day"
+            maint_cost = sum(p.s.storage.attr[b].degradation.maintenance_cost_per_kwh[day*i] for i in 1:batt_replace_count)
+            replacement_costs[mth] = maint_cost
+
+            residual_factor = 1 - (p.s.financial.analysis_years*12/mth - floor(p.s.financial.analysis_years*12/mth))
+            residual_value = p.s.storage.attr[b].degradation.maintenance_cost_per_kwh[end]*residual_factor
+            residual_values[mth] = residual_value
         end
 
-        # linearize the product of bmth & m[:dvStorageEnergy][b]
-        M = p.s.storage.attr[b].max_kwh  # the big M
-        @variable(m, 0 <= bmth_BkWh[months])
-        @constraint(m, [mth in months], bmth_BkWh[mth] <= m[:dvStorageEnergy][b])
-        @constraint(m, [mth in months], bmth_BkWh[mth] <= M * bmth[mth])
-        @constraint(m, [mth in months], bmth_BkWh[mth] >= m[:dvStorageEnergy][b] - M*(1-bmth[mth]))
+        # create replacement cost expression for objective
+        @expression(m, degr_cost, sum(replacement_costs[mth] * m[:dvSOHChangeTimesEnergy][mth] for mth in months))
 
-        # add replacment cost to objective
-        @expression(m, degr_cost,
-            sum(c[mth] * bmth_BkWh[mth] for mth in months)
-        )
+        # create residual value expression for objective
+        @expression(m, residual_value, sum(residual_values[mth] * m[:dvSOHChangeTimesEnergy][mth] for mth in months))
 
     elseif strategy == "augmentation"
 
         @expression(m, degr_cost,
             sum(
-                p.s.storage.attr[b].degradation.maintenance_cost_per_kwh[d-1] * (SOH[d-1] - SOH[d])
+                p.s.storage.attr[b].degradation.maintenance_cost_per_kwh[d-1] * (m[:SOH][d-1] - m[:SOH][d])
                 for d in days[2:end]
             )
         )
-        # add augmentation cost to objective
-        # maintenance_cost_per_kwh must have length == length(days) - 1, i.e. starts on day 2
+
+        # No lifetime based residual value assigned to battery under the augmentation strategy
+        @expression(m, residual_value, 0.0)
     else
         throw(@error("Battery maintenance strategy $strategy is not supported. Choose from augmentation and replacement."))
     end
-
-    @objective(m, Min, m[:Costs] + m[:degr_cost])
 
     # NOTE adding to Costs expression does not modify the objective function
  end

src/core/electric_utility.jl

@@ -8,7 +8,7 @@
     
     # Single Outage Modeling Inputs (Outage Modeling Option 1)
     outage_start_time_step::Int=0,  # for modeling a single outage, with critical load spliced into the baseline load ...
-    outage_end_time_step::Int=0,  # ... utiltity production_factor = 0 during the outage
+    outage_end_time_step::Int=0,  # ... utility production_factor = 0 during the outage
         
     # Multiple Outage Modeling Inputs (Outage Modeling Option 2): minimax the expected outage cost,
     # with max taken over outage start time, expectation taken over outage duration
@@ -115,7 +115,7 @@ struct ElectricUtility
     emissions_factor_SO2_decrease_fraction::Real
     emissions_factor_PM25_decrease_fraction::Real
     outage_start_time_step::Int  # for modeling a single outage, with critical load spliced into the baseline load ...
-    outage_end_time_step::Int  # ... utiltity production_factor = 0 during the outage
+    outage_end_time_step::Int  # ... utility production_factor = 0 during the outage
     allow_simultaneous_export_import::Bool  # if true the site has two meters (in effect)
     # next 5 variables below used for minimax the expected outage cost,
     # with max taken over outage start time, expectation taken over outage duration
@@ -147,7 +147,7 @@ struct ElectricUtility
         net_metering_limit_kw::Real = 0, # Upper limit on the total capacity of technologies that can participate in net metering agreement.
         interconnection_limit_kw::Real = 1.0e9,
         outage_start_time_step::Int=0,  # for modeling a single outage, with critical load spliced into the baseline load ...
-        outage_end_time_step::Int=0,  # ... utiltity production_factor = 0 during the outage
+        outage_end_time_step::Int=0,  # ... utility production_factor = 0 during the outage
         allow_simultaneous_export_import::Bool=true,  # if true the site has two meters (in effect)
         # next 5 variables below used for minimax the expected outage cost,
         # with max taken over outage start time, expectation taken over outage duration