Merge pull request #463 from NREL/develop

v0.48.2 Dec 2024

Author: adfarth

CHANGELOG.md

@@ -25,6 +25,20 @@ Classify the change according to the following categories:
     ### Deprecated
     ### Removed
 
+
+## v0.48.2
+### Added
+- Battery residual value if choosing replacement strategy for degradation
+- Add new **ElectricStorage** parameters **max_duration_hours** and **min_duration_hours** to bound the energy duration of battery storage
+### Changed
+- Revised the battery degradation model, refactoring some methods to increase model-building efficiency and reformulating indicator constraints as big-M constraints with smaller big-M's to reduce solve time.
+- Edited several documentation entries and docstrings for clarity.
+### Removed
+- 80% scaling of battery maintenance costs when using augmentation strategy
+### Fixed 
+- Fixed conditions for which a warning is presented indicating that the wholesale benefit threshold is met. 
+- When setting **thermal_production_series_mmbtu_per_hour** output in **ExistingBoiler**, sum over heating loads instead of time steps
+
 ## v0.48.1
 ### Changed
 - Replace all `1/p.s.settings.time_steps_per_hour` with `p.hours_per_time_step` for simplicity/consistency
@@ -48,6 +62,7 @@ Classify the change according to the following categories:
 - 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
 - In `src/core/existing_chiller.jl`, Added new atttribute **retire_in_optimal** to the **ExistingChiller** struct
 - Financial output **initial_capital_costs_after_incentives_without_macrs** which has "net year one" CapEx after incentives except for MACRS, which helps with users defining their own "simple payback period"
+
 ### Changed
 - Improve the full test suite reporting with a verbose summary table, and update the structure to reflect long-term open-source solver usage.
 - Removed MacOS from the runner list and just run with Windows OS, since MacOS commonly freezes and gets cancelled. We have not seen Windows OS pass while other OS's fail.
@@ -64,6 +79,7 @@ Classify the change according to the following categories:
 - An issue with setup_boiler_inputs in reopt_inputs.jl.
 - Fuel costs in proforma.jl were not consistent with the optimization costs, so that was corrected so that they are only added to the offtaker cashflows and not the owner/developer cashflows for third party.
 
+
 ## v0.47.2
 ### Fixed
 - Increased the big-M bound on maximum net metering benefit to prevent artificially low export benefits.

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.3"
 manifest_format = "2.0"
-project_hash = "2e3f73051e60a5c2ffca2d998f6043d51c7f6c2a"
+project_hash = "6787ce6711da7a7bf012a2defdfff479e23ec7c6"
 
 [[deps.AbstractFFTs]]
 deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
