I'm new to both Python and Pyomo. Nonetheless, I've made the below code following my previous question. However, if I run it using glpk solver, it says
RuntimeError: Cannot write legal LP file. Objective 'cost' has nonlinear terms that are not quadratic.
I think the formulation is an LP. I'm guessing it is the calculation of parameter "demandsperyear"? But, it shouldn't be the cause as it is a parameter, not a variable, right?
If i use cbc, I've got errors that I cannot understand.
from pyomo.environ import *
m = AbstractModel()
m.techs = Set()
m.products = Set()
m.suppliers = Set()
m.demands = Set()
m.years = Set()
m.initial_demands = Param(m.demands,m.products)
m.growths = Param(m.demands,m.products)
m.yields = Param(m.techs,m.products,m.suppliers)
m.CO2emissions = Param(m.techs,m.suppliers)
m.TAC = Param(m.techs,m.suppliers)
m.max_installed = Param(m.techs,m.suppliers)
m.installed_techs = Var(m.techs, m.suppliers, m.years, within=Binary)
m.number_installed_techs= Var(m.techs, m.suppliers, m.years, bounds=(0,None))
m.totalTAC = Var(bounds=(None,None))
m.totalCO2 = Var(bounds=(None,None))
m.totalTACperyear = Var(m.years, bounds=(None, None))
m.totalCO2peryear = Var(m.years, bounds=(None, None))
def total_TAC(m):
return m.totalTAC == sum(m.totalTACperyear[year] for year in m.years)
m.cost = Objective(rule=total_TAC, sense=minimize)
def total_CO2(m):
return m.totalCO2 == sum(m.totalCO2peryear[year] for year in m.years)
m.totalCO2_cons = Constraint(rule=total_CO2)
def total_TACperyear(m,year):
return m.totalTACperyear[year] == sum(m.number_installed_techs[tech,supplier,year]*m.TAC[tech,supplier] for tech in m.techs for supplier in m.suppliers)
m.TACperyear_cons = Constraint(m.years, rule=total_TACperyear)
def total_CO2peryear(m,year):
return m.totalCO2peryear[year] == sum(m.number_installed_techs[tech,supplier,year]*m.CO2emissions[tech,supplier] for tech in m.techs for supplier in m.suppliers)
m.CO2peryear_cons = Constraint(m.years, rule=total_CO2peryear)
def demandperyear(m,demand,product,year):
if year == 1:
expression = m.initial_demands[demand,product]
return expression
expression = m.initial_demands[demand,product]*(1 + m.growths[demand,product]*year)
return expression
m.demandsperyear = Param(m.demands, m.products,m.years, initialize=demandperyear)
def fulfill_demands(m,supplier,demand,product,year):
if supplier == demand:
expression = sum(m.number_installed_techs[tech,supplier,year]*m.yields[tech,product,supplier] for tech in m.techs) >= m.demandsperyear[demand,product,year]
return expression
return Constraint.Skip
m.fulfill_demands_cons = Constraint(m.suppliers,m.demands,m.products,m.years, rule=fulfill_demands)
def max_installed(m,tech,supplier):
expression = sum(m.number_installed_techs[tech,supplier,year] for year in m.years) <= m.max_installed[tech,supplier]
return expression
m.max_installed_cons = Constraint(m.techs,m.suppliers, rule=max_installed)
def accum_installed(m,tech,supplier,year):
if year == 1:
expression = m.number_installed_techs[tech,supplier,year] == m.installed_techs[tech,supplier,year]
return expression
expression = m.number_installed_techs[tech,supplier,year] == m.number_installed_techs[tech,supplier,year-1] + m.installed_techs[tech,supplier,year]
return expression
m.accum_cons = Constraint(m.techs,m.suppliers,m.years, rule=accum_installed)
# parameter
data = {None: {
'techs' : {None: ['tek_A', 'tek_B', 'tek_C']},
'products' : {None: ['prod_A', 'prod_B', 'prod_C']},
'suppliers' : {None: ['daerah_A', 'daerah_B', 'daerah_C']},
'demands' : {None: ['daerah_A', 'daerah_B', 'daerah_C']},
'years' : {None: [1, 2, 3, 4, 5,6,7,8,9,10]},
'initial_demands': {('daerah_A','prod_A'): 10.0, ('daerah_A','prod_B'): 10.0, ('daerah_A','prod_C'): 10.0,
('daerah_B','prod_A'): 20.0, ('daerah_B','prod_B'): 20.0, ('daerah_B','prod_C'): 20.0,
('daerah_C','prod_A'): 12.0, ('daerah_C','prod_B'): 12.0, ('daerah_C','prod_C'): 12.0},
'growths': {('daerah_A','prod_A'): 0.02, ('daerah_A','prod_B'): 0.02, ('daerah_A','prod_C'): 0.02,
('daerah_B','prod_A'): 0.02, ('daerah_B','prod_B'): 0.02, ('daerah_B','prod_C'): 0.02,
('daerah_C','prod_A'): 0.02, ('daerah_C','prod_B'): 0.02, ('daerah_C','prod_C'): 0.02},
'max_installed' : {('tek_A','daerah_A'): 100.0, ('tek_A','daerah_B'): 100.0, ('tek_A','daerah_C'): 100.0,
('tek_B','daerah_A'): 100.0, ('tek_B','daerah_B'): 100.0, ('tek_B','daerah_C'): 100.0,
('tek_C','daerah_A'): 100.0, ('tek_C','daerah_B'): 100.0, ('tek_C','daerah_C'): 100.0},
'TAC' : {('tek_A','daerah_A'): 100.0, ('tek_A','daerah_B'): 100.0, ('tek_A','daerah_C'): 100.0,
('tek_B','daerah_A'): 100.0, ('tek_B','daerah_B'): 100.0, ('tek_B','daerah_C'): 100.0,
('tek_C','daerah_A'): 100.0, ('tek_C','daerah_B'): 100.0, ('tek_C','daerah_C'): 100.0},
'CO2emissions' : {('tek_A','daerah_A'): 20.0, ('tek_A','daerah_B'): 20.0, ('tek_A','daerah_C'): 20.0,
('tek_B','daerah_A'): 20.0, ('tek_B','daerah_B'): 20.0, ('tek_B','daerah_C'): 20.0,
('tek_C','daerah_A'): 20.0, ('tek_C','daerah_B'): 20.0, ('tek_C','daerah_C'): 20.0},
'yields': {('tek_A','prod_A','daerah_A'):8.0, ('tek_A','prod_A','daerah_B'):8.0, ('tek_A','prod_A','daerah_C'):8.0,
('tek_A','prod_B','daerah_A'):8.0, ('tek_A','prod_B','daerah_B'):8.0, ('tek_A','prod_B','daerah_C'):8.0,
('tek_A','prod_C','daerah_A'):8.0, ('tek_A','prod_C','daerah_B'):8.0, ('tek_A','prod_C','daerah_C'):8.0,
('tek_B','prod_A','daerah_A'):8.0, ('tek_B','prod_A','daerah_B'):8.0, ('tek_B','prod_A','daerah_C'):8.0,
('tek_B','prod_B','daerah_A'):8.0, ('tek_B','prod_B','daerah_B'):8.0, ('tek_B','prod_B','daerah_C'):8.0,
('tek_B','prod_C','daerah_A'):8.0, ('tek_B','prod_C','daerah_B'):8.0, ('tek_B','prod_C','daerah_C'):8.0,
('tek_C','prod_A','daerah_A'):8.0, ('tek_C','prod_A','daerah_B'):8.0, ('tek_C','prod_A','daerah_C'):8.0,
('tek_C','prod_B','daerah_A'):8.0, ('tek_C','prod_B','daerah_B'):8.0, ('tek_C','prod_B','daerah_C'):8.0,
('tek_C','prod_C','daerah_A'):8.0, ('tek_C','prod_C','daerah_B'):8.0, ('tek_C','prod_C','daerah_C'):8.0}
}}
instance = m.create_instance(data)
instance.pprint()
solver = SolverFactory('glpk')
hasil = solver.solve(instance)
for tech in instance.installed_techs.values():
print(tech.value)
for number_of_techs_installed in instance.number_installed_techs.values():
print(number_of_techs_installed.value)
for TACperyear in instance.totalTACperyear.values():
print(TACperyear.value)
for CO2peryear in instance.totalCO2peryear.values():
print(CO2peryear.value)
print('Total TAC', instance.totalTAC())
print('Total CO2', instance.totalCO2())
If I have to calculate the "demandsperyear" outside of the Pyomo, how should i do that?
Thanks
To get something I had to solve, as @user2357112 suggested, the problem on total_TAC:
def total_TAC(m):
# return m.totalTAC == sum(m.totalTACperyear[year] for year in m.years)
return sum(m.totalTACperyear[year] for year in m.years)
Then I calculated demandsperyear before creating the instance:
def demandperyear_outside(demand, product, year, initial_demands, growths):
if year == 1:
return initial_demands[demand, product]
# I assume compound growth for the next years, maybe you want to change this
return initial_demands[demand, product] * (1 + growths[demand, product] * year)
demandsperyear_data = {
(demand, product, year): demandperyear_outside(
demand, product, year, data[None]['initial_demands'], data[None]['growths'])
for demand in data[None]['demands'][None]
for product in data[None]['products'][None]
for year in data[None]['years'][None]
}
which I add to the data:
data[None]['demandsperyear'] = demandsperyear_data
With this I change your m.demandsperyear to:
m.demandsperyear = Param(m.demands, m.products, m.years, initialize=data[None]['demandsperyear'])
I'm getting a Total CO2 (1400) but sadly a None for Total TAC.
Update for my second attempt: I added a constraint for totalTAC as follows:
def total_TAC_constraint(m):
return m.totalTAC == sum(m.totalTACperyear[year] for year in m.years)
m.totalTAC_cons = Constraint(rule=total_TAC_constraint)
and I'm now getting results:
Total TAC 7000.0
Total CO2 1400.0