From 467b87c97d32486e95afcf592cd12e07d381f2d8 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 14 Jul 2025 12:21:28 +0200 Subject: [PATCH 01/18] Include discount rates in AbstractInvData --- src/structures/investment_data.jl | 61 ++++++++++++++++++++++++++++++- 1 file changed, 59 insertions(+), 2 deletions(-) diff --git a/src/structures/investment_data.jl b/src/structures/investment_data.jl index 2fd815f..8d8f14a 100644 --- a/src/structures/investment_data.jl +++ b/src/structures/investment_data.jl @@ -30,6 +30,7 @@ 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, @@ -37,7 +38,25 @@ function NoStartInvData( 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 From f30dfc30d4ba5b45d39948194e643237c767073f Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 14 Jul 2025 12:22:18 +0200 Subject: [PATCH 02/18] Calculation of Capital Recovery Factor (CRF), add util to retrieve cumulative periods --- src/crf.jl | 15 +++++++++++++++ src/utils.jl | 13 +++++++++++++ 2 files changed, 28 insertions(+) create mode 100644 src/crf.jl diff --git a/src/crf.jl b/src/crf.jl new file mode 100644 index 0000000..4e8f812 --- /dev/null +++ b/src/crf.jl @@ -0,0 +1,15 @@ +""" + 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 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, 𝒯ᴵⁿᵛ) + return (disc_rate * (1 + disc_rate)^(remaining_horizon)) / ((1 + disc_rate)^(remaining_horizon) - 1) +end \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index 48201de..6b80b99 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -145,3 +145,16 @@ 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 From 697b456d8f33041cb41df6160e7262fc2388a49a Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 14 Jul 2025 12:23:28 +0200 Subject: [PATCH 03/18] Include annualisation in UnlimitedLife elements, include tests --- src/EnergyModelsInvestments.jl | 1 + src/model.jl | 7 ++++- test/runtests.jl | 4 +++ test/test_crf.jl | 56 ++++++++++++++++++++++++++++++++++ test/utils.jl | 4 +++ 5 files changed, 71 insertions(+), 1 deletion(-) create mode 100644 test/test_crf.jl diff --git a/src/EnergyModelsInvestments.jl b/src/EnergyModelsInvestments.jl index f511b7b..b8f404a 100644 --- a/src/EnergyModelsInvestments.jl +++ b/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") diff --git a/src/model.jl b/src/model.jl index 204f534..8a4be94 100644 --- a/src/model.jl +++ b/src/model.jl @@ -232,7 +232,12 @@ 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(𝒯ᴵⁿᵛ) + @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum(capex_val[t] * CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration 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 ∈ 𝒯ᴵⁿᵛ diff --git a/test/runtests.jl b/test/runtests.jl index 057738d..66e4664 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test_crf.jl b/test/test_crf.jl new file mode 100644 index 0000000..0118694 --- /dev/null +++ b/test/test_crf.jl @@ -0,0 +1,56 @@ + +@testset "UnlimitedLife - CRF" begin + # Creation and solving of the model + demand = StrategicProfile([10,30,30,20]) + inv_data = NoStartInvData( + FixedProfile(1000), + FixedProfile(30), + ContinuousInvestment(FixedProfile(0), FixedProfile(10)), + 0.07 + ) + m, para = simple_model(;demand=demand, inv_data=inv_data) + + # Extraction of required data + n = para[:node] + 𝒯 = para[:T] + 𝒯ᴵⁿᵛ = strat_periods(𝒯) + inv_data = para[:inv_data] + capex = StrategicProfile([1,1,1,0])*1e4 + + # Test the Annualised Capital Cost + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + + ## Test for sp2 + t_indx = 2 + capex_sp = StrategicProfile([0, 1, 0, 0])*1e4 + CRF = EMI.CRF(inv_data, collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ) + @testset "Calculation CRF" begin + @test isapprox(CRF, (0.07 * 1.07^30)/(1.07^30-1)) + end + + annualised_capex_sp = [ + sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ + ] + pv_annualised_capex_sp = present_value(annualised_capex_sp[t_indx:end], 0.07, 10) + @testset "Check annualised costs allocation and value" begin + @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 + @test all(annualised_capex_sp[t_indx:end] .== CRF * 1e4 * 10) # from sp2 onwards, annual costs allocated + + @test sum(pv_annualised_capex_sp) > capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]] # present value of sum of annual costs is higher than capex as it includes return costs + end + + vector_capex = [ + StrategicProfile([1, 0, 0, 0])*1e4, + StrategicProfile([0, 1, 0, 0])*1e4, + StrategicProfile([0, 0, 1, 0])*1e4, + StrategicProfile([0, 0, 0, 0])*1e4, + ] + annualised_capex = [ + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ + ] + @testset "Check with results" begin + @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) + end + +end \ No newline at end of file diff --git a/test/utils.jl b/test/utils.jl index 3fb201d..5cee2b2 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -127,3 +127,7 @@ function variables(m, n, 𝒯) @variable(m, surplus[𝒯] ≥ 0) @variable(m, deficit[𝒯] ≥ 0) end + +function present_value(vector, r, period_duration) + return [vector[i] * (1+r)^(-period_duration*(i - 1)) for i in 1:length(vector)] +end \ No newline at end of file From 83ffce4fd9b70b71fd7a7a7d78edecbe8a733438 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 14 Jul 2025 14:05:09 +0200 Subject: [PATCH 04/18] Include annualisation of StudyLIfe, included test --- src/model.jl | 26 ++++++++++----- test/runtests.jl | 22 ++++++------- test/test_crf.jl | 82 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 112 insertions(+), 18 deletions(-) diff --git a/src/model.jl b/src/model.jl index 8a4be94..28debbb 100644 --- a/src/model.jl +++ b/src/model.jl @@ -251,14 +251,26 @@ 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(𝒯ᴵⁿᵛ) + + @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], 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 = 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 ∈ 𝒯ᴵⁿᵛ diff --git a/test/runtests.jl b/test/runtests.jl index 66e4664..52b4713 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,16 +13,16 @@ const TEST_ATOL = 1e-6 include("utils.jl") -@testset "Investments" begin - @testset "Investments | Investment modes" begin - include("test_invest.jl") - end +# @testset "Investments" begin +# @testset "Investments | Investment modes" begin +# include("test_invest.jl") +# end - @testset "Investments | Lifetime" begin - include("test_lifetime.jl") - end +# @testset "Investments | Lifetime" begin +# include("test_lifetime.jl") +# end - @testset "Investments | Annualised" begin - include("test_crf.jl") - end -end +# @testset "Investments | Annualised" begin +# include("test_crf.jl") +# end +# end diff --git a/test/test_crf.jl b/test/test_crf.jl index 0118694..8b0c632 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -53,4 +53,86 @@ @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) end +end + + +@testset "StudyLife - CRF" begin + # Creation and solving of the model + demand = StrategicProfile([10,10,30,35]) + inv_data = NoStartInvData( + FixedProfile(1000), + FixedProfile(40), + ContinuousInvestment(FixedProfile(0), FixedProfile(15)), + StudyLife(FixedProfile(20)), + 0.07 + ) + m, para = simple_model(;demand=demand, inv_data=inv_data) + + # Extraction of required data + n = para[:node] + 𝒯 = para[:T] + 𝒯ᴵⁿᵛ = strat_periods(𝒯) + inv_data = para[:inv_data] + + capex = StrategicProfile([ + 10 , + 5 , + 15, + 5 , + ])*1e3 + capex_disc = StrategicProfile([ + EMI.set_capex_discounter( + remaining(t_inv, 𝒯ᴵⁿᵛ), + EMI.lifetime(inv_data, t_inv), EMI.get_discount_rate(inv_data) + ) for t_inv ∈ 𝒯ᴵⁿᵛ + ]) + + # Test the discount calculation + capex_explicit = StrategicProfile([ + 10 * (1 + 1/1.07^20), + 5 * (1 + (1/1.07^20 - 0.5 * 1/1.07^30)), + 15, + 5 * (1 - 0.5 * 1/1.07^10), + ])*1e3 + @testset "Discounted Capex calculations" + @test all(isapprox.([capex[t_inv] * capex_disc[t_inv] for t_inv in 𝒯ᴵⁿᵛ], [capex_explicit[t_inv] for t_inv in 𝒯ᴵⁿᵛ])) + end + + # Test the Annualised Capital Cost + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + + ## Test for sp2 + t_indx =2 + capex_sp = StrategicProfile([0, 5 * (1 + (1/1.07^20 - 0.5 * 1/1.07^30)), 0, 0])*1e3 + CRF = EMI.CRF(inv_data, collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ) + @testset "Calculation CRF" begin + @test isapprox(CRF, (0.07 * 1.07^30)/(1.07^30-1)) + end + + annualised_capex_sp = [ + sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ + ] + pv_annualised_capex_sp = present_value(annualised_capex_sp[t_indx:end], 0.07, 10) + @testset "Check annualised costs allocation and value" begin + @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 + @test all(isapprox.(annualised_capex_sp[t_indx:end], CRF * 1e3 * 5 * 10)) # from sp2 onwards, annual costs allocated + + @test sum(pv_annualised_capex_sp) > capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]] # present value of sum of annual costs is higher than capex as it includes return costs + end + + # Check all annualisation of capex + vector_capex = [ + StrategicProfile([10* (1 + 1/1.07^20), 0, 0, 0])*1e3, + StrategicProfile([0, 5* (1 + (1/1.07^20 - 0.5 * 1/1.07^30)), 0, 0])*1e3, + StrategicProfile([0, 0, 15, 0])*1e3, + StrategicProfile([0, 0, 0, 5* (1 - 0.5 * 1/1.07^10)])*1e3, + ] + annualised_capex = [ + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ + ] + @testset "Check with results" begin + @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) + end + end \ No newline at end of file From fb2034b1c970273e3d875faefb58c8b4e1aef4a5 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 14 Jul 2025 14:09:29 +0200 Subject: [PATCH 05/18] Fix small bug --- test/runtests.jl | 22 +++++++++++----------- test/test_crf.jl | 8 ++------ 2 files changed, 13 insertions(+), 17 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 52b4713..66e4664 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,16 +13,16 @@ const TEST_ATOL = 1e-6 include("utils.jl") -# @testset "Investments" begin -# @testset "Investments | Investment modes" begin -# include("test_invest.jl") -# end +@testset "Investments" begin + @testset "Investments | Investment modes" begin + include("test_invest.jl") + end -# @testset "Investments | Lifetime" begin -# include("test_lifetime.jl") -# end + @testset "Investments | Lifetime" begin + include("test_lifetime.jl") + end -# @testset "Investments | Annualised" begin -# include("test_crf.jl") -# end -# end + @testset "Investments | Annualised" begin + include("test_crf.jl") + end +end diff --git a/test/test_crf.jl b/test/test_crf.jl index 8b0c632..3b6cda5 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -94,7 +94,7 @@ end 15, 5 * (1 - 0.5 * 1/1.07^10), ])*1e3 - @testset "Discounted Capex calculations" + @testset "Discounted Capex calculations" begin @test all(isapprox.([capex[t_inv] * capex_disc[t_inv] for t_inv in 𝒯ᴵⁿᵛ], [capex_explicit[t_inv] for t_inv in 𝒯ᴵⁿᵛ])) end @@ -115,8 +115,6 @@ end pv_annualised_capex_sp = present_value(annualised_capex_sp[t_indx:end], 0.07, 10) @testset "Check annualised costs allocation and value" begin @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 - @test all(isapprox.(annualised_capex_sp[t_indx:end], CRF * 1e3 * 5 * 10)) # from sp2 onwards, annual costs allocated - @test sum(pv_annualised_capex_sp) > capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]] # present value of sum of annual costs is higher than capex as it includes return costs end @@ -129,10 +127,8 @@ end ] annualised_capex = [ sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) - for i in 1:4, t_inv in 𝒯ᴵⁿᵛ - ] + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) end - end \ No newline at end of file From 2c1bcad053bdacf5407574b53693681e21b7d1ba Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 14 Jul 2025 14:33:26 +0200 Subject: [PATCH 06/18] Include annualisation of PeriodLife --- src/model.jl | 27 ++++++++++++++----- test/test_crf.jl | 70 +++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 89 insertions(+), 8 deletions(-) diff --git a/src/model.jl b/src/model.jl index 28debbb..8455f98 100644 --- a/src/model.jl +++ b/src/model.jl @@ -286,14 +286,27 @@ 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(𝒯ᴵⁿᵛ) + + @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], 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 = 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) diff --git a/test/test_crf.jl b/test/test_crf.jl index 3b6cda5..e547d81 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -98,7 +98,7 @@ end @test all(isapprox.([capex[t_inv] * capex_disc[t_inv] for t_inv in 𝒯ᴵⁿᵛ], [capex_explicit[t_inv] for t_inv in 𝒯ᴵⁿᵛ])) end - # Test the Annualised Capital Cost + # Test the Annualised Capital Cost Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) ## Test for sp2 @@ -131,4 +131,72 @@ end @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) end +end + +@testset "PeriodLife - CRF" begin + # Creation and solving of the model + inv_data = NoStartInvData( + FixedProfile(1000), + FixedProfile(40), + ContinuousInvestment(FixedProfile(0), FixedProfile(15)), + PeriodLife(FixedProfile(20)), + 0.07 + ) + demand = StrategicProfile([5,10,15,15]) + m, para = simple_model(;inv_data,demand) + + # Extraction of required data + n = para[:node] + 𝒯 = para[:T] + 𝒯ᴵⁿᵛ = strat_periods(𝒯) + inv_data = para[:inv_data] + invest = StrategicProfile([5, 10, 15, 15])*1e3 + capex_explicit = invest * (1 - 0.5 * 1/1.07^10) + + capex = StrategicProfile([ + 5 , + 10 , + 15, + 15 , + ])*1e3 + capex_disc = StrategicProfile([ + EMI.set_capex_discounter( + duration_strat(t_inv), + EMI.lifetime(inv_data, t_inv), EMI.get_discount_rate(inv_data)) for t_inv ∈ 𝒯ᴵⁿᵛ]) + + @testset "Discounted Capex calculations" begin + @test all(isapprox.([capex[t_inv] * capex_disc[t_inv] for t_inv in 𝒯ᴵⁿᵛ], [capex_explicit[t_inv] for t_inv in 𝒯ᴵⁿᵛ])) + end + + # Test the Annualised Capital Cost + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + + ## Test for sp2 + t_indx =2 + capex_sp = StrategicProfile([0, 10 * (1 - 0.5 * 1/1.07^10), 0, 0])*1e3 + CRF = EMI.CRF(inv_data, collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ) + + annualised_capex_sp = [ + sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ + ] + pv_annualised_capex_sp = present_value(annualised_capex_sp[t_indx:end], 0.07, 10) + @testset "Check annualised costs allocation and value" begin + @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 + @test sum(pv_annualised_capex_sp) > capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]] # present value of sum of annual costs is higher than capex as it includes return costs + end + + # Check all annualisation of capex + vector_capex = [ + StrategicProfile([5 * (1 - 0.5 * 1/1.07^10), 0, 0, 0])*1e3, + StrategicProfile([0, 10 * (1 - 0.5 * 1/1.07^10), 0, 0])*1e3, + StrategicProfile([0, 0, 15 * (1 - 0.5 * 1/1.07^10), 0])*1e3, + StrategicProfile([0, 0, 0, 15* (1 - 0.5 * 1/1.07^10)])*1e3, + ] + annualised_capex = [ + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] + @testset "Check with results" begin + @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) + end + end \ No newline at end of file From 8142f5fe57abfa7f827470eaccfb5bb532904ad1 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 14 Jul 2025 19:50:56 +0200 Subject: [PATCH 07/18] Fix annuity calculation by removing * t.duration --- src/model.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/model.jl b/src/model.jl index 8455f98..8e55826 100644 --- a/src/model.jl +++ b/src/model.jl @@ -234,7 +234,7 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat capex_val = set_capex_value(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ) if has_discount_rate(inv_data) Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ) - @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum(capex_val[t] * CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv])) + @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum(capex_val[t] * CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv])) else @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == capex_val[t_inv]) end @@ -261,7 +261,7 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat ]) Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ) - @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum((capex_val[t] * capex_disc[t]) * CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv])) + @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum((capex_val[t] * capex_disc[t]) * CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv])) else capex_disc = StrategicProfile([ set_capex_discounter( @@ -296,7 +296,7 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat ]) Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ) - @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum((capex_val[t] * capex_disc[t]) * CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv])) + @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum((capex_val[t] * capex_disc[t]) * CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv])) else capex_disc = StrategicProfile([ set_capex_discounter( From 4b91ca73e3b05d91ad0369676402699da1236217 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 14 Jul 2025 19:54:58 +0200 Subject: [PATCH 08/18] Improve calculation present value --- test/test_crf.jl | 35 ++++++++++++++--------------------- test/utils.jl | 4 ++-- 2 files changed, 16 insertions(+), 23 deletions(-) diff --git a/test/test_crf.jl b/test/test_crf.jl index e547d81..483aa9f 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -29,14 +29,13 @@ end annualised_capex_sp = [ - sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ + sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ ] - pv_annualised_capex_sp = present_value(annualised_capex_sp[t_indx:end], 0.07, 10) + pv = present_value(annualised_capex_sp[t_indx], 0.07, TS.remaining(collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ)) @testset "Check annualised costs allocation and value" begin @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 - @test all(annualised_capex_sp[t_indx:end] .== CRF * 1e4 * 10) # from sp2 onwards, annual costs allocated - - @test sum(pv_annualised_capex_sp) > capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]] # present value of sum of annual costs is higher than capex as it includes return costs + @test all(annualised_capex_sp[t_indx:end] .== CRF * 1e4) # from sp2 onwards, annual costs allocated + @test isapprox(pv, capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]]) end vector_capex = [ @@ -46,9 +45,8 @@ StrategicProfile([0, 0, 0, 0])*1e4, ] annualised_capex = [ - sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) - for i in 1:4, t_inv in 𝒯ᴵⁿᵛ - ] + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) end @@ -105,17 +103,13 @@ end t_indx =2 capex_sp = StrategicProfile([0, 5 * (1 + (1/1.07^20 - 0.5 * 1/1.07^30)), 0, 0])*1e3 CRF = EMI.CRF(inv_data, collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ) - @testset "Calculation CRF" begin - @test isapprox(CRF, (0.07 * 1.07^30)/(1.07^30-1)) - end - annualised_capex_sp = [ - sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ + sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ ] - pv_annualised_capex_sp = present_value(annualised_capex_sp[t_indx:end], 0.07, 10) + pv = present_value(annualised_capex_sp[t_indx], 0.07, TS.remaining(collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ)) @testset "Check annualised costs allocation and value" begin @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 - @test sum(pv_annualised_capex_sp) > capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]] # present value of sum of annual costs is higher than capex as it includes return costs + @test isapprox(pv, capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]]) end # Check all annualisation of capex @@ -126,7 +120,7 @@ end StrategicProfile([0, 0, 0, 5* (1 - 0.5 * 1/1.07^10)])*1e3, ] annualised_capex = [ - sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) @@ -175,14 +169,13 @@ end t_indx =2 capex_sp = StrategicProfile([0, 10 * (1 - 0.5 * 1/1.07^10), 0, 0])*1e3 CRF = EMI.CRF(inv_data, collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ) - annualised_capex_sp = [ - sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ + sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ ] - pv_annualised_capex_sp = present_value(annualised_capex_sp[t_indx:end], 0.07, 10) + pv = present_value(annualised_capex_sp[t_indx], 0.07, TS.remaining(collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ)) @testset "Check annualised costs allocation and value" begin @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 - @test sum(pv_annualised_capex_sp) > capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]] # present value of sum of annual costs is higher than capex as it includes return costs + @test isapprox(pv, capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]]) end # Check all annualisation of capex @@ -193,7 +186,7 @@ end StrategicProfile([0, 0, 0, 15* (1 - 0.5 * 1/1.07^10)])*1e3, ] annualised_capex = [ - sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * t.duration for t in Tᶜᵘᵐ[t_inv]) + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) diff --git a/test/utils.jl b/test/utils.jl index 5cee2b2..11c91f2 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -128,6 +128,6 @@ function variables(m, n, 𝒯) @variable(m, deficit[𝒯] ≥ 0) end -function present_value(vector, r, period_duration) - return [vector[i] * (1+r)^(-period_duration*(i - 1)) for i in 1:length(vector)] +function present_value(annualised_cost, r, duration) + return annualised_cost * (1 - (1 + r)^(-duration))/r end \ No newline at end of file From 02d0b6e4e58e2199af792604d87a6bd8c10c5cc2 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 10:16:01 +0200 Subject: [PATCH 09/18] Fix Annuity calculation for periods, include these in the constraints for UnlimitedLIfe, modify and check tests --- src/crf.jl | 15 +++++++++++++-- src/model.jl | 4 +++- test/test_crf.jl | 35 ++++++++++++++++++++--------------- test/utils.jl | 5 +++-- 4 files changed, 39 insertions(+), 20 deletions(-) diff --git a/src/crf.jl b/src/crf.jl index 4e8f812..11bfcea 100644 --- a/src/crf.jl +++ b/src/crf.jl @@ -2,7 +2,10 @@ 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 `𝒯ᴵⁿᵛ`. +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 @@ -11,5 +14,13 @@ over its economic lifetime. function CRF(inv_data::AbstractInvData, t_inv, 𝒯ᴵⁿᵛ) disc_rate = get_discount_rate(inv_data) remaining_horizon = remaining(t_inv, 𝒯ᴵⁿᵛ) - return (disc_rate * (1 + disc_rate)^(remaining_horizon)) / ((1 + disc_rate)^(remaining_horizon) - 1) + 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 \ No newline at end of file diff --git a/src/model.jl b/src/model.jl index 8e55826..4e1e157 100644 --- a/src/model.jl +++ b/src/model.jl @@ -234,7 +234,9 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat capex_val = set_capex_value(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ) if has_discount_rate(inv_data) Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ) - @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum(capex_val[t] * CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv])) + 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 diff --git a/test/test_crf.jl b/test/test_crf.jl index 483aa9f..9a462b7 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -15,27 +15,32 @@ 𝒯 = para[:T] 𝒯ᴵⁿᵛ = strat_periods(𝒯) inv_data = para[:inv_data] - capex = StrategicProfile([1,1,1,0])*1e4 - - # Test the Annualised Capital Cost - Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + capex = StrategicProfile([1,1,1,0])*1e4 ## Test for sp2 t_indx = 2 + t_inv = collect(𝒯ᴵⁿᵛ)[t_indx] + capex_sp = StrategicProfile([0, 1, 0, 0])*1e4 - CRF = EMI.CRF(inv_data, collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ) + CRF = EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) @testset "Calculation CRF" begin @test isapprox(CRF, (0.07 * 1.07^30)/(1.07^30-1)) end + annuity_capex = capex_sp[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities + period_annuity_capex = annuity_capex * EMI.set_period_annuity(inv_data, t_inv) # compute the value of single payments at the beginning of each period + pv = present_value(period_annuity_capex, 0.07, t_inv.duration, 3) + @testset "Check present value of annuity and period annuity" begin + @test isapprox(pv_annuity, capex_sp[t_inv]) + @test isapprox(pv, capex_sp[t_inv]) + end - annualised_capex_sp = [ - sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ - ] - pv = present_value(annualised_capex_sp[t_indx], 0.07, TS.remaining(collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ)) - @testset "Check annualised costs allocation and value" begin - @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 - @test all(annualised_capex_sp[t_indx:end] .== CRF * 1e4) # from sp2 onwards, annual costs allocated - @test isapprox(pv, capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]]) + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + annuity_capex = Dict(t_inv => capex_sp[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) for t_inv in 𝒯ᴵⁿᵛ) + period_annuity_capex = Dict(t_inv => annuity_capex[t_inv] * EMI.set_period_annuity(inv_data, t_inv) for t_inv in 𝒯ᴵⁿᵛ) + cap_capex = Dict(t_inv => sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ) + @testset "Check assignment of period_annuities" begin + @test cap_capex[first(collect(𝒯ᴵⁿᵛ))] == 0 # sp1 has 0 cost from investments in sp2 + @test all([cap_capex[t] == period_annuity_capex[t_inv] for t in collect(𝒯ᴵⁿᵛ)[t_indx:end]]) # from sp2 onwards, annual costs allocated end vector_capex = [ @@ -45,9 +50,9 @@ StrategicProfile([0, 0, 0, 0])*1e4, ] annualised_capex = [ - sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) + sum((vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * EMI.set_period_annuity(inv_data, t)) for t in Tᶜᵘᵐ[t_inv]) for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] - @testset "Check with results" begin + @testset "Check explicit calculations with model's results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) end diff --git a/test/utils.jl b/test/utils.jl index 11c91f2..4a62d60 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -128,6 +128,7 @@ function variables(m, n, 𝒯) @variable(m, deficit[𝒯] ≥ 0) end -function present_value(annualised_cost, r, duration) - return annualised_cost * (1 - (1 + r)^(-duration))/r +function present_value(annuity_period, r, period_duration, n_periods) + pv = sum([annuity_period/(1+r)^(period_duration*i) for i in 0:(n_periods-1)]) + return pv end \ No newline at end of file From 031d8b77bfb942ab60d3562deffc99b8863c435a Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 10:43:38 +0200 Subject: [PATCH 10/18] Fix Annuity calculation for periods, include these into StudyLife and modify test --- src/model.jl | 17 +++++++++++++---- test/test_crf.jl | 47 ++++++++++++++++++++++++++++------------------- 2 files changed, 41 insertions(+), 23 deletions(-) diff --git a/src/model.jl b/src/model.jl index 4e1e157..64c4fb1 100644 --- a/src/model.jl +++ b/src/model.jl @@ -262,8 +262,10 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat ) for t_inv ∈ 𝒯ᴵⁿᵛ ]) Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ) - - @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum((capex_val[t] * capex_disc[t]) * CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv])) + + 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( @@ -323,6 +325,8 @@ 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 ∈ 𝒯ᴵⁿᵛ) + Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ) + for t_inv ∈ 𝒯ᴵⁿᵛ # Extract the values lifetime_val = lifetime(inv_data, t_inv) @@ -335,8 +339,13 @@ 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 = StrategicProfile([set_capex_discounter(duration_strat(t_inv), lifetime_val, get_discount_rate(inv_data)) for t_inv in 𝒯ᴵⁿᵛ]) + @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 diff --git a/test/test_crf.jl b/test/test_crf.jl index 9a462b7..5868300 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -24,9 +24,11 @@ capex_sp = StrategicProfile([0, 1, 0, 0])*1e4 CRF = EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) @testset "Calculation CRF" begin - @test isapprox(CRF, (0.07 * 1.07^30)/(1.07^30-1)) + @test isapprox(CRF, (0.07 * 1.07^(30-1))/(1.07^30-1)) end - annuity_capex = capex_sp[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities + + annuity_capex = capex_sp[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH + pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities period_annuity_capex = annuity_capex * EMI.set_period_annuity(inv_data, t_inv) # compute the value of single payments at the beginning of each period pv = present_value(period_annuity_capex, 0.07, t_inv.duration, 3) @testset "Check present value of annuity and period annuity" begin @@ -58,7 +60,6 @@ end - @testset "StudyLife - CRF" begin # Creation and solving of the model demand = StrategicProfile([10,10,30,35]) @@ -89,32 +90,40 @@ end EMI.lifetime(inv_data, t_inv), EMI.get_discount_rate(inv_data) ) for t_inv ∈ 𝒯ᴵⁿᵛ ]) - - # Test the discount calculation capex_explicit = StrategicProfile([ 10 * (1 + 1/1.07^20), 5 * (1 + (1/1.07^20 - 0.5 * 1/1.07^30)), 15, 5 * (1 - 0.5 * 1/1.07^10), ])*1e3 + @testset "Discounted Capex calculations" begin @test all(isapprox.([capex[t_inv] * capex_disc[t_inv] for t_inv in 𝒯ᴵⁿᵛ], [capex_explicit[t_inv] for t_inv in 𝒯ᴵⁿᵛ])) end - # Test the Annualised Capital Cost - Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + # Test for sp2 + t_indx = 2 + t_inv = collect(𝒯ᴵⁿᵛ)[t_indx] - ## Test for sp2 - t_indx =2 - capex_sp = StrategicProfile([0, 5 * (1 + (1/1.07^20 - 0.5 * 1/1.07^30)), 0, 0])*1e3 - CRF = EMI.CRF(inv_data, collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ) - annualised_capex_sp = [ - sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ - ] - pv = present_value(annualised_capex_sp[t_indx], 0.07, TS.remaining(collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ)) - @testset "Check annualised costs allocation and value" begin - @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 - @test isapprox(pv, capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]]) + capex_sp = StrategicProfile([0, 5, 0, 0])*1e3 + CRF = EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) + + annuity_capex = capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH + pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities + period_annuity_capex = annuity_capex * EMI.set_period_annuity(inv_data, t_inv) # compute the value of single payments at the beginning of each period + pv = present_value(period_annuity_capex, 0.07, t_inv.duration, 3) + @testset "Check present value of annuity and period annuity" begin + @test isapprox(pv_annuity, (capex_sp[t_inv] * capex_disc[t_inv])) + @test isapprox(pv, capex_sp[t_inv]* capex_disc[t_inv]) + end + + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + annuity_capex = Dict(t_inv => capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) for t_inv in 𝒯ᴵⁿᵛ) + period_annuity_capex = Dict(t_inv => annuity_capex[t_inv] * EMI.set_period_annuity(inv_data, t_inv) for t_inv in 𝒯ᴵⁿᵛ) + cap_capex = Dict(t_inv => sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ) + @testset "Check assignment of period_annuities" begin + @test cap_capex[first(collect(𝒯ᴵⁿᵛ))] == 0 # sp1 has 0 cost from investments in sp2 + @test all([cap_capex[t] == period_annuity_capex[t_inv] for t in collect(𝒯ᴵⁿᵛ)[t_indx:end]]) # from sp2 onwards, annual costs allocated end # Check all annualisation of capex @@ -125,7 +134,7 @@ end StrategicProfile([0, 0, 0, 5* (1 - 0.5 * 1/1.07^10)])*1e3, ] annualised_capex = [ - sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * EMI.set_period_annuity(inv_data, t) for t in Tᶜᵘᵐ[t_inv]) for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) From 6ea596a8837328ef8d50afc34a6e6daa518e89e0 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 10:52:24 +0200 Subject: [PATCH 11/18] Include fixed annuity into PeriodLife and modify test --- src/model.jl | 6 ++++-- test/test_crf.jl | 43 ++++++++++++++++++++++++++++--------------- 2 files changed, 32 insertions(+), 17 deletions(-) diff --git a/src/model.jl b/src/model.jl index 64c4fb1..18b1c95 100644 --- a/src/model.jl +++ b/src/model.jl @@ -299,8 +299,10 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat ) for t_inv ∈ 𝒯ᴵⁿᵛ ]) Tᶜᵘᵐ = get_cumulative_periods(𝒯ᴵⁿᵛ) - - @constraint(m, [t_inv ∈ 𝒯ᴵⁿᵛ], var_capex[t_inv] == sum((capex_val[t] * capex_disc[t]) * CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv])) + + 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( diff --git a/test/test_crf.jl b/test/test_crf.jl index 5868300..73c3705 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -159,8 +159,8 @@ end 𝒯ᴵⁿᵛ = strat_periods(𝒯) inv_data = para[:inv_data] invest = StrategicProfile([5, 10, 15, 15])*1e3 + capex_explicit = invest * (1 - 0.5 * 1/1.07^10) - capex = StrategicProfile([ 5 , 10 , @@ -176,20 +176,32 @@ end @test all(isapprox.([capex[t_inv] * capex_disc[t_inv] for t_inv in 𝒯ᴵⁿᵛ], [capex_explicit[t_inv] for t_inv in 𝒯ᴵⁿᵛ])) end - # Test the Annualised Capital Cost - Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + # Test for sp2 + t_indx = 2 + t_inv = collect(𝒯ᴵⁿᵛ)[t_indx] - ## Test for sp2 - t_indx =2 - capex_sp = StrategicProfile([0, 10 * (1 - 0.5 * 1/1.07^10), 0, 0])*1e3 - CRF = EMI.CRF(inv_data, collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ) - annualised_capex_sp = [ - sum(capex_sp[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ - ] - pv = present_value(annualised_capex_sp[t_indx], 0.07, TS.remaining(collect(𝒯ᴵⁿᵛ)[t_indx], 𝒯ᴵⁿᵛ)) - @testset "Check annualised costs allocation and value" begin - @test first(annualised_capex_sp) == 0 # sp1 has 0 cost from investments in sp2 - @test isapprox(pv, capex_sp[collect(𝒯ᴵⁿᵛ)[t_indx]]) + capex_sp = StrategicProfile([0, 10, 0, 0])*1e3 + CRF = EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) + + annuity_capex = capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH + @testset "Check discounted capex with explicit capex" begin + @test isapprox(capex_sp[t_inv] * capex_disc[t_inv], capex_explicit[t_inv]) + end + pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities + period_annuity_capex = annuity_capex * EMI.set_period_annuity(inv_data, t_inv) # compute the value of single payments at the beginning of each period + pv = present_value(period_annuity_capex, 0.07, t_inv.duration, 3) + @testset "Check present value of annuity and period annuity" begin + @test isapprox(pv_annuity, (capex_sp[t_inv] * capex_disc[t_inv])) + @test isapprox(pv, capex_sp[t_inv]* capex_disc[t_inv]) + end + + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + annuity_capex = Dict(t_inv => capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) for t_inv in 𝒯ᴵⁿᵛ) + period_annuity_capex = Dict(t_inv => annuity_capex[t_inv] * EMI.set_period_annuity(inv_data, t_inv) for t_inv in 𝒯ᴵⁿᵛ) + cap_capex = Dict(t_inv => sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ) + @testset "Check assignment of period_annuities" begin + @test cap_capex[first(collect(𝒯ᴵⁿᵛ))] == 0 # sp1 has 0 cost from investments in sp2 + @test all([cap_capex[t] == period_annuity_capex[t_inv] for t in collect(𝒯ᴵⁿᵛ)[t_indx:end]]) # from sp2 onwards, annual costs allocated end # Check all annualisation of capex @@ -200,10 +212,11 @@ end StrategicProfile([0, 0, 0, 15* (1 - 0.5 * 1/1.07^10)])*1e3, ] annualised_capex = [ - sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) for t in Tᶜᵘᵐ[t_inv]) + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * EMI.set_period_annuity(inv_data, t) for t in Tᶜᵘᵐ[t_inv]) for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) end +end end \ No newline at end of file From 10a4ddd17b605b62b3882da5c612e56a355588d7 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 11:27:56 +0200 Subject: [PATCH 12/18] Include annuity calculation for RollingLIfe - Shorter Lifetime than period, add test --- src/model.jl | 15 +++++++--- test/test_crf.jl | 71 ++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 79 insertions(+), 7 deletions(-) diff --git a/src/model.jl b/src/model.jl index 18b1c95..9696052 100644 --- a/src/model.jl +++ b/src/model.jl @@ -326,10 +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) @@ -342,8 +343,14 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat # to account for the required reinvestments if lifetime_val < duration_strat(t_inv) if has_discount_rate(inv_data) - capex_disc = StrategicProfile([set_capex_discounter(duration_strat(t_inv), lifetime_val, get_discount_rate(inv_data)) for t_inv in 𝒯ᴵⁿᵛ]) - @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])) + 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) diff --git a/test/test_crf.jl b/test/test_crf.jl index 73c3705..5b77cd7 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -184,9 +184,6 @@ end CRF = EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) annuity_capex = capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH - @testset "Check discounted capex with explicit capex" begin - @test isapprox(capex_sp[t_inv] * capex_disc[t_inv], capex_explicit[t_inv]) - end pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities period_annuity_capex = annuity_capex * EMI.set_period_annuity(inv_data, t_inv) # compute the value of single payments at the beginning of each period pv = present_value(period_annuity_capex, 0.07, t_inv.duration, 3) @@ -219,4 +216,72 @@ end end end + +@testset "RollingLife - Shorter Lifetime - CRD" begin + # Creation and solving of the model + inv_data = NoStartInvData( + FixedProfile(1000), + FixedProfile(40), + ContinuousInvestment(FixedProfile(0), FixedProfile(15)), + RollingLife(FixedProfile(5)), + 0.07 + ) + demand = StrategicProfile([5,10,15,15]) + m, para = simple_model(;inv_data,demand) + + # Extraction of required data + n = para[:node] + 𝒯 = para[:T] + 𝒯ᴵⁿᵛ = strat_periods(𝒯) + + capex = StrategicProfile([5, 10, 15, 15]) * 1e3 + capex_explicit = capex * (1 + 1/(1.07)^5) + capex_disc = StrategicProfile([ + EMI.set_capex_discounter( + duration_strat(t_inv), + EMI.lifetime(inv_data, t_inv), EMI.get_discount_rate(inv_data)) for t_inv ∈ 𝒯ᴵⁿᵛ]) + @testset "Discounted Capex calculations" begin + @test all(isapprox.([capex[t_inv] * capex_disc[t_inv] for t_inv in 𝒯ᴵⁿᵛ], [capex_explicit[t_inv] for t_inv in 𝒯ᴵⁿᵛ])) + end + + # Test for sp2 + t_indx = 2 + t_inv = collect(𝒯ᴵⁿᵛ)[t_indx] + + capex_sp = StrategicProfile([0, 10, 0, 0])*1e3 + CRF = EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) + + annuity_capex = capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH + pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities + period_annuity_capex = annuity_capex * EMI.set_period_annuity(inv_data, t_inv) # compute the value of single payments at the beginning of each period + pv = present_value(period_annuity_capex, 0.07, t_inv.duration, 3) + @testset "Check present value of annuity and period annuity" begin + @test isapprox(pv_annuity, (capex_sp[t_inv] * capex_disc[t_inv])) + @test isapprox(pv, capex_sp[t_inv]* capex_disc[t_inv]) + end + + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + annuity_capex = Dict(t_inv => capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) for t_inv in 𝒯ᴵⁿᵛ) + period_annuity_capex = Dict(t_inv => annuity_capex[t_inv] * EMI.set_period_annuity(inv_data, t_inv) for t_inv in 𝒯ᴵⁿᵛ) + cap_capex = Dict(t_inv => sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ) + @testset "Check assignment of period_annuities" begin + @test cap_capex[first(collect(𝒯ᴵⁿᵛ))] == 0 # sp1 has 0 cost from investments in sp2 + @test all([cap_capex[t] == period_annuity_capex[t_inv] for t in collect(𝒯ᴵⁿᵛ)[t_indx:end]]) # from sp2 onwards, annual costs allocated + end + + # Check all annualisation of capex + vector_capex = [ + StrategicProfile([5, 0, 0, 0]) * 1e3 * (1 + 1/(1.07)^5), + StrategicProfile([0, 10, 0, 0]) * 1e3 * (1 + 1/(1.07)^5), + StrategicProfile([0, 0, 15, 0]) * 1e3 * (1 + 1/(1.07)^5), + StrategicProfile([0, 0, 0, 15]) * 1e3 * (1 + 1/(1.07)^5), + ] + annualised_capex = [ + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * EMI.set_period_annuity(inv_data, t) for t in Tᶜᵘᵐ[t_inv]) + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] + @testset "Check with results" begin + @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) + end + + end \ No newline at end of file From 830a3674511b306a6d9dc3c17cbc08f027a93a2e Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 12:07:46 +0200 Subject: [PATCH 13/18] Include annuity calculation for RollingLIfe - Equal Lifetime than period, add test --- src/model.jl | 10 ++++++++- test/test_crf.jl | 58 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 67 insertions(+), 1 deletion(-) diff --git a/src/model.jl b/src/model.jl index 9696052..7ed5190 100644 --- a/src/model.jl +++ b/src/model.jl @@ -360,7 +360,15 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat # 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 diff --git a/test/test_crf.jl b/test/test_crf.jl index 5b77cd7..b29d088 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -282,6 +282,64 @@ end @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) end +end + +@testset "RollingLife - Equal Lifetime - CRD" begin + # Creation and solving of the model + inv_data = NoStartInvData( + FixedProfile(1000), + FixedProfile(40), + ContinuousInvestment(FixedProfile(0), FixedProfile(15)), + RollingLife(FixedProfile(10)), + 0.07 + ) + demand = StrategicProfile([5,10,15,15]) + m, para = simple_model(;inv_data,demand) + + # Extraction of required data + n = para[:node] + 𝒯 = para[:T] + 𝒯ᴵⁿᵛ = strat_periods(𝒯) + + capex = StrategicProfile([5, 10, 15, 15]) * 1e3 + capex_explicit = capex + + # Test for sp2 + t_indx = 2 + t_inv = collect(𝒯ᴵⁿᵛ)[t_indx] + + capex_sp = StrategicProfile([0, 10, 0, 0])*1e3 + CRF = EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) + annuity_capex = capex_sp[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH + pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities + period_annuity_capex = annuity_capex * EMI.set_period_annuity(inv_data, t_inv) # compute the value of single payments at the beginning of each period + pv = present_value(period_annuity_capex, 0.07, t_inv.duration, 3) + @testset "Check present value of annuity and period annuity" begin + @test isapprox(pv_annuity, (capex_sp[t_inv])) + @test isapprox(pv, capex_sp[t_inv]) + end + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + annuity_capex = Dict(t_inv => capex_sp[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) for t_inv in 𝒯ᴵⁿᵛ) + period_annuity_capex = Dict(t_inv => annuity_capex[t_inv] * EMI.set_period_annuity(inv_data, t_inv) for t_inv in 𝒯ᴵⁿᵛ) + cap_capex = Dict(t_inv => sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ) + @testset "Check assignment of period_annuities" begin + @test cap_capex[first(collect(𝒯ᴵⁿᵛ))] == 0 # sp1 has 0 cost from investments in sp2 + @test all([cap_capex[t] == period_annuity_capex[t_inv] for t in collect(𝒯ᴵⁿᵛ)[t_indx:end]]) # from sp2 onwards, annual costs allocated + end + + # Check all annualisation of capex + vector_capex = [ + StrategicProfile([5, 0, 0, 0]) * 1e3, + StrategicProfile([0, 10, 0, 0]) * 1e3, + StrategicProfile([0, 0, 15, 0]) * 1e3, + StrategicProfile([0, 0, 0, 15]) * 1e3, + ] + annualised_capex = [ + sum(vector_capex[i][t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * EMI.set_period_annuity(inv_data, t) for t in Tᶜᵘᵐ[t_inv]) + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] + @testset "Check with results" begin + @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) + end end \ No newline at end of file From 0c9d845c174dffde8090e69fafc2e3ab3073ab5c Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 12:44:28 +0200 Subject: [PATCH 14/18] create get_capex_disc to calculate the capex_disc and rem_dict used in longer lifetimes of RollingLife --- src/utils.jl | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/utils.jl b/src/utils.jl index 6b80b99..416204a 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -158,3 +158,36 @@ function get_cumulative_periods(𝒯::TS.AbstractStratPers) 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 From a9c7fb54d4486d6919c97982fda8d71d591ff616 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 12:45:00 +0200 Subject: [PATCH 15/18] Include annuity calculation for RollingLife - Longer Lifetime, add test --- src/model.jl | 42 ++++++++++-------------------- test/test_crf.jl | 66 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+), 29 deletions(-) diff --git a/src/model.jl b/src/model.jl index 7ed5190..073cea1 100644 --- a/src/model.jl +++ b/src/model.jl @@ -373,36 +373,20 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat # 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 - 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) + 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, 𝒯ᴵⁿᵛ) + @show capex_disc + 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 + @show period_annuity_capex - # 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) + @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 end end for (t_inv_rem, t_inv_vec) ∈ rem_dict diff --git a/test/test_crf.jl b/test/test_crf.jl index b29d088..31381fa 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -342,4 +342,70 @@ end @testset "Check with results" begin @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) end +end + +@testset "RollingLife - Longer Lifetime - CRD" begin + # Creation and solving of the model + inv_data = NoStartInvData( + FixedProfile(1000), + FixedProfile(40), + ContinuousInvestment(FixedProfile(0), FixedProfile(15)), + RollingLife(FixedProfile(20)), + 0.07 + ) + demand = StrategicProfile([5,10,15,15]) + m, para = simple_model(;inv_data,demand) + + # Extraction of required data + n = para[:node] + 𝒯 = para[:T] + 𝒯ᴵⁿᵛ = strat_periods(𝒯) + + capex = StrategicProfile([5, 5, 10, 5]) * 1e3 + capex_explicit = StrategicProfile([5, 5, 10, 5 * (1 - 0.5*1/(1.07)^10)]) * 1e3 + rem_dict = Dict(t_inv => eltype(𝒯ᴵⁿᵛ)[] for t_inv ∈ 𝒯ᴵⁿᵛ) + capex_disc = StrategicProfile([EMI.get_capex_disc(20, 0.07, rem_dict, t_inv, t_inv, 𝒯ᴵⁿᵛ)[1] for t_inv ∈ 𝒯ᴵⁿᵛ]) + + @testset "Discounted Capex calculations" begin + @test all(isapprox.([capex[t_inv] * capex_disc[t_inv] for t_inv in 𝒯ᴵⁿᵛ], [capex_explicit[t_inv] for t_inv in 𝒯ᴵⁿᵛ])) + end + + # Test for sp2 + t_indx = 2 + t_inv = collect(𝒯ᴵⁿᵛ)[t_indx] + + capex_sp = StrategicProfile([0, 5, 0, 0])*1e3 + CRF = EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) + + annuity_capex = capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) # calculate the annuity extended for the entire TH + pv_annuity = present_value(annuity_capex, 0.07, 1, 30) # compute the present value of the annuities + period_annuity_capex = annuity_capex * EMI.set_period_annuity(inv_data, t_inv) # compute the value of single payments at the beginning of each period + pv = present_value(period_annuity_capex, 0.07, t_inv.duration, 3) + @testset "Check present value of annuity and period annuity" begin + @test isapprox(pv_annuity, (capex_sp[t_inv])) + @test isapprox(pv, capex_sp[t_inv]) + end + + Tᶜᵘᵐ = EMI.get_cumulative_periods(𝒯ᴵⁿᵛ) + annuity_capex = Dict(t_inv => capex_sp[t_inv] * capex_disc[t_inv] * EMI.CRF(inv_data, t_inv, 𝒯ᴵⁿᵛ) for t_inv in 𝒯ᴵⁿᵛ) + period_annuity_capex = Dict(t_inv => annuity_capex[t_inv] * EMI.set_period_annuity(inv_data, t_inv) for t_inv in 𝒯ᴵⁿᵛ) + cap_capex = Dict(t_inv => sum(period_annuity_capex[t] for t in Tᶜᵘᵐ[t_inv]) for t_inv in 𝒯ᴵⁿᵛ) + @testset "Check assignment of period_annuities" begin + @test cap_capex[first(collect(𝒯ᴵⁿᵛ))] == 0 # sp1 has 0 cost from investments in sp2 + @test all([cap_capex[t] == period_annuity_capex[t_inv] for t in collect(𝒯ᴵⁿᵛ)[t_indx:end]]) # from sp2 onwards, annual costs allocated + end + + # Check all annualisation of capex + vector_capex = [ + StrategicProfile([5, 0, 0, 0]) * 1e3, + StrategicProfile([0, 5, 0, 0]) * 1e3, + StrategicProfile([0, 0, 10, 0]) * 1e3, + StrategicProfile([0, 0, 0, 5]) * 1e3, + ] + annualised_capex = [ + sum(vector_capex[i][t] * capex_disc[t] * EMI.CRF(inv_data, t, 𝒯ᴵⁿᵛ) * EMI.set_period_annuity(inv_data, t) for t in Tᶜᵘᵐ[t_inv]) + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] + @testset "Check with results" begin + @test all(isapprox.(sum(annualised_capex, dims=1), value.(m[:cap_capex]))) + end end \ No newline at end of file From b443e2b0746ec88d016200202fe3c95a8091e317 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 12:52:10 +0200 Subject: [PATCH 16/18] Add proof of concept --- test/test_crf.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/test/test_crf.jl b/test/test_crf.jl index 31381fa..c43239f 100644 --- a/test/test_crf.jl +++ b/test/test_crf.jl @@ -1,3 +1,27 @@ +@testset "Proof of concept - CRF" begin + demand = StrategicProfile([10,30,30,20]) + + inv_data1 = NoStartInvData( + FixedProfile(1000), + FixedProfile(30), + ContinuousInvestment(FixedProfile(0), FixedProfile(10)), + 0.07 # riskier technology + ) + m1, para1 = simple_model(;demand=demand, inv_data=inv_data1) + + inv_data2 = NoStartInvData( + FixedProfile(1000), + FixedProfile(30), + ContinuousInvestment(FixedProfile(0), FixedProfile(10)), + 0.02 + ) + m2, para2 = simple_model(;demand=demand, inv_data=inv_data2) + + @testset "Comparing cap_capex" begin + @test all(value.(m1[:cap_capex]).data .> value.(m2[:cap_capex]).data) + end + +end @testset "UnlimitedLife - CRF" begin # Creation and solving of the model From be94bf8c38a1510e6002be47ee81ce882e2905a4 Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Tue, 15 Jul 2025 12:54:17 +0200 Subject: [PATCH 17/18] remove @show --- src/model.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/model.jl b/src/model.jl index 073cea1..b0efa7a 100644 --- a/src/model.jl +++ b/src/model.jl @@ -376,11 +376,9 @@ function set_capacity_cost(m, element, inv_data, prefix, 𝒯ᴵⁿᵛ, disc_rat 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, 𝒯ᴵⁿᵛ) - @show capex_disc 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 - @show period_annuity_capex @constraint(m, var_capex[t_inv] == sum(period_annuity_capex_dict[t] for t in Tᶜᵘᵐ[t_inv])) else From e6199c421ac91be9f3c8bea962e762d89fb8f58d Mon Sep 17 00:00:00 2001 From: Raquel Alonso Date: Mon, 11 Aug 2025 15:08:26 +0200 Subject: [PATCH 18/18] Project.toml test --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index 81510f0..eb6c45e 100644 --- a/test/Project.toml +++ b/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"