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/crf.jl b/src/crf.jl new file mode 100644 index 0000000..11bfcea --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/model.jl b/src/model.jl index 204f534..b0efa7a 100644 --- a/src/model.jl +++ b/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 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 diff --git a/src/utils.jl b/src/utils.jl index 48201de..416204a 100644 --- a/src/utils.jl +++ b/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 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" 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..c43239f --- /dev/null +++ b/test/test_crf.jl @@ -0,0 +1,435 @@ +@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 + 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 for sp2 + t_indx = 2 + t_inv = collect(𝒯ᴵⁿᵛ)[t_indx] + + 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))/(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 + + 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 = [ + 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, 𝒯ᴵⁿᵛ) * EMI.set_period_annuity(inv_data, t)) for t in Tᶜᵘᵐ[t_inv]) + for i in 1:4, t_inv in 𝒯ᴵⁿᵛ] + @testset "Check explicit calculations with model's results" begin + @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 ∈ 𝒯ᴵⁿᵛ + ]) + 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 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] * 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([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, 𝒯ᴵⁿᵛ) * 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 + +@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 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 * (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, 𝒯ᴵⁿᵛ) * 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 + +@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 + +@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 + +@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 diff --git a/test/utils.jl b/test/utils.jl index 3fb201d..4a62d60 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -127,3 +127,8 @@ function variables(m, n, 𝒯) @variable(m, surplus[𝒯] ≥ 0) @variable(m, deficit[𝒯] ≥ 0) end + +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