solph._models.Model.solve()   B
last analyzed

Complexity

Conditions 5

Size

Total Lines 64
Code Lines 22

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
eloc 22
dl 0
loc 64
rs 8.8853
c 0
b 0
f 0
cc 5
nop 5

How to fix   Long Method   

Long Method

Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.

For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.

Commonly applied refactorings include:

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
    objective_weighting : array like (optional)
56
        Weights used for temporal objective function
57
        expressions. If nothing is passed, `timeincrement` will be used which
58
        is calculated from the freq length of the energy system timeindex or
59
        can be directly passed as a sequence.
60
    auto_construct : boolean
61
        If this value is true, the set, variables, constraints, etc. are added,
62
        automatically when instantiating the model. For sequential model
63
        building process set this value to False
64
        and use methods `_add_parent_block_sets`,
65
        `_add_parent_block_variables`, `_add_blocks`, `_add_objective`
66
67
    Attributes
68
    ----------
69
    timeincrement : sequence
70
        Time increments
71
    flows : dict
72
        Flows of the model
73
    name : str
74
        Name of the model
75
    es : solph.EnergySystem
76
        Energy system of the model
77
    meta : `pyomo.opt.results.results_.SolverResults` or None
78
        Solver results
79
    dual : `pyomo.core.base.suffix.Suffix` or None
80
        Store the dual variables of the model if pyomo suffix is set to IMPORT
81
    rc : `pyomo.core.base.suffix.Suffix` or None
82
        Store the reduced costs of the model if pyomo suffix is set to IMPORT
83
84
85
    **The following basic sets are created**:
86
87
    NODES
88
        A set with all nodes of the given energy system.
89
90
    TIMESTEPS
91
        A set with all timesteps of the given time horizon.
92
93
    PERIODS
94
        A set with all investment periods of the given time horizon.
95
96
    TIMEINDEX
97
        A set with all time indices of the given time horizon, whereby
98
        time indices are defined as a tuple consisting of the period and the
99
        timestep. E.g. (2, 10) would be timestep 10 (which is exactly the same
100
        as in the TIMESTEPS set) and which is in period 2.
101
102
    FLOWS
103
        A 2 dimensional set with all flows. Index: `(source, target)`
104
105
    **The following basic variables are created**:
106
107
    flow
108
        Flow from source to target indexed by FLOWS, TIMEINDEX.
109
        Note: Bounds of this variable are set depending on attributes of
110
        the corresponding flow object.
111
112
    """
113
114
    CONSTRAINT_GROUPS = [
115
        BusBlock,
116
        ConverterBlock,
117
        InvestmentFlowBlock,
118
        SimpleFlowBlock,
119
        NonConvexFlowBlock,
120
        InvestNonConvexFlowBlock,
121
    ]
122
123
    def __init__(self, energysystem, **kwargs):
124
        super().__init__()
125
126
        # Check root logger. Due to a problem with pyomo the building of the
127
        # model will take up to a 100 times longer if the root logger is set
128
        # to DEBUG
129
130
        if getLogger().level <= 10 and kwargs.get("debug", False) is False:
131
            msg = (
132
                "The root logger level is 'DEBUG'.\nDue to a communication "
133
                "problem between solph and the pyomo package,\nusing the "
134
                "DEBUG level will slow down the modelling process by the "
135
                "factor ~100.\nIf you need the debug-logging you can "
136
                "initialise the Model with 'debug=True`\nYou should only do "
137
                "this for small models. To avoid the slow-down use the "
138
                "logger\nfunction of oemof.tools (read docstring) or "
139
                "change the level of the root logger:\n\nimport logging\n"
140
                "logging.getLogger().setLevel(logging.INFO)"
141
            )
142
            raise LoggingError(msg)
143
144
        # ########################  Arguments #################################
145
146
        self.name = kwargs.get("name", type(self).__name__)
147
148
        self.es = energysystem
149
150
        if kwargs.get("timeincrement"):
151
            msg = "Resetting timeincrement from EnergySystem in Model."
152
            warnings.warn(msg, debugging.SuspiciousUsageWarning)
153
154
            self.timeincrement = kwargs.get("timeincrement")
155
        else:
156
            self.timeincrement = self.es.timeincrement
157
158
        self.objective_weighting = kwargs.get(
159
            "objective_weighting", self.timeincrement
160
        )
161
162
        self._constraint_groups = type(self).CONSTRAINT_GROUPS + kwargs.get(
163
            "constraint_groups", []
164
        )
165
166
        self._constraint_groups += [
167
            i
168
            for i in self.es.groups
169
            if hasattr(i, "CONSTRAINT_GROUP")
170
            and i not in self._constraint_groups
171
        ]
172
173
        self.flows = self.es.flows()
174
175
        self.solver_results = None
176
        self.dual = None
177
        self.rc = None
178
179
        if energysystem.periods is not None:
180
            self._set_discount_rate_with_warning()
181
        else:
182
            pass
183
184
        if kwargs.get("auto_construct", True):
185
            self._construct()
186
187
    def _construct(self):
188
        """Construct a Model by adding parent block sets and variables
