Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

degradation model cleanup #443

Draft
wants to merge 3 commits into
base: develop
Choose a base branch
from
Draft

degradation model cleanup #443

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7f627ec
degradation model cleanup
rathod-b Sep 19, 2024
6b20ef4
Update tests
rathod-b Sep 20, 2024
82617aa
Update runtests.jl
rathod-b Sep 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
129 changes: 74 additions & 55 deletions src/constraints/battery_degradation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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) +
Expand All @@ -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"
Expand All @@ -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
^
|
Expand All @@ -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
Expand Down
6 changes: 3 additions & 3 deletions src/core/electric_utility.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading
Loading