iiasa / message_ix

The integrated assessment and energy systems model MESSAGEix
https://docs.messageix.org
Apache License 2.0
116 stars 152 forks source link

COST_NODAL calculation error #517

Open zsalimian opened 2 years ago

zsalimian commented 2 years ago

Dear Sir,

I have developed a simple code with 2 technologies and one resource. Please see the results of model:

CAP: node_loc technology year_vtg year_act lvl mrg 0 World cycled_ppl 2010 2020 273.000000 0.000000e+00 1 World cycled_ppl 2010 2030 273.000000 0.000000e+00 2 World cycled_ppl 2020 2020 156.629630 0.000000e+00 3 World cycled_ppl 2020 2030 156.629630 0.000000e+00 4 World cycled_ppl 2020 2040 156.629630 0.000000e+00 5 World cycled_ppl 2030 2030 109.629630 0.000000e+00 6 World cycled_ppl 2030 2040 109.629630 0.000000e+00 7 World cycled_ppl 2040 2040 307.074074 0.000000e+00 8 World wind_ppl 2010 2020 0.000000 4.940656e-324 9 World wind_ppl 2010 2030 0.000000 4.940656e-324 10 World wind_ppl 2010 2040 0.000000 4.940656e-324 11 World wind_ppl 2020 2020 0.000000 6.353000e+02 12 World wind_ppl 2020 2030 0.000000 3.545000e+02 13 World wind_ppl 2020 2040 0.000000 4.940656e-324 14 World wind_ppl 2030 2030 0.000000 6.353000e+02 15 World wind_ppl 2030 2040 0.000000 3.545000e+02 16 World wind_ppl 2040 2040 0.000000 3.957000e+02

New_Cap

node_loc technology year_vtg lvl mrg 0 World cycled_ppl 2020 15.662963 0.000000e+00 1 World cycled_ppl 2030 10.962963 0.000000e+00 2 World cycled_ppl 2040 30.707407 0.000000e+00 3 World wind_ppl 2020 0.000000 4.940656e-324 4 World wind_ppl 2030 0.000000 4.940656e-324 5 World wind_ppl 2040 0.000000 0.000000e+00

The answers are correct. I calcuated manually and it has estimated correct CAP and CAP_new. But Cost_Nodal is calculated incorrectly. For instance, for year 2020 the results show that it constructed 15.66 MW. The cost of each MW is 632$.The total investment cost would be=9931$ The total capacity cost in 2020 is around 430 and the fixed cost is 10.32. The total fixed cost would be=4433$ The whole Cost_Nodal must be 14370$ in 2020. But as you see below, the model does not consider investment cost in results: node year lvl mrg 0 World 2020 4433.777778 0.0 1 World 2030 5565.155556 0.0 2 World 2040 5916.800000 0.0

Why this has happened?

###########################################

import pandas as pd
import ixmp
import message_ix
import matplotlib.pyplot as plt

from message_ix.utils import make_df

mp = ixmp.Platform()
scenario = message_ix.Scenario(mp, model='World',
                               scenario='baseline', version='new')
history = [2010]
model_horizon = [2020, 2030, 2040]
scenario.add_horizon(
    year=history + model_horizon,
    firstmodelyear=model_horizon[0],
)
year_df = scenario.vintage_and_active_years()
vintage_years, act_years = year_df['year_vtg'], year_df['year_act']
country = 'World'
scenario.add_spatial_sets({'country': country})
scenario.add_set("commodity", ["gas","wind","electricity"])
scenario.add_set("level", ["primary","primary_wind","final"])
scenario.add_set("emission", ["CO","CO2","CH4","N2O","NOx","SO2","PM"])
scenario.add_set("technology", ['cycled_ppl','wind_ppl'])
scenario.add_set("mode", "standard")
scenario.add_set("time","year")
scenario.add_set("type_emission",["greenhouse_gases","pollutans"])
scenario.add_set("grade", "1")
scenario.add_set("lvl_spatial", "World")
scenario.add_set("lvl_temporal", "year")
scenario.add_set("rating", ["firm","unrated"])
scenario.add_set("level_resource", "primary")
scenario.add_set("level_renewable", "primary_wind")