189
        as well as child blocks and variables to it.
190
        """
191
        self._add_parent_block_sets()
192
        self._add_parent_block_variables()
193
        self._add_child_blocks()
194
        self._add_objective()
195
196
    def _set_discount_rate_with_warning(self):
197
        """
198
        Sets the discount rate to the standard value and raises a warning.
199
        """
200
        self.discount_rate = 0.02
201
        msg = (
202
            f"By default, a discount_rate of {self.discount_rate} "
203
            f"is used for a multi-period model. "
204
            f"If you want to use another value, "
205
            f"you have to specify the `discount_rate` attribute."
206
        )
207
        warnings.warn(msg, debugging.SuspiciousUsageWarning)
208
209
    def _add_parent_block_sets(self):
210
        """Add all basic sets to the model, i.e. NODES, TIMESTEPS and FLOWS.
211
        Also create sets PERIODS and TIMEINDEX used for multi-period models.
212
        """
213
        self.nodes = list(self.es.nodes)
214
215
        # create set with all nodes
216
        self.NODES = po.Set(initialize=[n for n in self.nodes])
217
218
        if self.es.timeincrement is None:
219
            msg = (
220
                "The EnergySystem needs to have a valid 'timeincrement' "
221
                "attribute to build a model."
222
            )
223
            raise AttributeError(msg)
224
225
        # pyomo set for timesteps of optimization problem
226
        self.TIMESTEPS = po.Set(
227
            initialize=range(len(self.es.timeincrement)), ordered=True
228
        )
229
        self.TIMEPOINTS = po.Set(
230
            initialize=range(len(self.es.timeincrement) + 1), ordered=True
231
        )
232
233
        if self.es.periods is None:
234
            self.TIMEINDEX = po.Set(
235
                initialize=list(
236
                    zip(
237
                        [0] * len(self.es.timeincrement),
238
                        range(len(self.es.timeincrement)),
239
                    )
240
                ),
241
                ordered=True,
242
            )
243
            self.PERIODS = po.Set(initialize=[0])
244
        else:
245
            nested_list = [
246
                [k] * len(self.es.periods[k])
247
                for k in range(len(self.es.periods))
248
            ]
249
            flattened_list = [
250
                item for sublist in nested_list for item in sublist
251
            ]
252
            self.TIMEINDEX = po.Set(
253
                initialize=list(
254
                    zip(flattened_list, range(len(self.es.timeincrement)))
255
                ),
256
                ordered=True,
257
            )
258
            self.PERIODS = po.Set(
259
                initialize=sorted(list(set(range(len(self.es.periods)))))
260
            )
261
262
        # (Re-)Map timesteps to periods
263
        timesteps_in_period = {p: [] for p in self.PERIODS}
264
        for p, t in self.TIMEINDEX:
265
            timesteps_in_period[p].append(t)
266
        self.TIMESTEPS_IN_PERIOD = timesteps_in_period
267
268
        # Set up disaggregated timesteps from original timeseries
269
        self.TSAM_MODE = False
270
        if self.es.tsa_parameters is None:
271
            self.tsam_weighting = [1] * len(self.timeincrement)
272
        else:
273
            self.TSAM_MODE = True
274
275
            # Construct weighting from occurrences and order
276
            self.tsam_weighting = list(
277
                self.es.tsa_parameters[p]["occurrences"][k]
278
                for p in self.PERIODS
279
                for k in range(len(self.es.tsa_parameters[p]["occurrences"]))
280
                for _ in range(self.es.tsa_parameters[p]["timesteps"])
281
            )
282
            self.CLUSTERS = po.Set(
283
                initialize=list(
284
                    range(
285
                        sum(
286
                            len(self.es.tsa_parameters[p]["order"])
287
                            for p in self.PERIODS
288
                        )
289
                    )
290
                )
291
            )
292
            self.CLUSTERS_OFFSET = po.Set(
293
                initialize=list(
294
                    range(
295
                        sum(
296
                            len(self.es.tsa_parameters[p]["order"])
297
                            for p in self.PERIODS
298
                        )
299
                        + 1
300
                    )
301
                )
302
            )
303
            self.TYPICAL_CLUSTERS = po.Set(
304
                initialize=[
305
                    (p, i)
306
                    for p in self.PERIODS
307
                    for i in range(
308
                        len(self.es.tsa_parameters[p]["occurrences"])
309
                    )
310
                ]
311
            )
312
313
            self.TIMEINDEX_CLUSTER = self.get_cluster_index("order", 0)
314
            self.TIMEINDEX_TYPICAL_CLUSTER = self.get_cluster_index(
315
                "occurrences", 0
316
            )
317
            self.TIMEINDEX_TYPICAL_CLUSTER_OFFSET = self.get_cluster_index(
318
                "occurrences", 1
319
            )
320
321
        # previous timesteps
322
        previous_timesteps = [x - 1 for x in self.TIMESTEPS]
323
        previous_timesteps[0] = self.TIMESTEPS.last()
324
325
        self.previous_timesteps = dict(zip(self.TIMESTEPS, previous_timesteps))
326
327
        # pyomo set for all flows in the energy system graph
328
        self.FLOWS = po.Set(
329
            initialize=self.flows.keys(), ordered=True, dimen=2
330
        )
331
332
        self.BIDIRECTIONAL_FLOWS = po.Set(
333
            initialize=[k for (k, v) in self.flows.items() if v.bidirectional],
334
            ordered=True,
335
            dimen=2,
336
            within=self.FLOWS,
337
        )
338
339
        self.UNIDIRECTIONAL_FLOWS = po.Set(
340
            initialize=[
341
                k for (k, v) in self.flows.items() if not v.bidirectional
342
            ],
343
            ordered=True,
344
            dimen=2,
345
            within=self.FLOWS,
346
        )
347
348
    def _add_parent_block_variables(self):
349
        """Add the parent block variables, which is the `flow` variable,
