Annualised

src/EnergyModelsInvestments.jl

@@ -24,6 +24,7 @@ include(joinpath("structures", "investment_data.jl"))
 include(joinpath("structures", "legacy_constructors.jl"))
 
 # Core structure of the code
+include("crf.jl")
 include("model.jl")
 include("utils.jl")
 

src/crf.jl

@@ -0,0 +1,26 @@
+"""
+    CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ)
+
+Computes the Capital Recovery Factor (CRF) for the investment data `inv_data` 
+in the strategic period `t_inv`, given the set of investment periods `𝒯ᴵⁿᵛ`. 
+
+The formula of the annuity-due factor, assuming the payments start immediately, in the same period as the
+investment. 
+
+The CRF is calculated based on the discount rate of `inv_data` and the remaining horizon of `t_inv` 
+within `𝒯ᴵⁿᵛ`. It represents the annualised payment factor used to recover the investment cost 
+over its economic lifetime.
+"""
+function CRF(inv_data::AbstractInvData, t_inv, 𝒯ᴵⁿᵛ)
+    disc_rate = get_discount_rate(inv_data)
+    remaining_horizon = remaining(t_inv, 𝒯ᴵⁿᵛ)
+    annuity = (disc_rate * (1 + disc_rate)^(remaining_horizon-1)) / ((1 + disc_rate)^(remaining_horizon) - 1)
+    return annuity
+end
+
+function set_period_annuity(inv_data::AbstractInvData, t_inv)
+    disc_rate = get_discount_rate(inv_data)
+    d = t_inv.duration
+
+    return ((1 + disc_rate)^d - 1)/(disc_rate * (1 + disc_rate)^(d-1))
+end

src/model.jl

@@ -232,7 +232,14 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat
 
     # The capacity has an unlimited lifetime, one investment at the beginning of t_inv
     capex_val = set_capex_value(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ)
-    @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == capex_val[t_inv])
+    if has_discount_rate(inv_data)
+        Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ)
+        annuity_capex = @expression(m, [t_inv ∈ 𝒯ᴵⁿᵛ], capex_val[t_inv] * CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ))
+        period_annuity_capex = @expression(m, [t_inv ∈ 𝒯ᴵⁿᵛ], annuity_capex[t_inv] * set_period_annuity(inv_data, t_inv))
+        @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]))
+    else
+        @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == capex_val[t_inv])
+    end
 
     # Fix the binary variable
     for t_inv ∈ 𝒯ᴵⁿᵛ
@@ -246,14 +253,28 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat
 
     # The capacity is limited to the end of the study. Reinvestments are included
     # No capacity removed as there are reinvestments according to the study length
-    capex_disc = StrategicProfile([
-        set_capex_discounter(
-            remaining(t_inv, 𝒯ᴵⁿᵛ),
-            lifetime(inv_data, t_inv), disc_rate
-        ) for t_inv ∈ 𝒯ᴵⁿᵛ
-    ])
     capex_val = set_capex_value(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ)
-    @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == capex_val[t_inv] * capex_disc[t_inv])
+    if has_discount_rate(inv_data)
+        capex_disc = StrategicProfile([
+            set_capex_discounter(
+                remaining(t_inv, 𝒯ᴵⁿᵛ),
+                lifetime(inv_data, t_inv), get_discount_rate(inv_data)
+            ) for t_inv ∈ 𝒯ᴵⁿᵛ
+        ])
+        Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ)
+
+        annuity_capex = @expression(m, [t_inv ∈ 𝒯ᴵⁿᵛ], capex_val[t_inv] * capex_disc[t_inv] * CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ))
+        period_annuity_capex = @expression(m, [t_inv ∈ 𝒯ᴵⁿᵛ], annuity_capex[t_inv] * set_period_annuity(inv_data, t_inv))
+        @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]))
+    else
+        capex_disc = StrategicProfile([
+            set_capex_discounter(
+                remaining(t_inv, 𝒯ᴵⁿᵛ),
+                lifetime(inv_data, t_inv), disc_rate
+            ) for t_inv ∈ 𝒯ᴵⁿᵛ
+        ])
+        @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == capex_val[t_inv] * capex_disc[t_inv])
+    end
 
     # Fix the binary variable
     for t_inv ∈ 𝒯ᴵⁿᵛ
@@ -269,14 +290,29 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat
     # The capacity is limited to the current sp. It has to be removed in the next sp.
     # The capacity removal variable is corresponding to the removal of the capacity at the
     # end of the investment period. Hence, we have to enforce `var_rem[t_inv] == var_add[t_inv]`