#capacity factor########################################################################################
base_capacity_factor = {
    'node_loc': country,
    'year_vtg': vintage_years,
    'year_act': act_years,
    'time': 'year',
    'unit': '-',
}

capacity_factor = {
    'cycled_ppl': 1,
    'wind_ppl': 1,
    # 'grid':1
}

for tec, val in capacity_factor.items():
    df = make_df(base_capacity_factor, technology=tec, value=val)
    scenario.add_par('capacity_factor', df)
##########################################################################################

base_operation_factor = {
    'node_loc': country,
    'year_vtg': vintage_years,
    'year_act': act_years,
    'time': 'year',
    'unit': '-',
}
operation_factor = {
    'cycled_ppl': 0.675,
    'wind_ppl': 0.46,
    # 'grid':0.37
}

for tec, val in operation_factor.items():
    df = make_df(base_operation_factor, technology=tec, value=val)
    scenario.add_par('operation_factor', df)

##############################################################################################

base_inv_cost = {
    'node_loc': country,
    'year_vtg': model_horizon,
    'unit': 'USD/kW',
}

# Adding a new unit to the library
mp.add_unit('USD/kW')

costs = {
    'cycled_ppl': 632,
    'wind_ppl':[847.35, 820.85, 807.6],
    # 'grid':152.4,
}
################################################################################################3
fix_cost_pd = pd.DataFrame({
        'node_loc': country,
        'technology': 'cycled_ppl',
        'year_vtg':vintage_years,
        'year_act':act_years,
        'value': 10.32,
        'unit': '-',
    })

fix_cost_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("fix_cost", fix_cost_pd)

base_fix_cost = {
    'node_loc': country,
    'year_vtg': vintage_years,
    'year_act': act_years,
    'unit': 'USD/kWa',
}

costs = {

    'wind_ppl':[63.53,35.45,39.57,63.53,35.45,39.57,63.53,35.45,39.57],

}

for tec, val in costs.items():
    df = make_df(base_fix_cost, technology=tec, value=val)
    scenario.add_par('fix_cost', df)

#################################################################################################

base_var_cost = {
    'node_loc': country,
    'year_vtg': vintage_years,
    'year_act': act_years,
    'mode': 'standard',
    'time': 'year',
    'unit': 'USD/kWa',
}
# in $ / kWa (costs associatied to the degradation of equipment when the plant is functioning
# per unit of energy produced kW·year = 8760 kWh.
# Therefore this costs represents USD per 8760 kWh of energy). Do not confuse with fixed O&M units.

costs = {
    'cycled_ppl': 0,
    'wind_ppl': 0,
}

for tec, val in costs.items():
    df = make_df(base_var_cost, technology=tec, value=val)
    scenario.add_par('var_cost', df)

#################################################################################################

base = {
    'node_loc': country,
    'year_vtg': vintage_years,
    'year_act': act_years,
    'mode': 'standard',
    'time': 'year',
    'unit': '-',
}

base_input = make_df(base, node_origin=country, time_origin='year')
base_output = make_df(base, node_dest=country, time_dest='year')

cycled_out = make_df(base_output, technology='cycled_ppl', commodity='electricity',
                   level='final', value=1.)
scenario.add_par('output', cycled_out)

cycled_in = make_df(base_input, technology='cycled_ppl', commodity='gas',
                       level='primary', value=1.77)
scenario.add_par('input', cycled_in)

wind_out = make_df(base_output, technology='wind_ppl', commodity='electricity',
                    level='final', value=1.)
scenario.add_par('output', wind_out)

wind_in = make_df(base_input, technology='wind_ppl', commodity='wind',
                    level='primary_wind', value=1.)
scenario.add_par('input', wind_in)

################################################################################################

base_growth = {
    'node_loc': country,
    'year_act': model_horizon,
    'time': 'year',
    'unit': '-',
}

growth_activity_up = {
    'cycled_ppl': 0.1,
    'wind_ppl': 0.1,
    # 'grid':0.1,
}

for tec,val in growth_activity_up.items():
    df = make_df(base_growth, technology=tec, value=val)
    scenario.add_par('growth_activity_up', df)

################################################################################################
duration_period_pd = pd.DataFrame({
        # 'node': country,
        # 'commodity': 'gas',
        # 'grade':'1',
        'year': model_horizon,
        'value': 10,
        'unit': '-',
    })

