search

Pyomo / pyomo

An object-oriented algebraic modeling language in Python for structured optimization problems.
https://www.pyomo.org
Other
1.92k stars 494 forks source link

List index out of range error for GAMS solutions #2955

Open gidden opened 11 months ago

gidden commented 11 months ago

Summary

We run opt.solve(model) with a gams SolverFactory and see a list index out of range error. This is due to pyomo splitting lines in the GAMS.py solver file and expecting 3 elements separated by spaces, but only two exist.

Steps to reproduce the issue

Using

data.xlsx

from pyomo.environ import *
import pandas as pd
import matplotlib.pyplot as plt

# BASIC "TWO-FACTOR" PARAMETERS
tfl_params = {
    'b':{
            'Nuclear':.061,
            'Coal':.355,
            'Gas T':.275,
            'Gas CC':.031,
            'Bioenergy (Total)':.008,
            'Geothermal':-.607,
            'Hydro':-.527,
            'Wind Onshore':.039,
            'Wind Offshore':-.273,
            'Solar PV (Utility)':.414,
            'Solar CSP':-0.017
    },
    'rho':{
            'Nuclear':.435,
            'Coal':.844,
            'Gas T':.734,
            'Gas CC':.561,
            'Bioenergy (Total)':.571,
            'Geothermal':.859,
            'Hydro':1.181,
            'Wind Onshore':.773,
            'Wind Offshore':.121,
            'Solar PV (Utility)':.676,
            'Solar CSP':-0.346
    },
    'alpha':{
            'Nuclear':0.115,
            'Coal':2.536,
            'Gas T':.103,
            'Gas CC':.237,
            'Bioenergy (Total)':.056,
            'Geothermal':1.760,
            'Hydro':0.240,
            'Wind Onshore':.704,
            'Wind Offshore':.376,
            'Solar PV (Utility)':.428,
            'Solar CSP':.188
    },
    'nbr_init':{
            'Nuclear':770,
            'Coal':2161,
            'Gas T':3941,
            'Gas CC':3039,
            'Bioenergy (Total)':1436,
            'Geothermal':403,
            'Hydro':1275,
            'Wind Onshore':404073,
            'Wind Offshore':7191,
            'Solar PV (Utility)':11381,
            'Solar CSP':70
    },
    'unt_sizes':{
            ('Nuclear','s'): 0.6,
            ('Nuclear','m'): 1.1,
            ('Nuclear','l'): 1.5,
            ('Coal','s'): 0.45,
            ('Coal','m'): 0.9,
            ('Coal','l'): 1.2,
            ('Gas T','s'): 0.2,
            ('Gas T','m'): 0.5,
            ('Gas T','l'): 1.0,
            ('Gas CC','s'): 0.45,
            ('Gas CC','m'): 0.9,
            ('Gas CC','l'): 1.2,
            ('Bioenergy (Total)','s'): 0.045,
            ('Bioenergy (Total)','m'): 0.1,
            ('Bioenergy (Total)','l'): 0.5,
            ('Geothermal','s'): 0.05,
            ('Geothermal','m'): 0.1,
            ('Geothermal','l'): 0.5,
            ('Hydro','s'): 0.45,
            ('Hydro','m'): 0.9,
            ('Hydro','l'): 1.2,
            ('Wind Onshore','s'): 0.003,
            ('Wind Onshore','m'): 0.005,
            ('Wind Onshore','l'): 0.010,
            ('Wind Offshore','s'): 0.007,
            ('Wind Offshore','m'): 0.015,
            ('Wind Offshore','l'): 0.025,
            ('Solar CSP','s'): 0.05,
            ('Solar CSP','m'): 0.1,
            ('Solar CSP','l'): 0.5,
            ('Solar PV (Utility)','s'): 0.06,
            ('Solar PV (Utility)','m'): 0.2,
            ('Solar PV (Utility)','l'): 1.0,
    },
    'size_ref':{##
            'Nuclear':1.1,
            'Coal':0.590,
            'Gas T':0.154,
            'Gas CC':0.467,
            'Bioenergy (Total)':.049,
            'Geothermal':0.067,
            'Hydro':.455,
            'Wind Onshore':0.003,
            'Wind Offshore':0.007,
            'Solar PV (Utility)':0.060,
            'Solar CSP':0.084
    },
}

# READ DATA FROM EXCEL
sheet = ['Capacity',
        ]
df_newcapacity = pd.read_excel('./data.xlsx', 
                   sheet_name=sheet[0], header=[0], 
                   usecols='A:N', engine='openpyxl').fillna(0)
df_newcapacity = df_newcapacity.set_index('Technology')
df = df_newcapacity.copy()
df_bin = df.applymap(lambda x: 1 if x > 0 else 0)