-    capex_disc = StrategicProfile([
-        set_capex_discounter(
-            duration_strat(t_inv),
-            lifetime(inv_data, t_inv), disc_rate
-        ) for t_inv ∈ 𝒯ᴵⁿᵛ
-    ])
     capex_val = set_capex_value(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ)
-    @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == capex_val[t_inv] * capex_disc[t_inv])
+    if has_discount_rate(inv_data)
+        capex_disc = StrategicProfile([
+            set_capex_discounter(
+                duration_strat(t_inv),
+                lifetime(inv_data, t_inv), get_discount_rate(inv_data)
+            ) for t_inv ∈ 𝒯ᴵⁿᵛ
+        ])
+        Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ)
+
+        annuity_capex = @expression(m, [t_inv ∈ 𝒯ᴵⁿᵛ], capex_val[t_inv] * capex_disc[t_inv] * CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ))
+        period_annuity_capex = @expression(m, [t_inv ∈ 𝒯ᴵⁿᵛ], annuity_capex[t_inv] * set_period_annuity(inv_data, t_inv))
+        @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]))
+    else
+        capex_disc = StrategicProfile([
+            set_capex_discounter(
+                duration_strat(t_inv),
+                lifetime(inv_data, t_inv), disc_rate
+            ) for t_inv ∈ 𝒯ᴵⁿᵛ
+        ])
+        @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == capex_val[t_inv] * capex_disc[t_inv])
+    end
+
     @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_rem[t_inv] == var_add[t_inv])
 end
 function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rate, ::RollingLife)
@@ -290,8 +326,11 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat
 
     # Initialize a dictionary for the removal of capacity
     rem_dict = Dict(t_inv => eltype(𝒯ᴵⁿᵛ)[] for t_inv ∈ 𝒯ᴵⁿᵛ)
+
+    period_annuity_capex_dict = Dict{TS.StrategicPeriod, Any}()
+    Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ)
 
-    for t_inv ∈ 𝒯ᴵⁿᵛ
+    for t_inv ∈ sort(collect(𝒯ᴵⁿᵛ))
         # Extract the values
         lifetime_val = lifetime(inv_data, t_inv)
 
@@ -303,48 +342,49 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat
         # If lifetime is shorter than the sp duration, we apply the method for PeriodLife
         # to account for the required reinvestments
         if lifetime_val < duration_strat(t_inv)
-            capex_disc = set_capex_discounter(duration_strat(t_inv), lifetime_val, disc_rate)
-            @constraint(m, var_capex[t_inv] == capex_val[t_inv] * capex_disc)
+            if has_discount_rate(inv_data)
+                capex_disc = set_capex_discounter(duration_strat(t_inv), lifetime_val, get_discount_rate(inv_data))
+
+                annuity_capex = @expression(m, capex_val[t_inv] * capex_disc * CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ))
+                period_annuity_capex = @expression(m, annuity_capex * set_period_annuity(inv_data, t_inv))
+                period_annuity_capex_dict[t_inv] = period_annuity_capex
+
+                @constraint(m, var_capex[t_inv] == sum(period_annuity_capex_dict[t] for t in Tᶜᵘᵐ[t_inv]))
+                #@constraint(m, var_capex[t_inv] == sum(capex_val[t] * capex_disc[t] * CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]))
+            else
+                capex_disc = set_capex_discounter(duration_strat(t_inv), lifetime_val, disc_rate)
+                @constraint(m, var_capex[t_inv] == capex_val[t_inv] * capex_disc)
+            end
             push!(rem_dict[t_inv_rem], t_inv)
 
         # If lifetime is equal to sp duration we only need to invest once and there is no
         # rest value. The invested capacity is removed at the end of the investment period
         elseif lifetime_val == duration_strat(t_inv)
-            @constraint(m, var_capex[t_inv] == capex_val[t_inv])
+            if has_discount_rate(inv_data)
+                annuity_capex = @expression(m, capex_val[t_inv] * CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ))
+                period_annuity_capex = @expression(m, annuity_capex * set_period_annuity(inv_data, t_inv))
+                period_annuity_capex_dict[t_inv] = period_annuity_capex
+
+                @constraint(m, var_capex[t_inv] == sum(period_annuity_capex_dict[t] for t in Tᶜᵘᵐ[t_inv]))
+            else
+                @constraint(m, var_capex[t_inv] == capex_val[t_inv])
+            end
             push!(rem_dict[t_inv_rem], t_inv)
 
         # If lifetime is longer than sp duration, the capacity can roll over to the next sp
         elseif lifetime_val > duration_strat(t_inv)