duration_period_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("duration_period", duration_period_pd)

#################################################################################################
base_duration_time={

    'year': model_horizon,
    'unit': '-',
}

duration_t = {
    'year': 1,
    # 'seasons': 0.25,
    # 'day': 1/365,
}
#
for t, v in duration_t.items():
    df = make_df(base_duration_time, time=t, value=v)
    scenario.add_par('duration_time', df)

#######################################################################################################

base_activity = {
    'node_loc': country,
    'year_act': history,
    'mode': 'standard',
    'time': 'year',
    'unit': 'GWa',
}
old_activity = {
    'cycled_ppl':  184.182,
    # 'grid':184.182,
    'wind_ppl':0,
}

for tec, val in old_activity.items():
    df = make_df(base_activity, technology=tec, value=val)
    scenario.add_par('historical_activity', df)
############################################################################################################
base_historical_capacity = {
    'node_loc': country,
    'year_vtg': history,
    'unit': 'GW',
}

old_capacity = {
     'cycled_ppl':27.3,
     # 'grid':27.3,
     'wind_ppl':0,
}

for tec, val in old_capacity.items():
    df = make_df(base_historical_capacity, technology=tec, value=val)
    scenario.add_par('historical_new_capacity', df)
############################################################################################################
scenario.add_par("interestrate", model_horizon, value=0, unit='-')
#######################################################################################################
construction_time_pd = pd.DataFrame({
        'node_loc': country,
        'technology': 'cycled_ppl',
        'year_vtg': vintage_years,
        'value': 1,
        'unit': '-',
    })
# #.round()
construction_time_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("construction_time", construction_time_pd)
########################################################################################################

########################################################################################################
bound_extraction_up_pd = pd.DataFrame({
        'node': country,
        'commodity': 'gas',
        'grade':'1',
        'year': model_horizon,
        'value': 100000000000,
        'unit': '-',
    })

bound_extraction_up_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("bound_extraction_up", bound_extraction_up_pd)

##########################################################################################################
resource_remaining_pd = pd.DataFrame({
        'node': country,
        'commodity': 'gas',
        'grade':'1',
        'year': model_horizon,
        'value': 0.5,
        'unit': '-',
    })

resource_remaining_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("resource_remaining", resource_remaining_pd)
#######################################################################################
resource_volume_pd = pd.DataFrame({
        'node': country,
        'commodity': 'gas',
        'grade':'1',
        'year': history,
        'value': 100000000000,
        'unit': '-',
    })

resource_volume_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("resource_volume", resource_volume_pd)

###############################################################################################
renewable_potential_pd = pd.DataFrame({
        'node': country,
        'commodity': 'wind',
        'grade':'1',
        'level': 'primary_wind',
        'year': model_horizon,

        'value': 1000000000000,
        'unit': 'GWa',
    })

renewable_potential_pd
# We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("renewable_potential", renewable_potential_pd)

################################################################################################

###########################################################################################################
base_min_utilization_factor = {
    'node_loc': country,
    'year_vtg': vintage_years,
    'year_act': act_years,
    'time': 'year',
    'unit': '-',
}
min_utilization_factor = {
    'cycled_ppl': 0.1,
    'wind_ppl': 0.1,
    # 'grid':0.1
}

for tec, val in min_utilization_factor.items():
    df = make_df(base_min_utilization_factor, technology=tec, value=val)
    scenario.add_par('min_utilization_factor', df)

#############################################################################################################
renewable_capacity_factor_pd = pd.DataFrame({
        'node': country,
        'commodity':'wind',
        'grade':'1',
        'level':'primary_wind',
        'year': model_horizon,
        'value':1,
        'unit': '-',
    })
#.round()
renewable_capacity_factor_pd
# We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("renewable_capacity_factor", renewable_capacity_factor_pd)
############################################################################################################

base_initial_new_capacity_up={
        'node_loc': country,
        'year_vtg': vintage_years,
        'unit': '-',
}

initial_new_capacity_up={
        'wind_ppl': 10,
        'cycled_ppl': 10,
        # 'grid': 10,
}