# LINEARIZED LEARNING AND SCALE-UP FACTOR
## EoS factor calculation
eos_factor = {}
for t in list(df.index):
    for s in ['s','m','l']:
        eos_factor.update({(t,s): (tfl_params.get('unt_sizes').get((t,s))
                                   /tfl_params.get('size_ref').get(t))**(tfl_params.get('rho').get(t))})

## Linearized learning calculation
xnbr_values  = [*range(0,31)]
xnbr_segment = [*range(1,31)]

xnbr_unit = {}
yinv_cost = {}
minv_cost = {}
cinv_cost = {}
yu_avg = {}
mu_avg = {}
cu_avg = {}

## Xnbr and Yinv values
for t in list(df.index):
    for x in xnbr_values:
        xnbr_unit.update({(t,x): 2**(x)})
        yinv_cost.update({(t,x): (((tfl_params.get('nbr_init').get(t)
                                   +(2**(x)))/tfl_params.get('nbr_init').get(t))**(-tfl_params.get('b').get(t)))
                         })
        yu_avg.update({(t,x): tfl_params.get('size_ref').get(t) * (((tfl_params.get('nbr_init').get(t)
                                   +(2**(x)))/tfl_params.get('nbr_init').get(t))**(tfl_params.get('alpha').get(t)))
                         })

## Linearized parameters
for t in list(df.index):
    for x in xnbr_segment:
        # values calculation
        minv_cost_value = ((yinv_cost.get((t,x)) - yinv_cost.get((t,x-1)))
                                  /(xnbr_unit.get((t,x)) - xnbr_unit.get((t,x-1))))
        cinv_cost_value = minv_cost_value * xnbr_unit.get((t,x-1)) + yinv_cost.get((t,x-1))

        mu_avg_value    = ((yu_avg.get((t,x)) - yu_avg.get((t,x-1)))
                                  /(xnbr_unit.get((t,x)) - xnbr_unit.get((t,x-1))))
        cu_avg_value    = mu_avg_value * xnbr_unit.get((t,x-1)) + yu_avg.get((t,x-1))

        # updating the dictionary
        minv_cost.update({(t,x): minv_cost_value})
        mu_avg.update({(t,x): mu_avg_value})
        cinv_cost.update({(t,x): cinv_cost_value})
        cu_avg.update({(t,x): cu_avg_value})

# Sizes portfolio optimization
## Declare model's name
model   = ConcreteModel('Optimal separated effects')

## Define sets ##
model.t  = Set(initialize=list(df.index), doc='Technology')
model.s  = Set(initialize=['s','m','l'], doc='Unit sizes')
model.y  = Set(initialize=list(df.columns), doc='Periods')
model.yy = Set(initialize=list(df.columns), doc='Periods')
model.x  = Set(initialize=xnbr_segment, doc='Linearized-x segment')
#map(str, list(df.columns))

## Define data and parameters ## 
model.nbr_unit_ref = Param(model.t, initialize=tfl_params.get('nbr_init'), doc='Units # set as reference')
model.unt_size_ref = Param(model.t, initialize=tfl_params.get('size_ref'), doc='Units size set as reference')
model.unt_sizes    = Param(model.t, model.s, initialize=tfl_params.get('unt_sizes'), doc='Units size data')
model.learning_par = Param(model.t, initialize=tfl_params.get('b'), doc='Units # data')
model.eos_par      = Param(model.t, initialize=tfl_params.get('rho'), doc='Units # data')
model.eos_f        = Param(model.t, model.s, initialize=eos_factor, doc='Units # data')
model.scaleup_par  = Param(model.t, initialize=tfl_params.get('alpha'), doc='Units # data')
model.new_capacity = Param(model.t, model.y, initialize=df.stack().to_dict(), doc='Units # data')
model.bin_new_cap  = Param(model.t, model.y, initialize=df_bin.stack().to_dict(), doc='Units # data')
model.minv         = Param(model.t, model.x, initialize=minv_cost, doc='Units # data')
model.mu           = Param(model.t, model.x, initialize=mu_avg, doc='Units # data')
model.cinv         = Param(model.t, model.x, initialize=cinv_cost, doc='Units # data')
model.cu           = Param(model.t, model.x, initialize=cu_avg, doc='Units # data')

## Define variables ##
model.CAPEX_TEC_IDX  = Var(model.t, model.y, bounds=(0.0,None), doc='Investment cost')
model.CAPEX_TEC_IDXr = Var(model.t, model.y, bounds=(0.0,None), doc='Investment cost index')
model.NBR_UNIT       = Var(model.t, model.s, model.y, bounds=(0.0,None), doc='Number of units')

### FROM HERE
## Define constraints ##
def capacitybalance_rule(model, t, y):
  return model.new_capacity[t,y] <= sum(model.NBR_UNIT[t,s,y] * model.unt_sizes[t,s] for s in model.s)
