Passed
Pull Request — dev (#821)
by Uwe
02:17
created

GenericCHPBlock.__init__()   A

Complexity

Conditions 1

Size

Total Lines 2
Code Lines 2

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
eloc 2
dl 0
loc 2
rs 10
c 0
b 0
f 0
cc 1
nop 3
1
# -*- coding: utf-8 -
2
3
"""
4
GenericCHP and associated individual constraints (blocks) and groupings.
5
6
SPDX-FileCopyrightText: Uwe Krien <[email protected]>
7
SPDX-FileCopyrightText: Simon Hilpert
8
SPDX-FileCopyrightText: Cord Kaldemeyer
9
SPDX-FileCopyrightText: Patrik Schönfeldt
10
SPDX-FileCopyrightText: FranziPl
11
SPDX-FileCopyrightText: jnnr
12
SPDX-FileCopyrightText: Stephan Günther
13
SPDX-FileCopyrightText: FabianTU
14
SPDX-FileCopyrightText: Johannes Röder
15
16
SPDX-License-Identifier: MIT
17
18
"""
19
20
import numpy as np
21
from oemof.network import network
22
from pyomo.core.base.block import ScalarBlock
23
from pyomo.environ import Binary
24
from pyomo.environ import Constraint
25
from pyomo.environ import NonNegativeReals
26
from pyomo.environ import Set
27
from pyomo.environ import Var
28
29
from oemof.solph._plumbing import sequence as solph_sequence
30
31
32
class GenericCHP(network.Transformer):
33
    r"""
34
    Component `GenericCHP` to model combined heat and power plants.
35
36
    Can be used to model (combined cycle) extraction or back-pressure turbines
37
    and used a mixed-integer linear formulation. Thus, it induces more
38
    computational effort than the `ExtractionTurbineCHP` for the
39
    benefit of higher accuracy.
40
41
    The full set of equations is described in:
42
    Mollenhauer, E., Christidis, A. & Tsatsaronis, G.
43
    Evaluation of an energy- and exergy-based generic modeling
44
    approach of combined heat and power plants
45
    Int J Energy Environ Eng (2016) 7: 167.
46
    https://doi.org/10.1007/s40095-016-0204-6
47
48
    For a general understanding of (MI)LP CHP representation, see:
49
    Fabricio I. Salgado, P.
50
    Short - Term Operation Planning on Cogeneration Systems: A Survey
51
    Electric Power Systems Research (2007)
52
    Electric Power Systems Research
53
    Volume 78, Issue 5, May 2008, Pages 835-848
54
    https://doi.org/10.1016/j.epsr.2007.06.001
55
56
    Note
57
    ----
58
    An adaption for the flow parameter `H_L_FG_share_max` has been made to
59
    set the flue gas losses at maximum heat extraction `H_L_FG_max` as share of
60
    the fuel flow `H_F` e.g. for combined cycle extraction turbines.
61
    The flow parameter `H_L_FG_share_min` can be used to set the flue gas
62
    losses at minimum heat extraction `H_L_FG_min` as share of
63
    the fuel flow `H_F` e.g. for motoric CHPs.
64
    The boolean component parameter `back_pressure` can be set to model
65
    back-pressure characteristics.
66
67
    Also have a look at the examples on how to use it.
68
69
    Parameters
70
    ----------
71
    fuel_input : dict
72
        Dictionary with key-value-pair of `oemof.Bus` and `oemof.Flow` object
73
        for the fuel input.
74
    electrical_output : dict
75
        Dictionary with key-value-pair of `oemof.Bus` and `oemof.Flow` object
76
        for the electrical output. Related parameters like `P_max_woDH` are
77
        passed as attributes of the `oemof.Flow` object.
78
    heat_output : dict
79
        Dictionary with key-value-pair of `oemof.Bus` and `oemof.Flow` object
80
        for the heat output. Related parameters like `Q_CW_min` are passed as
81
        attributes of the `oemof.Flow` object.
82
    Beta : list of numerical values
83
        Beta values in same dimension as all other parameters (length of
84
        optimization period).
85
    back_pressure : boolean
86
        Flag to use back-pressure characteristics. Set to `True` and
87
        `Q_CW_min` to zero for back-pressure turbines. See paper above for more
88
        information.
89
90
    Note
91
    ----
92
    The following sets, variables, constraints and objective parts are created
93
     * :py:class:`~oemof.solph.components._generic_chp.GenericCHPBlock`
94
95
    Examples
96
    --------
97
    >>> from oemof import solph
98
    >>> bel = solph.buses.Bus(label='electricityBus')
99
    >>> bth = solph.buses.Bus(label='heatBus')
100
    >>> bgas = solph.buses.Bus(label='commodityBus')
101
    >>> ccet = solph.components.GenericCHP(
102
    ...    label='combined_cycle_extraction_turbine',
103
    ...    fuel_input={bgas: solph.flows.Flow(
104
    ...        H_L_FG_share_max=[0.183])},
105
    ...    electrical_output={bel: solph.flows.Flow(
106
    ...        P_max_woDH=[155.946],
107
    ...        P_min_woDH=[68.787],
108
    ...        Eta_el_max_woDH=[0.525],
109
    ...        Eta_el_min_woDH=[0.444])},
110
    ...    heat_output={bth: solph.flows.Flow(
111
    ...        Q_CW_min=[10.552])},
112
    ...    Beta=[0.122], back_pressure=False)
113
    >>> type(ccet)
114
    <class 'oemof.solph.components._generic_chp.GenericCHP'>
115
    """
116
117
    def __init__(self, *args, **kwargs):
118
        super().__init__(*args, **kwargs)
119
120
        self.fuel_input = kwargs.get("fuel_input")
121
        self.electrical_output = kwargs.get("electrical_output")
122
        self.heat_output = kwargs.get("heat_output")
123
        self.Beta = solph_sequence(kwargs.get("Beta"))
124
        self.back_pressure = kwargs.get("back_pressure")
125
        self._alphas = None
126
127
        # map specific flows to standard API
128
        fuel_bus = list(self.fuel_input.keys())[0]
129
        fuel_flow = list(self.fuel_input.values())[0]
130
        fuel_bus.outputs.update({self: fuel_flow})
131
132
        self.outputs.update(kwargs.get("electrical_output"))
133
        self.outputs.update(kwargs.get("heat_output"))
134
135
    def _calculate_alphas(self):
136
        """
137
        Calculate alpha coefficients.
138
139
        A system of linear equations is created from passed capacities and
140
        efficiencies and solved to calculate both coefficients.
141
        """
142
        alphas = [[], []]
143
144
        eb = list(self.electrical_output.keys())[0]
145
146
        attrs = [
147
            self.electrical_output[eb].P_min_woDH,
148
            self.electrical_output[eb].Eta_el_min_woDH,
149
            self.electrical_output[eb].P_max_woDH,
150
            self.electrical_output[eb].Eta_el_max_woDH,
151
        ]
152
153
        length = [len(a) for a in attrs if not isinstance(a, (int, float))]
154
        max_length = max(length)
155
156
        if all(len(a) == max_length for a in attrs):
157
            if max_length == 0:
158
                max_length += 1  # increment dimension for scalars from 0 to 1
159
            for i in range(0, max_length):
160
                A = np.array(
161
                    [
162
                        [1, self.electrical_output[eb].P_min_woDH[i]],
163
                        [1, self.electrical_output[eb].P_max_woDH[i]],
164
                    ]
165
                )
166
                b = np.array(
167
                    [
168
                        self.electrical_output[eb].P_min_woDH[i]
169
                        / self.electrical_output[eb].Eta_el_min_woDH[i],
170
                        self.electrical_output[eb].P_max_woDH[i]
171
                        / self.electrical_output[eb].Eta_el_max_woDH[i],
172
                    ]
173
                )
174
                x = np.linalg.solve(A, b)
175
                alphas[0].append(x[0])
176
                alphas[1].append(x[1])
177
        else:
178
            error_message = (
179
                "Attributes to calculate alphas "
180
                + "must be of same dimension."
181
            )
182
            raise ValueError(error_message)
183
184
        self._alphas = alphas
185
186
    @property
187
    def alphas(self):
188
        """Compute or return the _alphas attribute."""
189
        if self._alphas is None:
190
            self._calculate_alphas()
191
        return self._alphas
192
193
    def constraint_group(self):
194
        return GenericCHPBlock
195
196
197
class GenericCHPBlock(ScalarBlock):
198
    r"""
199
    Block for the relation of the :math:`n` nodes with
200
    type class:`.GenericCHP`.
201
202
    **The following constraints are created:**
203
204
    .. _GenericCHP-equations1-10:
205
206
    .. math::
207
        &
208
        (1)\qquad \dot{H}_F(t) = fuel\ input \\
209
        &
210
        (2)\qquad \dot{Q}(t) = heat\ output \\
211
        &
212
        (3)\qquad P_{el}(t) = power\ output\\
213
        &
214
        (4)\qquad \dot{H}_F(t) = \alpha_0(t) \cdot Y(t) + \alpha_1(t) \cdot
215
        P_{el,woDH}(t)\\
216
        &
217
        (5)\qquad \dot{H}_F(t) = \alpha_0(t) \cdot Y(t) + \alpha_1(t) \cdot
218
        ( P_{el}(t) + \beta \cdot \dot{Q}(t) )\\
219
        &
220
        (6)\qquad \dot{H}_F(t) \leq Y(t) \cdot
221
        \frac{P_{el, max, woDH}(t)}{\eta_{el,max,woDH}(t)}\\
222
        &
223
        (7)\qquad \dot{H}_F(t) \geq Y(t) \cdot
224
        \frac{P_{el, min, woDH}(t)}{\eta_{el,min,woDH}(t)}\\
225
        &
226
        (8)\qquad \dot{H}_{L,FG,max}(t) = \dot{H}_F(t) \cdot
227
        \dot{H}_{L,FG,sharemax}(t)\\
228
        &
229
        (9)\qquad \dot{H}_{L,FG,min}(t) = \dot{H}_F(t) \cdot
230
        \dot{H}_{L,FG,sharemin}(t)\\
231
        &
232
        (10)\qquad P_{el}(t) + \dot{Q}(t) + \dot{H}_{L,FG,max}(t) +
233
        \dot{Q}_{CW, min}(t) \cdot Y(t) = / \leq \dot{H}_F(t)\\
234
235
    where :math:`= / \leq` depends on the CHP being back pressure or not.
236
237
    The coefficients :math:`\alpha_0` and :math:`\alpha_1`
238
    can be determined given the efficiencies maximal/minimal load:
239
240
    .. math::
241
        &
242
        \eta_{el,max,woDH}(t) = \frac{P_{el,max,woDH}(t)}{\alpha_0(t)
243
        \cdot Y(t) + \alpha_1(t) \cdot P_{el,max,woDH}(t)}\\
244
        &
245
        \eta_{el,min,woDH}(t) = \frac{P_{el,min,woDH}(t)}{\alpha_0(t)
246
        \cdot Y(t) + \alpha_1(t) \cdot P_{el,min,woDH}(t)}\\
247
248
249
    **For the attribute** :math:`\dot{H}_{L,FG,min}` **being not None**,
250
    e.g. for a motoric CHP, **the following is created:**
251
252
        **Constraint:**
253
254
    .. _GenericCHP-equations11:
255
256
    .. math::
257
        &
258
        (11)\qquad P_{el}(t) + \dot{Q}(t) + \dot{H}_{L,FG,min}(t) +
259
        \dot{Q}_{CW, min}(t) \cdot Y(t) \geq \dot{H}_F(t)\\[10pt]
260
261
    The symbols used are defined as follows (with Variables (V) and Parameters (P)):
262
263
    =============================== ======================= ==== =======================
264
    math. symbol                    attribute               type explanation
265
    =============================== ======================= ==== =======================
266
    :math:`\dot{H}_{F}`             `H_F[n,t]`              V    input of enthalpy
267
                                                                         through fuel input
268
    :math:`P_{el}`                  `P[n,t]`                V    provided
269
                                                                         electric power
270
    :math:`P_{el,woDH}`             `P_woDH[n,t]`           V    electric power without
271
                                                                         district heating
272
    :math:`P_{el,min,woDH}`         `P_min_woDH[n,t]`       P    min. electric power
273
                                                                         without district heating
274
    :math:`P_{el,max,woDH}`         `P_max_woDH[n,t]`       P    max. electric power
275
                                                                         without district heating
276
    :math:`\dot{Q}`                 `Q[n,t]`                V    provided heat
277
278
    :math:`\dot{Q}_{CW, min}`       `Q_CW_min[n,t]`         P    minimal therm. condenser
279
                                                                         load to cooling water
280
    :math:`\dot{H}_{L,FG,min}`      `H_L_FG_min[n,t]`       V    flue gas enthalpy loss
281
                                                                         at min heat extraction
282
    :math:`\dot{H}_{L,FG,max}`      `H_L_FG_max[n,t]`       V    flue gas enthalpy loss
283
                                                                         at max heat extraction
284
    :math:`\dot{H}_{L,FG,sharemin}` `H_L_FG_share_min[n,t]` P    share of flue gas loss
285
                                                                         at min heat extraction
286
    :math:`\dot{H}_{L,FG,sharemax}` `H_L_FG_share_max[n,t]` P    share of flue gas loss
287
                                                                         at max heat extraction
288
    :math:`Y`                       `Y[n,t]`                V    status variable
289
                                                                         on/off
290
    :math:`\alpha_0`                `n.alphas[0][n,t]`      P    coefficient
291
                                                                         describing efficiency
292
    :math:`\alpha_1`                `n.alphas[1][n,t]`      P    coefficient
293
                                                                         describing efficiency
294
    :math:`\beta`                   `Beta[n,t]`             P    power loss index
295
296
    :math:`\eta_{el,min,woDH}`      `Eta_el_min_woDH[n,t]`  P    el. eff. at min. fuel
297
                                                                         flow w/o distr. heating
298
    :math:`\eta_{el,max,woDH}`      `Eta_el_max_woDH[n,t]`  P    el. eff. at max. fuel
299
                                                                         flow w/o distr. heating
300
    =============================== ======================= ==== =======================
301
302
    """  # noqa: E501
303
    CONSTRAINT_GROUP = True
304
305
    def __init__(self, *args, **kwargs):
306
        super().__init__(*args, **kwargs)
307
308
    def _create(self, group=None):
309
        """
310
        Create constraints for GenericCHPBlock.
311
312
        Parameters
313
        ----------
314
        group : list
315
            List containing `GenericCHP` objects.
316
            e.g. groups=[ghcp1, gchp2,..]
317
        """
318
        m = self.parent_block()
319
320
        if group is None:
321
            return None
322
323
        self.GENERICCHPS = Set(initialize=[n for n in group])
324
325
        # variables
326
        self.H_F = Var(self.GENERICCHPS, m.TIMESTEPS, within=NonNegativeReals)
327
        self.H_L_FG_max = Var(
328
            self.GENERICCHPS, m.TIMESTEPS, within=NonNegativeReals
329
        )
330
        self.H_L_FG_min = Var(
331
            self.GENERICCHPS, m.TIMESTEPS, within=NonNegativeReals
332
        )
333
        self.P_woDH = Var(
334
            self.GENERICCHPS, m.TIMESTEPS, within=NonNegativeReals
335
        )
336
        self.P = Var(self.GENERICCHPS, m.TIMESTEPS, within=NonNegativeReals)
337
        self.Q = Var(self.GENERICCHPS, m.TIMESTEPS, within=NonNegativeReals)
338
        self.Y = Var(self.GENERICCHPS, m.TIMESTEPS, within=Binary)
339
340
        # constraint rules
341
        def _H_flow_rule(block, n, t):
342
            """Link fuel consumption to component inflow."""
343
            expr = 0
344
            expr += self.H_F[n, t]
345
            expr += -m.flow[list(n.fuel_input.keys())[0], n, t]
346
            return expr == 0
347
348
        self.H_flow = Constraint(
349
            self.GENERICCHPS, m.TIMESTEPS, rule=_H_flow_rule
350
        )
351
352
        def _Q_flow_rule(block, n, t):
353
            """Link heat flow to component outflow."""
354
            expr = 0
355
            expr += self.Q[n, t]
356
            expr += -m.flow[n, list(n.heat_output.keys())[0], t]
357
            return expr == 0
358
359
        self.Q_flow = Constraint(
360
            self.GENERICCHPS, m.TIMESTEPS, rule=_Q_flow_rule
361
        )
362
363
        def _P_flow_rule(block, n, t):
364
            """Link power flow to component outflow."""
365
            expr = 0
366
            expr += self.P[n, t]
367
            expr += -m.flow[n, list(n.electrical_output.keys())[0], t]
368
            return expr == 0
369
370
        self.P_flow = Constraint(
371
            self.GENERICCHPS, m.TIMESTEPS, rule=_P_flow_rule
372
        )
373
374
        def _H_F_1_rule(block, n, t):
375
            """Set P_woDH depending on H_F."""
376
            expr = 0
377
            expr += -self.H_F[n, t]
378
            expr += n.alphas[0][t] * self.Y[n, t]
379
            expr += n.alphas[1][t] * self.P_woDH[n, t]
380
            return expr == 0
381
382
        self.H_F_1 = Constraint(
383
            self.GENERICCHPS, m.TIMESTEPS, rule=_H_F_1_rule
384
        )
385
386
        def _H_F_2_rule(block, n, t):
387
            """Determine relation between H_F, P and Q."""
388
            expr = 0
389
            expr += -self.H_F[n, t]
390
            expr += n.alphas[0][t] * self.Y[n, t]
391
            expr += n.alphas[1][t] * (self.P[n, t] + n.Beta[t] * self.Q[n, t])
392
            return expr == 0
393
394
        self.H_F_2 = Constraint(
395
            self.GENERICCHPS, m.TIMESTEPS, rule=_H_F_2_rule
396
        )
397
398
        def _H_F_3_rule(block, n, t):
399
            """Set upper value of operating range via H_F."""
400
            expr = 0
401
            expr += self.H_F[n, t]
402
            expr += -self.Y[n, t] * (
403
                list(n.electrical_output.values())[0].P_max_woDH[t]
404
                / list(n.electrical_output.values())[0].Eta_el_max_woDH[t]
405
            )
406
            return expr <= 0
407
408
        self.H_F_3 = Constraint(
409
            self.GENERICCHPS, m.TIMESTEPS, rule=_H_F_3_rule
410
        )
411
412
        def _H_F_4_rule(block, n, t):
413
            """Set lower value of operating range via H_F."""
414
            expr = 0
415
            expr += self.H_F[n, t]
416
            expr += -self.Y[n, t] * (
417
                list(n.electrical_output.values())[0].P_min_woDH[t]
418
                / list(n.electrical_output.values())[0].Eta_el_min_woDH[t]
419
            )
420
            return expr >= 0
421
422
        self.H_F_4 = Constraint(
423
            self.GENERICCHPS, m.TIMESTEPS, rule=_H_F_4_rule
424
        )
425
426
        def _H_L_FG_max_rule(block, n, t):
427
            """Set max. flue gas loss as share fuel flow share."""
428
            expr = 0
429
            expr += -self.H_L_FG_max[n, t]
430
            expr += (
431
                self.H_F[n, t]
432
                * list(n.fuel_input.values())[0].H_L_FG_share_max[t]
433
            )
434
            return expr == 0
435
436
        self.H_L_FG_max_def = Constraint(
437
            self.GENERICCHPS, m.TIMESTEPS, rule=_H_L_FG_max_rule
438
        )
439
440
        def _Q_max_res_rule(block, n, t):
441
            """Set maximum Q depending on fuel and electrical flow."""
442
            expr = 0
443
            expr += self.P[n, t] + self.Q[n, t] + self.H_L_FG_max[n, t]
444
            expr += list(n.heat_output.values())[0].Q_CW_min[t] * self.Y[n, t]
445
            expr += -self.H_F[n, t]
446
            # back-pressure characteristics or one-segment model
447
            if n.back_pressure is True:
448
                return expr == 0
449
            else:
450
                return expr <= 0
451
452
        self.Q_max_res = Constraint(
453
            self.GENERICCHPS, m.TIMESTEPS, rule=_Q_max_res_rule
454
        )
455
456
        def _H_L_FG_min_rule(block, n, t):
457
            """Set min. flue gas loss as fuel flow share."""
458
            # minimum flue gas losses e.g. for motoric CHPs
459
            if getattr(
460
                list(n.fuel_input.values())[0], "H_L_FG_share_min", None
461
            ):
462
                expr = 0
463
                expr += -self.H_L_FG_min[n, t]
464
                expr += (
465
                    self.H_F[n, t]
466
                    * list(n.fuel_input.values())[0].H_L_FG_share_min[t]
467
                )
468
                return expr == 0
469
            else:
470
                return Constraint.Skip
471
472
        self.H_L_FG_min_def = Constraint(
473
            self.GENERICCHPS, m.TIMESTEPS, rule=_H_L_FG_min_rule
474
        )
475
476
        def _Q_min_res_rule(block, n, t):
477
            """Set minimum Q depending on fuel and eletrical flow."""
478
            # minimum restriction for heat flows e.g. for motoric CHPs
479
            if getattr(
480
                list(n.fuel_input.values())[0], "H_L_FG_share_min", None
481
            ):
482
                expr = 0
483
                expr += self.P[n, t] + self.Q[n, t] + self.H_L_FG_min[n, t]
484
                expr += (
485
                    list(n.heat_output.values())[0].Q_CW_min[t] * self.Y[n, t]
486
                )
487
                expr += -self.H_F[n, t]
488
                return expr >= 0
489
            else:
490
                return Constraint.Skip
491
492
        self.Q_min_res = Constraint(
493
            self.GENERICCHPS, m.TIMESTEPS, rule=_Q_min_res_rule
494
        )
495
496
    def _objective_expression(self):
497
        r"""Objective expression for generic CHPs with no investment.
498
499
        Note: This adds nothing as variable costs are already
500
        added in the Block :class:`FlowBlock`.
501
        """
502
        if not hasattr(self, "GENERICCHPS"):
503
            return 0
504
505
        return 0
506