for tec,val in initial_new_capacity_up.items():
    df=make_df(base_initial_new_capacity_up,technology=tec,value=val)
    scenario.add_par('initial_new_capacity_up',df)

#############################################################################################################

base_growth_new_capacity_up= {
    'node_loc': country,
    'year_vtg': vintage_years,
    'unit': '-',
}

growth_new_capacity_up= {
    'wind_ppl': 0.1,
    'cycled_ppl': 0.1,
    # 'grid': 0.1,
}

for tec,val in growth_new_capacity_up.items():
    df=make_df(base_growth_new_capacity_up,technology=tec,value=val)
    scenario.add_par('growth_new_capacity_up',df)

###########################################################################################################
base_bound_activity_up={
        'node_loc': country,
        'year_act': act_years,
        'mode': 'standard',
        'time':'year',
        'unit': '-',
}

bound_activity_up={
        'wind_ppl': 10000000000000000,
        'cycled_ppl': 10000000000000000,
        # 'grid': 10000000000000000,
}

for tec,val in bound_activity_up.items():
    df=make_df(base_bound_activity_up,technology=tec,value=val)
    scenario.add_par('bound_activity_up',df)

##############################################################################################################
technical_lifetime1_pd = pd.DataFrame({
        'node_loc': country,
        'year_vtg': vintage_years,
        'technology':'cycled_ppl',
        'value': 30,
        'unit': '-',
    })
# #.round()
technical_lifetime1_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("technical_lifetime", technical_lifetime1_pd)

#############################################################################################################
technical_lifetime2_pd = pd.DataFrame({
        'node_loc': country,
        'year_vtg': vintage_years,
        'technology':'wind_ppl',
        'value': 20,
        'unit': '-',
    })
# #.round()
technical_lifetime2_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("technical_lifetime", technical_lifetime2_pd)

################################################################################################################
bound_total_capacity_up_pd = pd.DataFrame({
        'node_loc': country,
        'year_act': act_years,
        'technology':'cycled_ppl',
        'value': 100000000000000000,
        'unit': '-',
    })
bound_total_capacity_up_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("bound_total_capacity_up", bound_total_capacity_up_pd)

#####################################################################################################
bound_total_capacity_up_pd = pd.DataFrame({
        'node_loc': country,
        'year_act': act_years,
        'technology':'wind_ppl',
        'value': 100000000000000000,
        'unit': '-',
    })

bound_total_capacity_up_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("bound_total_capacity_up", bound_total_capacity_up_pd)
############################################################################################################
bound_new_capacity_up_pd = pd.DataFrame({
        'node_loc': country,
        'year_vtg': vintage_years,
        'technology':'wind_ppl',
        'value': 100000000000000000,
        'unit': '-',
    })

bound_new_capacity_up_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("bound_new_capacity_up", bound_new_capacity_up_pd)

############################################################################################################
bound_new_capacity_up_pd = pd.DataFrame({
        'node_loc': country,
        'year_vtg': vintage_years,
        'technology':'cycled_ppl',
        'value': 100000000000000000,
        'unit': '-',
    })

bound_new_capacity_up_pd
# # We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("bound_new_capacity_up", bound_new_capacity_up_pd)

##############################################################################################################
#
demand_profile = pd.Series([290.345, 363.56, 386.78],
                        index=pd.Index(model_horizon, name='Time'))
demand_profile.plot(title='Demand')

electricity_demand = pd.DataFrame({
        'node': country,
        'commodity': 'electricity',
        'level': 'final',
        'year': model_horizon,
        'time': 'year',
        'value': (demand_profile).round(),
        'unit': 'GWa',
    })

electricity_demand
# We use add_par for adding data to a MESSAGEix parameter
scenario.add_par("demand", electricity_demand)

################################################################################################################
from message_ix import log
log.info('version number prior to commit: {}'.format(scenario.version))

scenario.commit(comment='basic model of Westeros electrification')

log.info('version number prior committing to the database: {}'.format(scenario.version))
scenario.set_as_default()
scenario.solve()
scenario.var('OBJ')['lvl']
#df = scenario.var('PRICE_COMMODITY')
LauWien commented 2 years ago

Dear @zsalimian, thank you for reaching out. For better readability, please be so kind to mark code with triple backticks ```, described here.