Results differences between code versions

[HauHe]

For the European power sector model we us an OSeMOSYS code that contains custom equations for the ReserveMargin. Therefore, updates to the OSeMOSYS code cannot be directly used.

A few days ago I have updated the code using the osemosys_fast code from GitHub. However, when running OSeMBE with the updated code the results are notably different from the old code. The scenario I ran is the ECEMF WP1 NetZero scenario which is set-up to reach negative emissions in 2050, i.e., emissions are 0 in a year before 2050. This is done by using an EmissionPenalty. For my use case I have tested the impact that the availability of CCS has on the results. Before 2040 the model can't invest in CCS, only from 2040 onwards it may do so.

With the old code the model starts investing into biomass CCS directly in 2040. In 2040 a bit less than 10 time more than in the following years.
With the updated code however, the investments into BM CCS 2040 are 10 times higher than with the old code.

I guess this behaviour is not desired. I wonder if this might be caused by something in the updated formulation of discounting and or the issue of considering emissions penalties for negative emissions #99.

I attach the two osemosys versions. osemosys_osembe_ecemf.txt is the old code version. osemosys_osembe.txt the new code version.

osemosys_osembe_ecemf.txt
osemosys_osembe.txt

[willu47]

There are indeed a number of differences across the two files you shared and the current version in this repo. Taking just the objective function:

osembe_ecemf

minimize cost: sum{r in REGION, t in TECHNOLOGY, y in YEAR}
    (((((sum{yy in YEAR: y-yy < OperationalLife[r,t] && y-yy>=0}
        NewCapacity[r,t,yy])+ ResidualCapacity[r,t,y])*FixedCost[r,t,y] +
        sum{m in MODEperTECHNOLOGY[t], l in TIMESLICE}
            RateOfActivity[r,l,t,m,y]*YearSplit[l,y]*VariableCost[r,t,m,y]) / DiscountFactorMid[r,t,y] +
            CapitalCost[r,t,y] * NewCapacity[r,t,y] / DiscountFactor[r,t,y] +
            DiscountedTechnologyEmissionsPenalty[r,t,y] - DiscountedSalvageValue[r,t,y]) +
            sum{s in STORAGE} (CapitalCostStorage[r,s,y] * NewStorageCapacity[r,s,y] / DiscountFactor[r,t,y] -
            CapitalCostStorage[r,s,y] * NewStorageCapacity[r,s,y] / DiscountFactor[r,t,y]));

osembe

minimize cost: sum{r in REGION, t in TECHNOLOGY, y in YEAR}
    ((((sum{yy in YEAR: y-yy < OperationalLife[r,t] && y-yy>=0}
        NewCapacity[r,t,yy])+ ResidualCapacity[r,t,y]) * FixedCost[r,t,y] +
        sum{m in MODEperTECHNOLOGY[t], l in TIMESLICE}
            RateOfActivity[r,l,t,m,y] * YearSplit[l,y]*VariableCost[r,t,m,y]) / DiscountFactorMid[r,y] +
            CapitalCost[r,t,y] * NewCapacity[r,t,y] * CapitalRecoveryFactor[r,t] * PvAnnuity[r,t] / DiscountFactor[r,y] +
            DiscountedTechnologyEmissionsPenalty[r,t,y] - DiscountedSalvageValue[r,t,y]) +
            sum{r in REGION, s in STORAGE, y in YEAR} (CapitalCostStorage[r,s,y] * NewStorageCapacity[r,s,y] / DiscountFactorStorage[r,s,y] -
            - SalvageValueStorage[r,s,y] / ((1+DiscountRateStorage[r,s])^(max{yy in YEAR} max(yy)-min{yy in YEAR} min(yy)+1)));

fast

