1
|
|
|
# -*- coding: utf-8 -*- |
2
|
|
|
|
3
|
|
|
"""Solph Optimization Models. |
4
|
|
|
|
5
|
|
|
SPDX-FileCopyrightText: Uwe Krien <[email protected]> |
6
|
|
|
SPDX-FileCopyrightText: Simon Hilpert |
7
|
|
|
SPDX-FileCopyrightText: Cord Kaldemeyer |
8
|
|
|
SPDX-FileCopyrightText: gplssm |
9
|
|
|
SPDX-FileCopyrightText: Patrik Schönfeldt |
10
|
|
|
SPDX-FileCopyrightText: Saeed Sayadi |
11
|
|
|
SPDX-FileCopyrightText: Johannes Kochems |
12
|
|
|
SPDX-FileCopyrightText: Lennart Schürmann |
13
|
|
|
|
14
|
|
|
SPDX-License-Identifier: MIT |
15
|
|
|
|
16
|
|
|
""" |
17
|
|
|
import logging |
18
|
|
|
import warnings |
19
|
|
|
from logging import getLogger |
20
|
|
|
|
21
|
|
|
from oemof.tools import debugging |
22
|
|
|
from pyomo import environ as po |
23
|
|
|
from pyomo.core.plugins.transform.relax_integrality import RelaxIntegrality |
24
|
|
|
from pyomo.opt import SolverFactory |
25
|
|
|
|
26
|
|
|
from oemof.solph import processing |
27
|
|
|
from oemof.solph.buses._bus import BusBlock |
28
|
|
|
from oemof.solph.components._converter import ConverterBlock |
29
|
|
|
from oemof.solph.flows._invest_non_convex_flow_block import ( |
30
|
|
|
InvestNonConvexFlowBlock, |
31
|
|
|
) |
32
|
|
|
from oemof.solph.flows._investment_flow_block import InvestmentFlowBlock |
33
|
|
|
from oemof.solph.flows._non_convex_flow_block import NonConvexFlowBlock |
34
|
|
|
from oemof.solph.flows._simple_flow_block import SimpleFlowBlock |
35
|
|
|
|
36
|
|
|
|
37
|
|
|
class LoggingError(BaseException): |
38
|
|
|
"""Raised when the wrong logging level is used.""" |
39
|
|
|
|
40
|
|
|
pass |
41
|
|
|
|
42
|
|
|
|
43
|
|
|
class Model(po.ConcreteModel): |
44
|
|
|
"""An energy system model for operational and/or investment |
45
|
|
|
optimization. |
46
|
|
|
|
47
|
|
|
Parameters |
48
|
|
|
---------- |
49
|
|
|
energysystem : EnergySystem object |
50
|
|
|
Object that holds the nodes of an oemof energy system graph. |
51
|
|
|
constraint_groups : list |
52
|
|
|
Solph looks for these groups in the given energy system and uses them |
53
|
|
|
to create the constraints of the optimization problem. |
54
|
|
|
Defaults to `Model.CONSTRAINT_GROUPS` |
55
|
|
|
discount_rate : float or None |
56
|
|
|
The rate used for discounting in a multi-period model. |
57
|
|
|
A 2% discount rate needs to be defined as 0.02. |
58
|
|
|
objective_weighting : array like (optional) |
59
|
|
|
Weights used for temporal objective function |
60
|
|
|
expressions. If nothing is passed, `timeincrement` will be used which |
61
|
|
|
is calculated from the freq length of the energy system timeindex or |
62
|
|
|
can be directly passed as a sequence. |
63
|
|
|
auto_construct : boolean |
64
|
|
|
If this value is true, the set, variables, constraints, etc. are added, |
65
|
|
|
automatically when instantiating the model. For sequential model |
66
|
|
|
building process set this value to False |
67
|
|
|
and use methods `_add_parent_block_sets`, |
68
|
|
|
`_add_parent_block_variables`, `_add_blocks`, `_add_objective` |
69
|
|
|
|
70
|
|
|
Attributes |
71
|
|
|
---------- |
72
|
|
|
timeincrement : sequence |
73
|
|
|
Time increments |
74
|
|
|
flows : dict |
75
|
|
|
Flows of the model |
76
|
|
|
name : str |
77
|
|
|
Name of the model |
78
|
|
|
es : solph.EnergySystem |
79
|
|
|
Energy system of the model |
80
|
|
|
meta : `pyomo.opt.results.results_.SolverResults` or None |
81
|
|
|
Solver results |
82
|
|
|
dual : `pyomo.core.base.suffix.Suffix` or None |
83
|
|
|
Store the dual variables of the model if pyomo suffix is set to IMPORT |
84
|
|
|
rc : `pyomo.core.base.suffix.Suffix` or None |
85
|
|
|
Store the reduced costs of the model if pyomo suffix is set to IMPORT |
86
|
|
|
|
87
|
|
|
Note |
88
|
|
|
---- |
89
|
|
|
|
90
|
|
|
* The discount rate is only applicable for a multi-period model. |
91
|
|
|
* If you want to work with costs data in nominal terms, |
92
|
|
|
you should specify a discount rate. |
93
|
|
|
* By default, there is a discount rate of 2% in a multi-period model. |
94
|
|
|
* If you want to provide your costs data in real terms, |
95
|
|
|
just specify `discount_rate = 0`, i.e. effectively there will be |
96
|
|
|
no discounting. |
97
|
|
|
|
98
|
|
|
|
99
|
|
|
**The following basic sets are created**: |
100
|
|
|
|
101
|
|
|
NODES |
102
|
|
|
A set with all nodes of the given energy system. |
103
|
|
|
|
104
|
|
|
TIMESTEPS |
105
|
|
|
A set with all timesteps of the given time horizon. |
106
|
|
|
|
107
|
|
|
PERIODS |
108
|
|
|
A set with all investment periods of the given time horizon. |
109
|
|
|
|
110
|
|
|
TIMEINDEX |
111
|
|
|
A set with all time indices of the given time horizon, whereby |
112
|
|
|
time indices are defined as a tuple consisting of the period and the |
113
|
|
|
timestep. E.g. (2, 10) would be timestep 10 (which is exactly the same |
114
|
|
|
as in the TIMESTEPS set) and which is in period 2. |
115
|
|
|
|
116
|
|
|
FLOWS |
117
|
|
|
A 2 dimensional set with all flows. Index: `(source, target)` |
118
|
|
|
|
119
|
|
|
**The following basic variables are created**: |
120
|
|
|
|
121
|
|
|
flow |
122
|
|
|
Flow from source to target indexed by FLOWS, TIMEINDEX. |
123
|
|
|
Note: Bounds of this variable are set depending on attributes of |
124
|
|
|
the corresponding flow object. |
125
|
|
|
|
126
|
|
|
""" |
127
|
|
|
|
128
|
|
|
CONSTRAINT_GROUPS = [ |
129
|
|
|
BusBlock, |
130
|
|
|
ConverterBlock, |
131
|
|
|
InvestmentFlowBlock, |
132
|
|
|
SimpleFlowBlock, |
133
|
|
|
NonConvexFlowBlock, |
134
|
|
|
InvestNonConvexFlowBlock, |
135
|
|
|
] |
136
|
|
|
|
137
|
|
|
def __init__(self, energysystem, discount_rate=None, **kwargs): |
138
|
|
|
super().__init__() |
139
|
|
|
|
140
|
|
|
# Check root logger. Due to a problem with pyomo the building of the |
141
|
|
|
# model will take up to a 100 times longer if the root logger is set |
142
|
|
|
# to DEBUG |
143
|
|
|
|
144
|
|
|
if getLogger().level <= 10 and kwargs.get("debug", False) is False: |
145
|
|
|
msg = ( |
146
|
|
|
"The root logger level is 'DEBUG'.\nDue to a communication " |
147
|
|
|
"problem between solph and the pyomo package,\nusing the " |
148
|
|
|
"DEBUG level will slow down the modelling process by the " |
149
|
|
|
"factor ~100.\nIf you need the debug-logging you can " |
150
|
|
|
"initialise the Model with 'debug=True`\nYou should only do " |
151
|
|
|
"this for small models. To avoid the slow-down use the " |
152
|
|
|
"logger\nfunction of oemof.tools (read docstring) or " |
153
|
|
|
"change the level of the root logger:\n\nimport logging\n" |
154
|
|
|
"logging.getLogger().setLevel(logging.INFO)" |
155
|
|
|
) |
156
|
|
|
raise LoggingError(msg) |
157
|
|
|
|
158
|
|
|
# ######################## Arguments ################################# |
159
|
|
|
|
160
|
|
|
self.name = kwargs.get("name", type(self).__name__) |
161
|
|
|
|
162
|
|
|
self.es = energysystem |
163
|
|
|
|
164
|
|
|
if kwargs.get("timeincrement"): |
165
|
|
|
msg = "Resetting timeincrement from EnergySystem in Model." |
166
|
|
|
warnings.warn(msg, debugging.SuspiciousUsageWarning) |
167
|
|
|
|
168
|
|
|
self.timeincrement = kwargs.get("timeincrement") |
169
|
|
|
else: |
170
|
|
|
self.timeincrement = self.es.timeincrement |
171
|
|
|
|
172
|
|
|
self.objective_weighting = kwargs.get( |
173
|
|
|
"objective_weighting", self.timeincrement |
174
|
|
|
) |
175
|
|
|
|
176
|
|
|
self._constraint_groups = type(self).CONSTRAINT_GROUPS + kwargs.get( |
177
|
|
|
"constraint_groups", [] |
178
|
|
|
) |
179
|
|
|
|
180
|
|
|
self._constraint_groups += [ |
181
|
|
|
i |
182
|
|
|
for i in self.es.groups |
183
|
|
|
if hasattr(i, "CONSTRAINT_GROUP") |
184
|
|
|
and i not in self._constraint_groups |
185
|
|
|
] |
186
|
|
|
|
187
|
|
|
self.flows = self.es.flows() |
188
|
|
|
|
189
|
|
|
self.solver_results = None |
190
|
|
|
self.dual = None |
191
|
|
|
self.rc = None |
192
|
|
|
|
193
|
|
|
if discount_rate is not None: |
194
|
|
|
self.discount_rate = discount_rate |
195
|
|
|
elif energysystem.periods is not None: |
196
|
|
|
self._set_discount_rate_with_warning() |
197
|
|
|
else: |
198
|
|
|
pass |
199
|
|
|
|
200
|
|
|
if kwargs.get("auto_construct", True): |
201
|
|
|
self._construct() |
202
|
|
|
|
203
|
|
|
def _construct(self): |
204
|
|
|
"""Construct a Model by adding parent block sets and variables |
205
|
|
|
as well as child blocks and variables to it. |
206
|
|
|
""" |
207
|
|
|
self._add_parent_block_sets() |
208
|
|
|
self._add_parent_block_variables() |
209
|
|
|
self._add_child_blocks() |
210
|
|
|
self._add_objective() |
211
|
|
|
|
212
|
|
|
def _set_discount_rate_with_warning(self): |
213
|
|
|
""" |
214
|
|
|
Sets the discount rate to the standard value and raises a warning. |
215
|
|
|
""" |
216
|
|
|
self.discount_rate = 0.02 |
217
|
|
|
msg = ( |
218
|
|
|
f"By default, a discount_rate of {self.discount_rate} " |
219
|
|
|
f"is used for a multi-period model. " |
220
|
|
|
f"If you want to use another value, " |
221
|
|
|
f"you have to specify the `discount_rate` attribute." |
222
|
|
|
) |
223
|
|
|
warnings.warn(msg, debugging.SuspiciousUsageWarning) |
224
|
|
|
|
225
|
|
|
def _add_parent_block_sets(self): |
226
|
|
|
"""Add all basic sets to the model, i.e. NODES, TIMESTEPS and FLOWS. |
227
|
|
|
Also create sets PERIODS and TIMEINDEX used for multi-period models. |
228
|
|
|
""" |
229
|
|
|
self.nodes = list(self.es.nodes) |
230
|
|
|
|
231
|
|
|
# create set with all nodes |
232
|
|
|
self.NODES = po.Set(initialize=[n for n in self.nodes]) |
233
|
|
|
|
234
|
|
|
if self.es.timeincrement is None: |
235
|
|
|
msg = ( |
236
|
|
|
"The EnergySystem needs to have a valid 'timeincrement' " |
237
|
|
|
"attribute to build a model." |
238
|
|
|
) |
239
|
|
|
raise AttributeError(msg) |
240
|
|
|
|
241
|
|
|
# pyomo set for timesteps of optimization problem |
242
|
|
|
self.TIMESTEPS = po.Set( |
243
|
|
|
initialize=range(len(self.es.timeincrement)), ordered=True |
244
|
|
|
) |
245
|
|
|
self.TIMEPOINTS = po.Set( |
246
|
|
|
initialize=range(len(self.es.timeincrement) + 1), ordered=True |
247
|
|
|
) |
248
|
|
|
|
249
|
|
|
if self.es.periods is None: |
250
|
|
|
self.TIMEINDEX = po.Set( |
251
|
|
|
initialize=list( |
252
|
|
|
zip( |
253
|
|
|
[0] * len(self.es.timeincrement), |
254
|
|
|
range(len(self.es.timeincrement)), |
255
|
|
|
) |
256
|
|
|
), |
257
|
|
|
ordered=True, |
258
|
|
|
) |
259
|
|
|
self.PERIODS = po.Set(initialize=[0]) |
260
|
|
|
else: |
261
|
|
|
nested_list = [ |
262
|
|
|
[k] * len(self.es.periods[k]) |
263
|
|
|
for k in range(len(self.es.periods)) |
264
|
|
|
] |
265
|
|
|
flattened_list = [ |
266
|
|
|
item for sublist in nested_list for item in sublist |
267
|
|
|
] |
268
|
|
|
self.TIMEINDEX = po.Set( |
269
|
|
|
initialize=list( |
270
|
|
|
zip(flattened_list, range(len(self.es.timeincrement))) |
271
|
|
|
), |
272
|
|
|
ordered=True, |
273
|
|
|
) |
274
|
|
|
self.PERIODS = po.Set( |
275
|
|
|
initialize=sorted(list(set(range(len(self.es.periods))))) |
276
|
|
|
) |
277
|
|
|
|
278
|
|
|
# (Re-)Map timesteps to periods |
279
|
|
|
timesteps_in_period = {p: [] for p in self.PERIODS} |
280
|
|
|
for p, t in self.TIMEINDEX: |
281
|
|
|
timesteps_in_period[p].append(t) |
282
|
|
|
self.TIMESTEPS_IN_PERIOD = timesteps_in_period |
283
|
|
|
|
284
|
|
|
# previous timesteps |
285
|
|
|
previous_timesteps = [x - 1 for x in self.TIMESTEPS] |
286
|
|
|
previous_timesteps[0] = self.TIMESTEPS.last() |
287
|
|
|
|
288
|
|
|
self.previous_timesteps = dict(zip(self.TIMESTEPS, previous_timesteps)) |
289
|
|
|
|
290
|
|
|
# pyomo set for all flows in the energy system graph |
291
|
|
|
self.FLOWS = po.Set( |
292
|
|
|
initialize=self.flows.keys(), ordered=True, dimen=2 |
293
|
|
|
) |
294
|
|
|
|
295
|
|
|
self.BIDIRECTIONAL_FLOWS = po.Set( |
296
|
|
|
initialize=[k for (k, v) in self.flows.items() if v.bidirectional], |
297
|
|
|
ordered=True, |
298
|
|
|
dimen=2, |
299
|
|
|
within=self.FLOWS, |
300
|
|
|
) |
301
|
|
|
|
302
|
|
|
self.UNIDIRECTIONAL_FLOWS = po.Set( |
303
|
|
|
initialize=[ |
304
|
|
|
k for (k, v) in self.flows.items() if not v.bidirectional |
305
|
|
|
], |
306
|
|
|
ordered=True, |
307
|
|
|
dimen=2, |
308
|
|
|
within=self.FLOWS, |
309
|
|
|
) |
310
|
|
|
|
311
|
|
|
def _add_parent_block_variables(self): |
312
|
|
|
"""Add the parent block variables, which is the `flow` variable, |
313
|
|
|
indexed by FLOWS and TIMEINDEX.""" |
314
|
|
|
self.flow = po.Var(self.FLOWS, self.TIMESTEPS, within=po.Reals) |
315
|
|
|
|
316
|
|
|
for o, i in self.FLOWS: |
317
|
|
|
if self.flows[o, i].nominal_value is not None: |
318
|
|
|
if self.flows[o, i].fix[self.TIMESTEPS.at(1)] is not None: |
319
|
|
|
for t in self.TIMESTEPS: |
320
|
|
|
self.flow[o, i, t].value = ( |
321
|
|
|
self.flows[o, i].fix[t] |
322
|
|
|
* self.flows[o, i].nominal_value |
323
|
|
|
) |
324
|
|
|
self.flow[o, i, t].fix() |
325
|
|
|
else: |
326
|
|
|
for t in self.TIMESTEPS: |
327
|
|
|
self.flow[o, i, t].setub( |
328
|
|
|
self.flows[o, i].max[t] |
329
|
|
|
* self.flows[o, i].nominal_value |
330
|
|
|
) |
331
|
|
|
if not self.flows[o, i].nonconvex: |
332
|
|
|
for t in self.TIMESTEPS: |
333
|
|
|
self.flow[o, i, t].setlb( |
334
|
|
|
self.flows[o, i].min[t] |
335
|
|
|
* self.flows[o, i].nominal_value |
336
|
|
|
) |
337
|
|
|
elif (o, i) in self.UNIDIRECTIONAL_FLOWS: |
338
|
|
|
for t in self.TIMESTEPS: |
339
|
|
|
self.flow[o, i, t].setlb(0) |
340
|
|
|
else: |
341
|
|
|
if (o, i) in self.UNIDIRECTIONAL_FLOWS: |
342
|
|
|
for t in self.TIMESTEPS: |
343
|
|
|
self.flow[o, i, t].setlb(0) |
344
|
|
|
|
345
|
|
|
def _add_child_blocks(self): |
346
|
|
|
"""Method to add the defined child blocks for components that have |
347
|
|
|
been grouped in the defined constraint groups. This collects all the |
348
|
|
|
constraints from the buses, components and flows blocks |
349
|
|
|
and adds them to the model. |
350
|
|
|
""" |
351
|
|
|
for group in self._constraint_groups: |
352
|
|
|
block = group() |
353
|
|
|
self.add_component(str(block), block) |
354
|
|
|
|
355
|
|
|
# create constraints etc. related with block for all nodes |
356
|
|
|
# in the group |
357
|
|
|
block._create(group=self.es.groups.get(group)) |
358
|
|
|
|
359
|
|
|
def _add_objective(self, sense=po.minimize, update=False): |
360
|
|
|
"""Method to sum up all objective expressions from the child blocks |
361
|
|
|
that have been created. This method looks for `_objective_expression` |
362
|
|
|
attribute in the block definition and will call this method to add |
363
|
|
|
their return value to the objective function. |
364
|
|
|
""" |
365
|
|
|
if update: |
366
|
|
|
self.del_component("objective") |
367
|
|
|
|
368
|
|
|
expr = 0 |
369
|
|
|
|
370
|
|
|
for block in self.component_data_objects(): |
371
|
|
|
if hasattr(block, "_objective_expression"): |
372
|
|
|
expr += block._objective_expression() |
373
|
|
|
|
374
|
|
|
self.objective = po.Objective(sense=sense, expr=expr) |
375
|
|
|
|
376
|
|
|
def receive_duals(self): |
377
|
|
|
"""Method sets solver suffix to extract information about dual |
378
|
|
|
variables from solver. Shadow prices (duals) and reduced costs (rc) are |
379
|
|
|
set as attributes of the model. |
380
|
|
|
""" |
381
|
|
|
# shadow prices |
382
|
|
|
self.dual = po.Suffix(direction=po.Suffix.IMPORT) |
383
|
|
|
# reduced costs |
384
|
|
|
self.rc = po.Suffix(direction=po.Suffix.IMPORT) |
385
|
|
|
|
386
|
|
|
def results(self): |
387
|
|
|
"""Returns a nested dictionary of the results of this optimization. |
388
|
|
|
See the processing module for more information on results extraction. |
389
|
|
|
""" |
390
|
|
|
return processing.results(self) |
391
|
|
|
|
392
|
|
|
def solve(self, solver="cbc", solver_io="lp", **kwargs): |
393
|
|
|
r"""Takes care of communication with solver to solve the model. |
394
|
|
|
|
395
|
|
|
Parameters |
396
|
|
|
---------- |
397
|
|
|
solver : string |
398
|
|
|
solver to be used e.g. "cbc", "glpk", "gurobi", "cplex" |
399
|
|
|
solver_io : string |
400
|
|
|
pyomo solver interface file format: "lp", "python", "nl", etc. |
401
|
|
|
\**kwargs : keyword arguments |
402
|
|
|
Possible keys can be set see below: |
403
|
|
|
|
404
|
|
|
Other Parameters |
405
|
|
|
---------------- |
406
|
|
|
solve_kwargs : dict |
407
|
|
|
Other arguments for the pyomo.opt.SolverFactory.solve() method |
408
|
|
|
Example : {"tee":True} |
409
|
|
|
cmdline_options : dict |
410
|
|
|
Dictionary with command line options for solver e.g. |
411
|
|
|
{"mipgap":"0.01"} results in "--mipgap 0.01" |
412
|
|
|
\{"interior":" "} results in "--interior" |
413
|
|
|
\Gurobi solver takes numeric parameter values such as |
414
|
|
|
{"method": 2} |
415
|
|
|
""" |
416
|
|
|
solve_kwargs = kwargs.get("solve_kwargs", {}) |
417
|
|
|
solver_cmdline_options = kwargs.get("cmdline_options", {}) |
418
|
|
|
|
419
|
|
|
opt = SolverFactory(solver, solver_io=solver_io) |
420
|
|
|
# set command line options |
421
|
|
|
options = opt.options |
422
|
|
|
for k in solver_cmdline_options: |
423
|
|
|
options[k] = solver_cmdline_options[k] |
424
|
|
|
|
425
|
|
|
solver_results = opt.solve(self, **solve_kwargs) |
426
|
|
|
|
427
|
|
|
status = solver_results["Solver"][0]["Status"] |
428
|
|
|
termination_condition = solver_results["Solver"][0][ |
429
|
|
|
"Termination condition" |
430
|
|
|
] |
431
|
|
|
|
432
|
|
|
if status == "ok" and termination_condition == "optimal": |
433
|
|
|
logging.info("Optimization successful...") |
434
|
|
|
else: |
435
|
|
|
msg = ( |
436
|
|
|
"Optimization ended with status {0} and termination " |
437
|
|
|
"condition {1}" |
438
|
|
|
) |
439
|
|
|
warnings.warn( |
440
|
|
|
msg.format(status, termination_condition), UserWarning |
441
|
|
|
) |
442
|
|
|
self.es.results = solver_results |
443
|
|
|
self.solver_results = solver_results |
444
|
|
|
|
445
|
|
|
return solver_results |
446
|
|
|
|
447
|
|
|
def relax_problem(self): |
448
|
|
|
"""Relaxes integer variables to reals of optimization model self.""" |
449
|
|
|
relaxer = RelaxIntegrality() |
450
|
|
|
relaxer._apply_to(self) |
451
|
|
|
|
452
|
|
|
return self |
453
|
|
|
|