@@ -57,9 +57,9 @@ version = "1.21.4+0"
 
 [[deps.Bzip2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
+git-tree-sha1 = "8873e196c2eb87962a2048b3b8e08946535864a1"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
-version = "1.0.8+0"
+version = "1.0.8+2"
 
 [[deps.CEnum]]
 git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90"
@@ -513,9 +513,9 @@ version = "1.9.4+0"
 
 [[deps.MPICH_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "8a5b4d2220377d1ece13f49438d71ad20cf1ba83"
+git-tree-sha1 = "2ee75365ca243c1a39d467e35ffd3d4d32eef11e"
 uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4"
-version = "4.1.2+0"
+version = "4.1.2+1"
 
 [[deps.MPIPreferences]]
 deps = ["Libdl", "Preferences"]
@@ -525,9 +525,9 @@ version = "0.1.9"
 
 [[deps.MPItrampoline_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "6979eccb6a9edbbb62681e158443e79ecc0d056a"
+git-tree-sha1 = "8eeb3c73bbc0ca203d0dc8dad4008350bbe5797b"
 uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748"
-version = "5.3.1+0"
+version = "5.3.1+1"
 
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
@@ -569,9 +569,9 @@ version = "1.4.1"
 
 [[deps.MicrosoftMPI_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "a7023883872e52bc29bcaac74f19adf39347d2d5"
+git-tree-sha1 = "bc95bf4149bf535c09602e3acdf950d9b4376227"
 uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf"
-version = "10.1.4+0"
+version = "10.1.4+3"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
@@ -942,9 +942,9 @@ version = "100.700.100+0"
 
 [[deps.libpng_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
+git-tree-sha1 = "f7c281e9c61905521993a987d38b5ab1d4b53bef"
 uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.38+0"
+version = "1.6.38+1"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]

Project.toml

@@ -1,7 +1,7 @@
 name = "REopt"
 uuid = "d36ad4e8-d74a-4f7a-ace1-eaea049febf6"
 authors = ["Nick Laws", "Hallie Dunham <[email protected]>", "Bill Becker <[email protected]>", "Bhavesh Rathod <[email protected]>", "Alex Zolan <[email protected]>", "Amanda Farthing <[email protected]>"]
-version = "0.48.1"
+version = "0.48.2"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

src/constraints/battery_degradation.jl

@@ -7,36 +7,39 @@ 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
 
 
 function constrain_degradation_variables(m, p; b="ElectricStorage")
     days = 1:365*p.s.financial.analysis_years
     ts_per_day = 24 / p.hours_per_time_step
     ts_per_year = ts_per_day * 365
+    ts0 = Dict()
+    tsF = Dict()
     for d in days
-        ts0 = Int((ts_per_day * (d - 1) + 1) % ts_per_year)
-        tsF = Int(ts_per_day * d % ts_per_year)
-        if tsF == 0
-            tsF = Int(ts_per_day * 365)
+        ts0[d] = Int((ts_per_day * (d - 1) + 1) % ts_per_year)
+        tsF[d] = Int(ts_per_day * d % ts_per_year)
+        if tsF[d] == 0
+            tsF[d] = Int(ts_per_day * 365)
         end
-        @constraint(m, 
-            m[:Eavg][d] == 1/ts_per_day * sum(m[:dvStoredEnergy][b, ts] for ts in ts0:tsF)
-        )
-        @constraint(m,
-            m[:Eplus_sum][d] == 
-                p.hours_per_time_step * (
-                    sum(m[:dvProductionToStorage][b, t, ts] for t in p.techs.elec, ts in ts0:tsF) 
-                    + sum(m[:dvGridToStorage][b, ts] for ts in ts0:tsF)
-                )
-        )
-        @constraint(m,
-            m[:Eminus_sum][d] == p.hours_per_time_step * sum(m[:dvDischargeFromStorage][b, ts] for ts in ts0:tsF)
-        )
-        @constraint(m,
-            m[:EFC][d] == (m[:Eplus_sum][d] + m[:Eminus_sum][d]) / 2
-        )
     end
+    @constraint(m, [d in days],
+        m[:Eavg][d] == 1/ts_per_day * sum(m[:dvStoredEnergy][b, ts] for ts in ts0[d]:tsF[d])
+    )
+    @constraint(m, [d in days],
+        m[:Eplus_sum][d] == 
+            p.hours_per_time_step * (
+                sum(m[:dvProductionToStorage][b, t, ts] for t in p.techs.elec, ts in ts0[d]:tsF[d]) 
+                + sum(m[:dvGridToStorage][b, ts] for ts in ts0[d]:tsF[d])
+            )
+    )
+    @constraint(m, [d in days],
+        m[:Eminus_sum][d] == p.hours_per_time_step * sum(m[:dvDischargeFromStorage][b, ts] for ts in ts0[d]:tsF[d])
+    )
+    @constraint(m, [d in days],
+        m[:EFC][d] == (m[:Eplus_sum][d] + m[:Eminus_sum][d]) / 2
+    )
 end
 
 
@@ -47,33 +50,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 +87,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 +105,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 +124,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/constraints/storage_constraints.jl

@@ -21,6 +21,19 @@ function add_storage_size_constraints(m, p, b; _n="")
 	@constraint(m,
         m[Symbol("dvStoragePower"*_n)][b] <= p.s.storage.attr[b].max_kw
     )
+
+	# Constraint (4c)-3: Limit on ElectricStorage Energy Capacity based on Duration Hours
+    if p.s.storage.attr[b] isa ElectricStorage
+        @constraint(m,
+            m[Symbol("dvStorageEnergy"*_n)][b] <= m[Symbol("dvStoragePower"*_n)][b] * p.s.storage.attr[b].max_duration_hours
+        )
+
+        @constraint(m,
+            m[Symbol("dvStorageEnergy"*_n)][b] >= m[Symbol("dvStoragePower"*_n)][b] * p.s.storage.attr[b].min_duration_hours
+        )
+    end    
+
+
 end