我正在尝试创建一个模型来最小化电解槽的 LCOH。价格有每小时的价值，应根据生产的数量来决定是否值得生产。所有其他值最初都是恒定的。这个想法是，投资成本和可变成本除以生产总量，从而使 LCOH 下降到一定的边际价格。 LCOH 的计算方法是投资成本加上可变成本除以产量。投资成本一开始就固定，剩下的由生产时间决定。 如果我尝试将固定成本包含在分子中并使用决策变量的总和来表示可变成本，我会得到一个错误，即我的函数不是线性的。我尝试了两种不同的方法来定义我的目标函数，但都不起作用。有人可以帮我解决这个问题吗？
import pandas as pd
from pyomo.environ import *
# prices with values in €/MWh
el_prices =[-5.17, -1.07, -1.47, -5.08, -4.49, -5.4, -5.02, -1.3, -1.44, -1.09, -1.07, -1.07, -0.79, -0.27, 0.85, 23.53, 36.54, 46.03, 55.57, 54.95, 49.23, 44.99, 45.96, 35.0, 57.91, 51.67, 52.93, 44.09, 50.08, 69.72, 105.08, 140.64, 145.98, 147.05, 145.61, 143.35, 144.38, 143.76, 148.2, 155.34, 162.89, 170.0, 174.74, 164.46, 153.0, 141.67, 134.91, 124.22, 130.01, 120.0]
investment_cost_ely = 1500 # Beispiel-Fixkosten in Euro/kw
power_ely = 1000 # Beispiel-Leistung des Elektrolyseurs in kW
conversion_factor = 0.6 # Beispiel-Wirkungsgrad des Elektrolyseurs
number_years = 20
interest_rate = 0.06
a = (((1+interest_rate)**number_years)*interest_rate)/(((1+interest_rate)**number_years)-1)
annualized_cost = investment_cost_ely * power_ely * a
m3 = ConcreteModel()
m3.betriebsstunden = Var(hours, domain=Binary)
def gestehungskosten_rule(model):
return fixe_kosten + sum(m3.betriebsstunden[t] * power_ely * el_prices[t]) / sum(m3.betriebsstunden[t] * power_ely * conversion_factor for t in hours)
m3.gestehungskosten = Objective(rule=gestehungskosten_rule, sense=minimize)
我收到此错误
ERROR: Rule failed when generating expression for Objective gestehungskosten
with index None: NameError: name 't' is not defined
ERROR: Constructing component 'gestehungskosten' from data=None failed:
NameError: name 't' is not defined
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 10
8 def gestehungskosten_rule(model):
9 return fixe_kosten + sum(m3.betriebsstunden[t] * leistung * strompreise[t]) / sum(m3.betriebsstunden[t] * leistung * wirkungsgrad for t in hours)
---> 10 m3.gestehungskosten = Objective(rule=gestehungskosten_rule, sense=minimize)
12 m3.gestehungskosten.pprint()
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/core/base/block.py:568, in _BlockData.__setattr__(self, name, val)
563 if name not in self.__dict__:
564 if isinstance(val, Component):
565 #
566 # Pyomo components are added with the add_component method.
567 #
--> 568 self.add_component(name, val)
569 else:
570 #
571 # Other Python objects are added with the standard __setattr__
572 # method.
573 #
574 super(_BlockData, self).__setattr__(name, val)
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/core/base/block.py:1126, in _BlockData.add_component(self, name, val)
1118 logger.debug(
1119 "Constructing %s '%s' on %s from data=%s",
1120 val.__class__.__name__,
(...)
1123 str(data),
1124 )
1125 try:
-> 1126 val.construct(data)
1127 except:
1128 err = sys.exc_info()[1]
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/core/base/objective.py:318, in Objective.construct(self, data)
315 else:
316 # Bypass the index validation and create the member directly
317 for index in self.index_set():
--> 318 ans = self._setitem_when_not_present(index, rule(block, index))
319 if ans is not None:
320 ans.set_sense(self._init_sense(block, index))
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/core/base/initializer.py:438, in ScalarCallInitializer.__call__(self, parent, idx)
437 def __call__(self, parent, idx):
--> 438 return self._fcn(parent)
Cell In[7], line 9, in gestehungskosten_rule(model)
8 def gestehungskosten_rule(model):
----> 9 return fixe_kosten + sum(m3.betriebsstunden[t] * leistung * strompreise[t]) / sum(m3.betriebsstunden[t] * leistung * wirkungsgrad for t in hours)
NameError: name 't' is not defined
这是我的第二种方法
#define promo model second ry
m6 = ConcreteModel()
# variables
fix_costs = 100000
m6.betriebsstunden = Var(hours, domain=Binary)
def gestehungskosten_rule(model):
variable_costs = sum(m6.betriebsstunden[t] * power_ely * el_prices[t] for t in hours)
quantity = sum(m6.betriebsstunden[t] * power_ely * conversion_factor for t in hours)
return (fix_costs + variable_costs) / quantity
m6.gestehungskosten = Objective(rule=gestehungskosten_rule, sense=minimize)
# solve problem
solver = SolverFactory('glpk') # Anpassen des Solvers je nach Bedarf (glpk)
results = solver.solve(m6)
# print result
print("Minimale Gestehungskosten:", m6.gestehungskosten())
print("Anzahl der Betriebsstunden:", sum(m6.betriebsstunden[t].value for t in hours))
print("Stunden mit Produktion:")
for t in hours:
if m6.betriebsstunden[t].value == 1:
print("Stunde:", t, "Strompreis:", el_prices[t])
这里出现这个错误
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[6], line 3
1 # Solver konfigurieren und Problem lösen
2 solver = SolverFactory('glpk') # Anpassen des Solvers je nach Bedarf (glpk)
----> 3 results = solver.solve(m6)
5 # Ergebnisse ausgeben
6 print("Minimale Gestehungskosten:", m6.gestehungskosten())
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/opt/base/solvers.py:598, in OptSolver.solve(self, *args, **kwds)
594 try:
595 # we're good to go.
596 initial_time = time.time()
--> 598 self._presolve(*args, **kwds)
600 presolve_completion_time = time.time()
601 if self._report_timing:
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/opt/solver/shellcmd.py:222, in SystemCallSolver._presolve(self, *args, **kwds)
219 self._keepfiles = kwds.pop("keepfiles", False)
220 self._define_signal_handlers = kwds.pop('use_signal_handling', None)
--> 222 OptSolver._presolve(self, *args, **kwds)
224 #
225 # Verify that the input problems exists
226 #
227 for filename in self._problem_files:
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/opt/base/solvers.py:704, in OptSolver._presolve(self, *args, **kwds)
701 if self._problem_format:
702 write_start_time = time.time()
703 (self._problem_files, self._problem_format, self._smap_id) = (
--> 704 self._convert_problem(
705 args, self._problem_format, self._valid_problem_formats, **kwds
706 )
707 )
708 total_time = time.time() - write_start_time
709 if self._report_timing:
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/opt/base/solvers.py:756, in OptSolver._convert_problem(self, args, problem_format, valid_problem_formats, **kwds)
755 def _convert_problem(self, args, problem_format, valid_problem_formats, **kwds):
--> 756 return convert_problem(
757 args, problem_format, valid_problem_formats, self.has_capability, **kwds
758 )
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/opt/base/convert.py:97, in convert_problem(args, target_problem_type, valid_problem_types, has_capability, **kwds)
95 tmpkw = kwds
96 tmpkw['capabilities'] = has_capability
---> 97 problem_files, symbol_map = converter.apply(*tmp, **tmpkw)
98 return problem_files, ptype, symbol_map
100 msg = 'No conversion possible. Source problem type: %s. Valid target types: %s'
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/solvers/plugins/converter/model.py:78, in PyomoMIPConverter.apply(self, *args, **kwds)
70 symbol_map_id = instance.write(
71 problem_filename,
72 format=ProblemFormat.cpxlp,
(...)
75 **io_options
76 )
77 else:
---> 78 (problem_filename, symbol_map_id) = instance.write(
79 filename=problem_filename,
80 format=ProblemFormat.cpxlp,
81 solver_capability=capabilities,
82 io_options=io_options,
83 )
84 return (problem_filename,), symbol_map_id
85 else:
86 #
87 # I'm simply exposing a fatal issue with
(...)
90 # arguments that can be sent to the writer?
91 #
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/core/base/block.py:1937, in _BlockData.write(self, filename, format, solver_capability, io_options, int_marker)
1934 def solver_capability(x):
1935 return True
-> 1937 (filename, smap) = problem_writer(self, filename, solver_capability, io_options)
1938 smap_id = id(smap)
1939 if not hasattr(self, 'solutions'):
1940 # This is a bit of a hack. The write() method was moved
1941 # here from PyomoModel to support the solution of arbitrary
(...)
1946 # dependency (we only need it here because we store the
1947 # SymbolMap returned by the writer in the solutions).
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/repn/plugins/lp_writer.py:208, in LPWriter.__call__(self, model, filename, solver_capability, io_options)
205 io_options['allow_quadratic_constraint'] = qc
207 with open(filename, 'w', newline='') as FILE:
--> 208 info = self.write(model, FILE, **io_options)
209 return filename, info.symbol_map
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/repn/plugins/lp_writer.py:241, in LPWriter.write(self, model, ostream, **options)
237 # Pause the GC, as the walker that generates the compiled LP
238 # representation generates (and disposes of) a large number of
239 # small objects.
240 with PauseGC():
--> 241 return _LPWriter_impl(ostream, config).write(model)
File /opt/anaconda3/lib/python3.11/site-packages/pyomo/repn/plugins/lp_writer.py:378, in _LPWriter_impl.write(self, model)
376 repn = objective_visitor.walk_expression(obj.expr)
377 if repn.nonlinear is not None:
--> 378 raise ValueError(
379 f"Model objective ({obj.name}) contains nonlinear terms that "
380 "cannot be written to LP format"
381 )
382 if repn.constant or not (repn.linear or getattr(repn, 'quadratic', None)):
383 # Older versions of CPLEX (including 12.6) and all versions
384 # of GLPK (through 5.0) do not support constants in the
(...)
390 # objective, this will ensure we at least write out
391 # 0*ONE_VAR_CONSTANT.
392 repn.linear[id(ONE_VAR_CONSTANT)] = repn.constant
ValueError: Model objective (gestehungskosten) contains nonlinear terms that cannot be written to LP format
感谢您的帮助
如果你面临非线性约束，你可以尝试的是使用另一个Python包：docplex，并且该API允许非线性约束和目标，因为这依赖于约束编程
小示例位于 https://github.com/AlexFleischerParis/zoodocplex/blob/master/zoononlinear.py
from docplex.cp.model import CpoModel
mdl = CpoModel(name='buses')
nbbus40 = mdl.integer_var(0,1000,name='nbBus40')
nbbus30 = mdl.integer_var(0,1000,name='nbBus30')
mdl.add(nbbus40*40 + nbbus30*30 >= 300)
#non linear objective
mdl.minimize(mdl.exponent(nbbus40)*500 + mdl.exponent(nbbus30)*400)
msol=mdl.solve()
print(msol[nbbus40]," buses 40 seats")
print(msol[nbbus30]," buses 30 seats")