-            # Initialization of the the remaining lifetime
-            remaining_lifetime = lifetime_val
-            bool_lifetime = true
-
-            # Iteration to identify investment period in which the remaining lifetime is
-            # smaller than its duration
-            for sp ∈ 𝒯ᴵⁿᵛ
-                if sp ≥ t_inv
-                    if remaining_lifetime < duration_strat(sp)
-                        break
-                    end
-                    remaining_lifetime -= duration_strat(sp)
-                    t_inv_rem = sp
-                    if sp == last(𝒯ᴵⁿᵛ) && remaining_lifetime > 0
-                        bool_lifetime = false
-                    end
-                end
+            if has_discount_rate(inv_data)
+                r = get_discount_rate(inv_data)
+                capex_disc, rem_dict = get_capex_disc(lifetime_val, r, rem_dict, t_inv_rem, t_inv, 𝒯ᴵⁿᵛ)
+                annuity_capex = @expression(m, capex_val[t_inv] * capex_disc * CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ))
+                period_annuity_capex = @expression(m, annuity_capex * set_period_annuity(inv_data, t_inv))
+                period_annuity_capex_dict[t_inv] = period_annuity_capex
+
+                @constraint(m, var_capex[t_inv] == sum(period_annuity_capex_dict[t] for t in Tᶜᵘᵐ[t_inv]))
+            else
+                capex_disc, rem_dict = get_capex_disc(lifetime_val, disc_rate, rem_dict, t_inv_rem, t_inv, 𝒯ᴵⁿᵛ)
+                @constraint(m, var_capex[t_inv] == capex_val[t_inv] * capex_disc)
             end
-
-            # If the remaining life is larger than 0 at the end of the analysis horizon, we
-            # do not remove the capacity
-            bool_lifetime && push!(rem_dict[t_inv_rem], t_inv)
-
-            # Calculation of cost and rest value
-            capex_disc = (
-                1 -
-                (remaining_lifetime / lifetime_val) *
-                (1 + disc_rate)^(-(lifetime_val - remaining_lifetime))
-            )
-            @constraint(m, var_capex[t_inv] == capex_val[t_inv] * capex_disc)
         end
     end
     for (t_inv_rem, t_inv_vec) ∈ rem_dict

src/structures/investment_data.jl

@@ -30,14 +30,33 @@ struct NoStartInvData <: AbstractInvData
     max_inst::TimeProfile
     inv_mode::Investment
     life_mode::LifetimeMode
+    disc_rate::Union{Float64, Nothing} #Optional parameter
 end
 function NoStartInvData(
     capex_trans::TimeProfile,
     trans_max_inst::TimeProfile,
     inv_mode::Investment,
 )
 
-    return NoStartInvData(capex_trans, trans_max_inst, inv_mode, UnlimitedLife())
+    return NoStartInvData(capex_trans, trans_max_inst, inv_mode, UnlimitedLife(), nothing)
+end
+function NoStartInvData(
+    capex_trans::TimeProfile,
+    trans_max_inst::TimeProfile,
+    inv_mode::Investment,
+    disc_rate::Float64
+)
+
+    return NoStartInvData(capex_trans, trans_max_inst, inv_mode, UnlimitedLife(), disc_rate)
+end
+function NoStartInvData(
+    capex_trans::TimeProfile,
+    trans_max_inst::TimeProfile,
+    inv_mode::Investment,
+    life_mode::LifetimeMode
+)
+
+    return NoStartInvData(capex_trans, trans_max_inst, inv_mode, life_mode, nothing)
 end
 
 
@@ -58,14 +77,33 @@ struct StartInvData <: AbstractInvData
     initial::TimeProfile
     inv_mode::Investment
     life_mode::LifetimeMode
