Completed
Push — dev ( 91c66b...5706ca )
by Patrik
22s queued 17s
created

solph._models.Model.solve()   A

Complexity

Conditions 4

Size

Total Lines 54
Code Lines 20

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
eloc 20
dl 0
loc 54
rs 9.4
c 0
b 0
f 0
cc 4
nop 4

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
    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