350
        indexed by FLOWS and TIMEINDEX."""
351
        self.flow = po.Var(self.FLOWS, self.TIMESTEPS, within=po.Reals)
352
353
        for o, i in self.FLOWS:
354
            if self.flows[o, i].nominal_capacity is not None:
355
                if self.flows[o, i].fix[self.TIMESTEPS.at(1)] is not None:
356
                    for t in self.TIMESTEPS:
357
                        self.flow[o, i, t].value = (
358
                            self.flows[o, i].fix[t]
359
                            * self.flows[o, i].nominal_capacity
360
                        )
361
                        self.flow[o, i, t].fix()
362
                else:
363
                    for t in self.TIMESTEPS:
364
                        self.flow[o, i, t].setub(
365
                            self.flows[o, i].max[t]
366
                            * self.flows[o, i].nominal_capacity
367
                        )
368
                    if not self.flows[o, i].nonconvex:
369
                        for t in self.TIMESTEPS:
370
                            self.flow[o, i, t].setlb(
371
                                self.flows[o, i].min[t]
372
                                * self.flows[o, i].nominal_capacity
373
                            )
374
                    elif (o, i) in self.UNIDIRECTIONAL_FLOWS:
375
                        for t in self.TIMESTEPS:
376
                            self.flow[o, i, t].setlb(0)
377
            else:
378
                if (o, i) in self.UNIDIRECTIONAL_FLOWS:
379
                    for t in self.TIMESTEPS:
380
                        self.flow[o, i, t].setlb(0)
381
382
    def _add_child_blocks(self):
383
        """Method to add the defined child blocks for components that have
384
        been grouped in the defined constraint groups. This collects all the
385
        constraints from the buses, components and flows blocks
386
        and adds them to the model.
387
        """
388
        for group in self._constraint_groups:
389
            block = group()
390
            self.add_component(str(block), block)
391
392
            # create constraints etc. related with block for all nodes
393
            # in the group
394
            block._create(group=self.es.groups.get(group))
395
396
    def _add_objective(self, sense=po.minimize, update=False):
397
        """Method to sum up all objective expressions from the child blocks