model.capacitybalance = Constraint(model.t, model.y, rule=capacitybalance_rule, doc='Absolute x calculation')

def costlearn_rule(model, t, y, x):
  return model.CAPEX_TEC_IDXr[t,y] >= ((sum(model.NBR_UNIT[t,s,yy] for s in model.s for yy in model.yy if yy <= y)
                                        * model.minv[t,x]) + model.cinv[t,x])
model.costlearn = Constraint(model.t, model.y, model.x, rule=costlearn_rule, doc='Cost reduction from the learning effect')

def costhold_rule(model, t, y):
  if y == model.y.first():
    return model.CAPEX_TEC_IDX[t,y] == (model.bin_new_cap[t,y] * model.CAPEX_TEC_IDX[t,y] 
                                            + (1-model.bin_new_cap[t,y]) * 1)
  else:
    return model.CAPEX_TEC_IDX[t,y] == (model.bin_new_cap[t,y] * model.CAPEX_TEC_IDX[t,y] 
                                            + (1-model.bin_new_cap[t,y]) * model.CAPEX_TEC_IDX[t,model.y.prev(y)])
model.costhold = Constraint(model.t, model.y, rule=costhold_rule, doc='Cost hold at previous period level if no new capacity')

def costscale_rule(model, t, y):
  return (model.CAPEX_TEC_IDX[t,y] * model.new_capacity[t,y]) >= (sum(model.NBR_UNIT[t,s,y] * model.eos_f[t,s] for s in model.s)
                                                                  * model.CAPEX_TEC_IDXr[t,y]
                                                                  * model.unt_size_ref[t])
model.costscale = Constraint(model.t, model.y, rule=costscale_rule, doc='cost weight average calculation')

## Define objective functions ##
def objective_rule(model):
  return sum(model.CAPEX_TEC_IDX[t,y] * model.new_capacity[t,y] for t in model.t for y in model.y)
model.objective = Objective(rule=objective_rule, sense=minimize, doc='Objective function in squared absolute percentage term')
#  return sum(model.CAPEX_TEC_IDX[t,y] * model.new_capacity[t,y] for t in model.t for y in model.y)

if __name__ == '__main__':
    # This emulates what the pyomo command-line tools does
    from pyomo.opt import SolverFactory
    import pyomo.environ
    opt = SolverFactory("gams")
    opt.options['solver'] = 'minos'
    results = opt.solve(model)
    model.display()

obj_value = model.objective(value)
print('The objective value is', obj_value)

df_capex_tec_idx=[]
fig, ax = plt.subplots(figsize=(10,7), sharex=True)
for t in list(df.index):
    capex_tec_idx = []
    for y in list(df.columns):
        capex_tec_idx.append(model.CAPEX_TEC_IDX[t,y](value))
    df_capex_tec_idx.append(capex_tec_idx)
    plot_1 = ax.plot(list(df.columns), capex_tec_idx,
                     ls='-', mfc='white', label=t, zorder=2)
plt.show()

Error Message

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [3], in <cell line: 66>()
     70     opt = SolverFactory("gams")
     71     opt.options['solver'] = 'minos'
---> 72     results = opt.solve(model)
     73     model.display()
     75 obj_value = model.objective(value)

File ~\Anaconda3\lib\site-packages\pyomo\solvers\plugins\solvers\GAMS.py:853, in GAMSShell.solve(self, *args, **kwds)
    850         model_soln, stat_vars = self._parse_gdx_results(
    851             results_filename, statresults_filename)
    852     else:
--> 853         model_soln, stat_vars = self._parse_dat_results(
    854             results_filename, statresults_filename)
    855 finally:
    856     if not keepfiles:

File ~\Anaconda3\lib\site-packages\pyomo\solvers\plugins\solvers\GAMS.py:1236, in GAMSShell._parse_dat_results(self, results_filename, statresults_filename)
   1233 for line in results_text.splitlines()[1:]:
   1234     # line = line.replace('GAMS_OBJECTIVE', 'GAMS_OBJECTIVE ')
   1235     items = line.split()
-> 1236     model_soln[items[0]] = (items[1], items[2])
   1238 return model_soln, stat_vars

IndexError: list index out of range

Information on your system

Pyomo version: 6.4.2 Python version: 3.9.12 Operating system: windows How Pyomo was installed (PyPI, conda, source): conda Solver (if applicable): GAMS using minos or ipopt (observed in both)

Additional information

We were able to address this by patching GAMS.py line 1234 with

line = line.replace("GAMS_OBJECTIVE", "GAMS_OBJECTIVE ")

effectively forcing the 3 elements. Of course a complete hack, not a fix at all - just to show that a solution is actually found.

gidden commented 11 months ago

cc @ywpratama