minimize cost: sum{r in REGION, t in TECHNOLOGY, y in YEAR}
    ((((sum{yy in YEAR: y-yy < OperationalLife[r,t] && y-yy>=0}
        NewCapacity[r,t,yy])+ ResidualCapacity[r,t,y]) * FixedCost[r,t,y] +
        sum{m in MODEperTECHNOLOGY[t], l in TIMESLICE}
            RateOfActivity[r,l,t,m,y] * YearSplit[l,y] * VariableCost[r,t,m,y])  /  DiscountFactorMid[r,y] +
            CapitalCost[r,t,y] * NewCapacity[r,t,y] * CapitalRecoveryFactor[r,t] * PvAnnuity[r,t]  /  DiscountFactor[r,y] +
            DiscountedTechnologyEmissionsPenalty[r,t,y] - DiscountedSalvageValue[r,t,y]) +
            sum{r in REGION, s in STORAGE, y in YEAR} (CapitalCostStorage[r,s,y]  *  NewStorageCapacity[r,s,y]  /  DiscountFactorStorage[r,s,y] -
            - SalvageValueStorage[r,s,y]  /  ((1+DiscountRateStorage[r,s])^(max{yy in YEAR} max(yy)-min{yy in YEAR} min(yy)+1)));

To be honest, none of these objective functions are correct. The ECEMF correctly includes a technology-specific version of DiscountFactorMid but omits the correction to storage which is in fast and osembe. osembe and fast include the updated capacity recovery factor * pv annuity, but omit a technology-specific discountfactor.

[willu47]

I've create #102 to update the DiscountFactor and DiscountFactorMid parameters to technology-specific versions (FYI @robertodawid)

[HauHe]

Thanks @willu47 for taking the time to have a look at this!
The objective function osembe and fast should be, and from what I see are, the same.
When comparing the code I see that in the ecemf version (the older osemosys version) the equation for the CapitalRecoveryFactor and CapitalRecoveryFactorMid are almost the same as the ones of the DiscountFactor and DiscountFactorMid in the osembe and fast. However, in the former the DiscountRate is defined per region and technology, while in the latter only per region.

[HauHe]

One thing I missed to mention in my description of the issue above, I removed the non-negativity condition for DiscountedTechnologyEmissionsPenalty and ModelPeriodEmissions, however i did so in both codes.

[willu47]

@robertodawid's appendix here describes the implementation.

Discount factor is a model generic value and the CRF is specific to the tech:
In short: The general or global DR is used to a) discount fixed and variable operating costs and b) discount investment cost and salvage value, then calculate the present value of the total system cost.
The individual discount rate is only used to calculate the annual payments (CRF). it is also used to calculate the alvage value. All these payments (annuity) are discounted back to the investment year using the global discount rate.

[robertodawid]

@HauHe, I think there is confusion in the ecemf model. The capital recovery factor allows the calculation of a series of equal payments (annuity), and the discount rate is used to discount cash flows to their present value. It would be a good idea to update the naming from CapitalRecoveryFactor to DiscountFactor.

[HauHe]

@robertodawid, I fully agree that the naming in the latest version of the osemosys code makes more sense.

However, this is not the point. If you have a look into the code that I shared, you will see that in the old code the equation for the CapitalRecoveryFactor are exactly the equations of the DiscountFactor in the latest code. Changing the name of a parameter won't change the results.

[robertodawid]

Hi, I have an error with Otoole (dev version) that led me to check the formulation to calculate the Capital Investment (result_package.py). Specifically, the lines that assign the value of a discount rate (lines: 295-334)

if "DiscountRateIdv" in self.keys():
   discount_rate = self["DiscountRateIdv"]
 else:
   discount_rate = self["DiscountRate"]

According to the actual formulation, it checks if there is or is not an Individual Discount rate (DiscountRateIdv), and a discount rate is assigned to disocunt_rate. Then, it calculates the capital investment and other factors like the CRF and PVannuity. This is okay if there is no DiscountRateIdv, the global discount rate is used as a discount_rate . However, this is not the case when having Individual Discount Rates, the 'DiscountRateIdv' is used ONLY in calculating the CapitalRecoveryFactor (CRF).
With the current formulation, if there is DiscountRateIdv, this rate will be used to calculate the other factors, such as the present value of an annuity, which is incorrect. Happy to talk about this @trevorb1

[trevorb1]

Hello @robertodawid !! 🎉🎉🎉

I believe there is an active PR in the otoole repository right now that may address this; this comment specifically goes over what the PR covers (given that the scope of the PR expanded a bit and the name is a little outdated). Included in this are fixing issues surrounding how DiscountRate and DiscountRateIdv are handled in the otoole result calculations.

Apologies for not merging the PR yet; reviewing it and getting it merged has been on my todo list for a while now. Would you be able to check the code updates in the PR (ie. the write-defaults branch) to see if this addresses your comment?

(Side note; not sure that this otoole PR will address the issue described in this issue ticket!)