398
        that have been created. This method looks for `_objective_expression`
399
        attribute in the block definition and will call this method to add
400
        their return value to the objective function.
401
        """
402
        if update:
403
            self.del_component("objective")
404
405
        expr = 0
406
407
        for block in self.component_data_objects():
408
            if hasattr(block, "_objective_expression"):
409
                expr += block._objective_expression()
410
411
        self.objective = po.Objective(sense=sense, expr=expr)
412
413
    def receive_duals(self):
414
        """Method sets solver suffix to extract information about dual
415
        variables from solver. Shadow prices (duals) and reduced costs (rc) are
416
        set as attributes of the model.
417
        """
418
        # shadow prices
419
        self.dual = po.Suffix(direction=po.Suffix.IMPORT)
420
        # reduced costs
421
        self.rc = po.Suffix(direction=po.Suffix.IMPORT)
422
423
    def results(self):
424
        """Returns a nested dictionary of the results of this optimization.
425
        See the processing module for more information on results extraction.
426
        """
427
        return processing.results(self)
428
429
    def solve(
430
        self, solver="cbc", solver_io="lp", allow_nonoptimal=False, **kwargs
431
    ):
432
        r"""Takes care of communication with solver to solve the model.
433
434
        Parameters
435
        ----------
436
        solver : string
437
            solver to be used e.g. "cbc", "glpk", "gurobi", "cplex"
438
        solver_io : string
439
            pyomo solver interface file format: "lp", "python", "nl", etc.
440
        allow_nonoptimal : bool
441
            False: If no optimal solution is found, an error will be risen.
442
            True: If no optimal solution is found, there will be a warning.
443
        \**kwargs : keyword arguments
444
            Possible keys can be set see below:
445
446
        Other Parameters
447
        ----------------
448
        solve_kwargs : dict
449
            Other arguments for the pyomo.opt.SolverFactory.solve() method
450
            Example : {"tee":True}
451
        cmdline_options : dict
452
            Dictionary with command line options for solver e.g.
453
            {"mipgap":"0.01"} results in "--mipgap 0.01"
454
            \{"interior":" "} results in "--interior"
455
            \Gurobi solver takes numeric parameter values such as
456
            {"method": 2}
457
        """
458
        solve_kwargs = kwargs.get("solve_kwargs", {})
459
        solver_cmdline_options = kwargs.get("cmdline_options", {})
460
        opt = SolverFactory(solver, solver_io=solver_io)
461
462
        # set command line options
463
        options = opt.options
464
        for k in solver_cmdline_options:
465
            options[k] = solver_cmdline_options[k]
466
467
        solver_results = opt.solve(self, **solve_kwargs)
468
469
        status = solver_results.Solver.Status
470
        termination_condition = solver_results.Solver.Termination_condition
471
472
        self.es.results = solver_results
473
        self.solver_results = solver_results
474
475
        if status == "ok" and termination_condition == "optimal":
476
            logging.info("Optimization successful...")
477
        else:
478
            msg = (
479
                f"The solver did not return an optimal solution. "
480
                f"Instead the optimization ended with\n "
481
                f"      - status: {status}\n"
482
                f"       - termination condition: {termination_condition}"
483
            )
484
485
            if allow_nonoptimal:
486
                warnings.warn(
487
                    msg.format(status, termination_condition), UserWarning
488
                )
489
            else:
490
                raise RuntimeError(msg)
491
492
        return solver_results
493
494
    def relax_problem(self):
495
        """Relaxes integer variables to reals of optimization model self."""
496
        relaxer = RelaxIntegrality()
497
        relaxer._apply_to(self)
498
499
        return self
500
501
    def get_timestep_from_tsam_timestep(self, p, ik, g):
502
        """Return original timestep from cluster-based timestep"""
503
        t = (
504
            p * len(self.TIMESTEPS_IN_PERIOD[p])
505
            + ik * self.es.tsa_parameters[p]["timesteps"]
506
            + g
507
        )
508
        return t
509
510
    def get_cluster_index(self, cluster_type, offset):
511
        """
512
        Return cluster index for original or typical periods with or
513
        without offset
514
        """
515
        return [
516
            (p, k, t)
517
            for p in range(len(self.es.tsa_parameters))
518
            for k in range(len(self.es.tsa_parameters[p][cluster_type]))
519
            for t in range(self.es.tsa_parameters[p]["timesteps"] + offset)
520
        ]
521