+    disc_rate::Union{Float64, Nothing} #Optional parameter
+end
+function StartInvData(
+    capex_trans::TimeProfile,
+    trans_max_inst::TimeProfile,
+    initial::TimeProfile,
+    inv_mode::Investment,
+)
+    return StartInvData(capex_trans, trans_max_inst, initial, inv_mode, UnlimitedLife(), nothing)
+end
+function StartInvData(
+    capex_trans::TimeProfile,
+    trans_max_inst::TimeProfile,
+    initial::TimeProfile,
+    inv_mode::Investment,
+    disc_rate::Float64
+)
+    return StartInvData(capex_trans, trans_max_inst, initial, inv_mode, UnlimitedLife(), disc_rate)
 end
 function StartInvData(
     capex_trans::TimeProfile,
     trans_max_inst::TimeProfile,
     initial::TimeProfile,
     inv_mode::Investment,
+    life_mode::LifetimeMode
 )
-    return StartInvData(capex_trans, trans_max_inst, initial, inv_mode, UnlimitedLife())
+    return StartInvData(capex_trans, trans_max_inst, initial, inv_mode, life_mode, nothing)
 end
 
 """
@@ -165,3 +203,22 @@ investment period `t_inv`.
 invest_capacity(inv_data::AbstractInvData) = invest_capacity(investment_mode(inv_data))
 invest_capacity(inv_data::AbstractInvData, t_inv) =
     invest_capacity(investment_mode(inv_data), t_inv)
+
+"""
+   get_discount_rate(inv_data::AbstractInvData)
+
+Returns the discount rate of the investment data `inv_data`. If the element has not an associated `disc_rate` it 
+returns `nothing`.
+"""
+get_discount_rate(inv_data::AbstractInvData) = inv_data.disc_rate
+
+"""
+    has_discount_rate(inv_data::AbstractInvData)
+
+Returns `true` if the investment data `inv_data` has an associated `disc_rate`. Otherwise, it returns `false`.
+"""
+
+function has_discount_rate(inv_data::AbstractInvData)
+    disc_rate = get_discount_rate(inv_data)
+    return isa(disc_rate, Float64)
+end

src/utils.jl

@@ -145,3 +145,49 @@ function set_capex_discounter(years, lifetime, disc_rate)
         ((N_inv * lifetime - years) / lifetime) * (1 + disc_rate)^(-years)
     return capex_disc
 end
+
+"""
+    get_cumulative_periods(𝒯::AbstractStratPers)
+
+Given a collection of strategic periods `𝒯`, returns a dictionary mapping each period `t` in `𝒯` 
+to a vector of all periods in `𝒯` up to and including `t`.
+This is used to retrieve the cumulative set of periods leading up to each strategic period.
+"""
+function get_cumulative_periods(𝒯::TS.AbstractStratPers)
+    chunks_t_inv = collect(collect(ts) for ts in chunk(Iterators.reverse(𝒯), 𝒯.ts.len))
+    chunks_t_inv_dict = Dict(t_inv => first(filter(c -> c[1] == t_inv, chunks_t_inv)) for t_inv in 𝒯)
+    return chunks_t_inv_dict
+end
+
+function get_capex_disc(lifetime_val, disc_rate, rem_dict, t_inv_rem, t_inv, 𝒯ᴵⁿᵛ)
+    # Initialization of the remaining lifetime
+    remaining_lifetime = lifetime_val
+    bool_lifetime = true
+
+    # Iteration to identify investment period in which the remaining lifetime is 
+    # smaller than its duration
+    for sp ∈ 𝒯ᴵⁿᵛ
+        if sp ≥ t_inv
+            if remaining_lifetime < duration_strat(sp)
+                break
+            end
+            remaining_lifetime -= duration_strat(sp)
+            t_inv_rem = sp
+            if sp == last(𝒯ᴵⁿᵛ) && remaining_lifetime > 0
+                bool_lifetime = false
+            end
+        end
+    end
+    # If the remaining life is larger than 0 at the end of the analysis horizon, we
+    # do not remove the capacity
+    bool_lifetime && push!(rem_dict[t_inv_rem], t_inv)
+
+     # Calculation of cost and rest value
+    capex_disc = (
+        1 -
+        (remaining_lifetime / lifetime_val) *
+        (1 + disc_rate)^(-(lifetime_val - remaining_lifetime))
+    )
+
+    return capex_disc, rem_dict
+end

test/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+EnergyModelsInvestments = "fca3f8eb-b383-437d-8e7b-aac76bb2004f"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

test/runtests.jl

@@ -21,4 +21,8 @@ include("utils.jl")
     @testset "Investments | Lifetime" begin
         include("test_lifetime.jl")
     end
+
+    @testset "Investments | Annualised" begin
+        include("test_crf.jl")
+    end
 end