2d_constraints.gms 157 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
$ontext
This file is part of Backbone.

Backbone is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

Backbone is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with Backbone.  If not, see <http://www.gnu.org/licenses/>.
$offtext

18

19
* =============================================================================
20
* --- Constraint Equation Definitions -----------------------------------------
21
22
23
24
* =============================================================================

* --- Energy Balance ----------------------------------------------------------

25
26
27
q_balance(gn(grid, node), msft(m, s, f, t)) // Energy/power balance dynamics solved using implicit Euler discretization
    ${  not p_gn(grid, node, 'boundAll')
        } ..
28
29
30
31

    // The left side of the equation is the change in the state (will be zero if the node doesn't have a state)
    + p_gn(grid, node, 'energyStoredPerUnitOfState')${gn_state(grid, node)} // Unit conversion between v_state of a particular node and energy variables (defaults to 1, but can have node based values if e.g. v_state is in Kelvins and each node has a different heat storage capacity)
        * [
32
33
            + v_state(grid, node, s, f+df_central(f,t), t)                   // The difference between current
            - v_state(grid, node, s+ds_state(grid,node,s,t), f+df(f,t+dt(t)), t+dt(t))       // ... and previous state of the node
34
35
36
37
38
39
40
41
            ]

    =E=

    // The right side of the equation contains all the changes converted to energy terms
    + p_stepLength(m, f, t) // Multiply with the length of the timestep to convert power into energy
        * (
            // Self discharge out of the model boundaries
42
            - p_gn(grid, node, 'selfDischargeLoss')${ gn_state(grid, node) }
43
                * v_state(grid, node, s, f+df_central(f,t), t) // The current state of the node
44
45

            // Energy diffusion from this node to neighbouring nodes
46
            - sum(gnn_state(grid, node, to_node),
47
                + p_gnn(grid, node, to_node, 'diffCoeff')
48
                    * v_state(grid, node, s, f+df_central(f,t), t)
49
50
51
                ) // END sum(to_node)

            // Energy diffusion from neighbouring nodes to this node
52
            + sum(gnn_state(grid, from_node, node),
53
                + p_gnn(grid, from_node, node, 'diffCoeff')
54
                    * v_state(grid, from_node, s, f+df_central(f,t), t) // Incoming diffusion based on the state of the neighbouring node
55
56
57
                ) // END sum(from_node)

            // Controlled energy transfer, applies when the current node is on the left side of the connection
58
            - sum(gn2n_directional(grid, node, node_),
59
                + (1 - p_gnn(grid, node, node_, 'transferLoss')) // Reduce transfer losses
60
                    * v_transfer(grid, node, node_, s, f, t)
61
                + p_gnn(grid, node, node_, 'transferLoss') // Add transfer losses back if transfer is from this node to another node
62
                    * v_transferRightward(grid, node, node_, s, f, t)
63
64
65
                ) // END sum(node_)

            // Controlled energy transfer, applies when the current node is on the right side of the connection
66
            + sum(gn2n_directional(grid, node_, node),
67
                + v_transfer(grid, node_, node, s, f, t)
68
                - p_gnn(grid, node_, node, 'transferLoss') // Reduce transfer losses if transfer is from another node to this node
69
                    * v_transferRightward(grid, node_, node, s, f, t)
70
71
72
73
                ) // END sum(node_)

            // Interactions between the node and its units
            + sum(gnuft(grid, node, unit, f, t),
74
                + v_gen(grid, node, unit, s, f, t) // Unit energy generation and consumption
75
                )
76
77

            // Spilling energy out of the endogenous grids in the model
78
            - v_spill(grid, node, s, f, t)${node_spill(node)}
79
80

            // Power inflow and outflow timeseries to/from the node
81
            + ts_influx_(grid, node, f, t, s)   // Incoming (positive) and outgoing (negative) absolute value time series
82
83

            // Dummy generation variables, for feasibility purposes
84
85
            + vq_gen('increase', grid, node, s, f, t) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
            - vq_gen('decrease', grid, node, s, f, t) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
86
    ) // END * p_stepLength
87
;
88
89

* --- Reserve Demand ----------------------------------------------------------
90
91
// NOTE! Currently, there are multiple identical instances of the reserve balance equation being generated for each forecast branch even when the reserves are committed and identical between the forecasts.
// NOTE! This could be solved by formulating a new "ft_reserves" set to cover only the relevant forecast-time steps, but it would possibly make the reserves even more confusing.
92

93
94
q_resDemand(restypeDirectionGroup(restype, up_down, group), sft(s, f, t))
    ${  ord(t) < tSolveFirst + p_groupReserves(group, restype, 'reserve_length')
95
        and not [ restypeReleasedForRealization(restype)
96
                  and sft_realized(s, f, t)]
97
        } ..
98

99
100
    // Reserve provision by capable units on this group
    + sum(gnuft(grid, node, unit, f, t)${gnGroup(grid, node, group) and nuRescapable(restype, up_down, node, unit)},
101
        + v_reserve(restype, up_down, node, unit, s, f+df_reserves(node, restype, f, t), t)
102
103
104
105
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                ] // END * v_reserve
106
        ) // END sum(gnuft)
107

108
    // Reserve provision from other reserve categories when they can be shared
109
    + sum((gnuft(grid, node, unit, f, t), restype_)${gnGroup(grid, node, group) and p_nuRes2Res(node, unit, restype_, up_down, restype)},
110
        + v_reserve(restype_, up_down, node, unit, s, f+df_reserves(node, restype_, f, t), t)
111
            * p_nuRes2Res(node, unit, restype_, up_down, restype)
112
113
114
115
116
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                    * p_nuReserves(node, unit, restype_, 'reserveReliability')
                ] // END * v_reserve
117
        ) // END sum(gnuft)
118

119
120
121
122
123
    // Reserve provision to this group via transfer links
    + sum(gn2n_directional(grid, node_, node)${gnGroup(grid, node, group)
                                               and not gnGroup(grid, node_, group)
                                               and restypeDirectionNodeNode(restype, up_down, node_, node)
                                               },
124
        + (1 - p_gnn(grid, node_, node, 'transferLoss') )
125
            * v_resTransferRightward(restype, up_down, node_, node, s, f+df_reserves(node_, restype, f, t), t) // Reserves from another node - reduces the need for reserves in the node
126
        ) // END sum(gn2n_directional)
127
128
129
130
    + sum(gn2n_directional(grid, node, node_)${gnGroup(grid, node, group)
                                               and not gnGroup(grid, node_, group)
                                               and restypeDirectionNodeNode(restype, up_down, node_, node)
                                               },
131
        + (1 - p_gnn(grid, node, node_, 'transferLoss') )
132
            * v_resTransferLeftward(restype, up_down, node, node_, s, f+df_reserves(node_, restype, f, t), t) // Reserves from another node - reduces the need for reserves in the node
133
134
135
136
137
        ) // END sum(gn2n_directional)

    =G=

    // Demand for reserves
138
139
    + ts_reserveDemand(restype, up_down, group, f, t)${p_groupReserves(group, restype, 'use_time_series')}
    + p_groupReserves(group, restype, up_down)${not p_groupReserves(group, restype, 'use_time_series')}
140

141
    // Reserve demand increase because of units
142
    + sum(gnuft(grid, node, unit, f, t)${gnGroup(grid, node, group) and p_nuReserves(node, unit, restype, 'reserve_increase_ratio')}, // Could be better to have 'reserve_increase_ratio' separately for up and down directions
143
        + sum(gnu(grid, node, unit), v_gen(grid, node, unit, s, f, t)) // Reserve sets and variables are currently lacking the grid dimension...
144
145
146
            * p_nuReserves(node, unit, restype, 'reserve_increase_ratio')
        ) // END sum(nuft)

147
148
149
150
151
    // Reserve provisions to other groups via transfer links
    + sum(gn2n_directional(grid, node, node_)${gnGroup(grid, node, group)
                                               and not gnGroup(grid, node_, group)
                                               and restypeDirectionNodeNode(restype, up_down, node, node_)
                                               },   // If trasferring reserves to another node, increase your own reserves by same amount
152
        + v_resTransferRightward(restype, up_down, node, node_, s, f+df_reserves(node, restype, f, t), t)
153
        ) // END sum(gn2n_directional)
154
155
156
157
    + sum(gn2n_directional(grid, node_, node)${gnGroup(grid, node, group)
                                               and not gnGroup(grid, node_, group)
                                               and restypeDirectionNodeNode(restype, up_down, node, node_)
                                               },   // If trasferring reserves to another node, increase your own reserves by same amount
158
        + v_resTransferLeftward(restype, up_down, node_, node, s, f+df_reserves(node, restype, f, t), t)
159
160
161
        ) // END sum(gn2n_directional)

    // Reserve demand feasibility dummy variables
162
163
    - vq_resDemand(restype, up_down, group, s, f+df_reservesGroup(group, restype, f, t), t)
    - vq_resMissing(restype, up_down, group, s, f+df_reservesGroup(group, restype, f, t), t)${ft_reservesFixed(group, restype, f+df_reservesGroup(group, restype, f, t), t)}
164
;
165

166
167
168
169
* --- N-1 Reserve Demand ----------------------------------------------------------
// NOTE! Currently, there are multiple identical instances of the reserve balance equation being generated for each forecast branch even when the reserves are committed and identical between the forecasts.
// NOTE! This could be solved by formulating a new "ft_reserves" set to cover only the relevant forecast-time steps, but it would possibly make the reserves even more confusing.

170
q_resDemandLargestInfeedUnit(grid, restypeDirectionNode(restype, 'up', node), unit_fail(unit_), sft(s, f, t))
171
    ${  ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')
172
        and gn(grid, node)
173
174
175
        and not [ restypeReleasedForRealization(restype)
            and ft_realized(f, t)
            ]
176
        and p_nuReserves(node, unit_, restype, 'portion_of_infeed_to_reserve')
177
        } ..
178

179
180
    // Reserve provision by capable units on this node excluding the failing one
    + sum(nuft(node, unit, f, t)${nuRescapable(restype, 'up', node, unit) and (ord(unit_) ne ord(unit))},
181
        + v_reserve(restype, 'up', node, unit, s, f+df_reserves(node, restype, f, t), t)
182
183
184
185
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                ] // END * v_reserve
186
187
188
        ) // END sum(nuft)

    // Reserve provision from other reserve categories when they can be shared
189
    + sum((nuft(node, unit, f, t), restype_)${p_nuRes2Res(node, unit, restype_, 'up', restype)},
190
        + v_reserve(restype_, 'up', node, unit, s, f+df_reserves(node, restype_, f, t), t)
191
            * p_nuRes2Res(node, unit, restype_, 'up', restype)
192
193
194
195
196
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                    * p_nuReserves(node, unit, restype_, 'reserveReliability')
                ] // END * v_reserve
197
198
199
200
201
        ) // END sum(nuft)

    // Reserve provision to this node via transfer links
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, 'up', node_, node)},
        + (1 - p_gnn(grid, node_, node, 'transferLoss') )
202
            * v_resTransferRightward(restype, 'up', node_, node, s, f+df_reserves(node_, restype, f, t), t) // Reserves from another node - reduces the need for reserves in the node
203
204
205
        ) // END sum(gn2n_directional)
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, 'up', node_, node)},
        + (1 - p_gnn(grid, node, node_, 'transferLoss') )
206
            * v_resTransferLeftward(restype, 'up', node, node_, s, f+df_reserves(node_, restype, f, t), t) // Reserves from another node - reduces the need for reserves in the node
207
208
209
210
        ) // END sum(gn2n_directional)

    =G=

211
    // Demand for reserves due to a large unit that could fail
212
    + v_gen(grid,node,unit_,s,f,t) * p_nuReserves(node, unit_, restype, 'portion_of_infeed_to_reserve')
213
214

    // Reserve provisions to another nodes via transfer links
215
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, 'up', node, node_)},   // If trasferring reserves to another node, increase your own reserves by same amount
216
        + v_resTransferRightward(restype, 'up', node, node_, s, f+df_reserves(node, restype, f, t), t)
217
        ) // END sum(gn2n_directional)
218
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, 'up', node, node_)},   // If trasferring reserves to another node, increase your own reserves by same amount
219
        + v_resTransferLeftward(restype, 'up', node_, node, s, f+df_reserves(node, restype, f, t), t)
220
221
222
        ) // END sum(gn2n_directional)

    // Reserve demand feasibility dummy variables
223
224
*    - vq_resDemand(restype, 'up', node, s, f+df_reserves(node, restype, f, t), t)
*    - vq_resMissing(restype, 'up', node, s, f+df_reserves(node, restype, f, t), t)${ft_reservesFixed(node, restype, f+df_reserves(node, restype, f, t), t)}
225
;
226

227
228
229
230
231
q_rateOfChangeOfFrequencyUnit(group, unit_fail(unit_), sft(s, f, t))
    ${  p_groupPolicy(group, 'defaultFrequency')
        and p_groupPolicy(group, 'ROCOF')
        and uft(unit_, f, t)
        and sum(gnu_output(grid, node, unit_)${gnGroup(grid, node, group)}, 1) // only units with output capacity 'inside the group'
232
233
        } ..

234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
    // Kinetic/rotational energy in the system
    + p_groupPolicy(group, 'ROCOF')*2
        * [
            + sum(gnu_output(grid, node, unit)${   ord(unit) ne ord(unit_)
                                                   and gnGroup(grid, node, group)
                                                   and gnuft(grid, node, unit, f, t)
                                                   },
                + p_gnu(grid, node, unit, 'inertia')
                    * p_gnu(grid ,node, unit, 'unitSizeMVA')
                    * [
                        + v_online_LP(unit, s, f+df_central(f,t), t)
                            ${uft_onlineLP(unit, f, t)}
                        + v_online_MIP(unit, s, f+df_central(f,t), t)
                            ${uft_onlineMIP(unit, f, t)}
                        + v_gen(grid, node, unit, s, f, t)${not uft_online(unit, f, t)}
                            / p_gnu(grid, node, unit, 'unitSizeGen')
                        ] // * p_gnu
                ) // END sum(gnu_output)
            ] // END * p_groupPolicy
253
254
255

    =G=

256
257
258
259
260
261
    // Demand for kinetic/rotational energy due to a large unit that could fail
    + p_groupPolicy(group, 'defaultFrequency')
        * sum(gnu_output(grid, node, unit_)${   gnGroup(grid, node, group)
                                                },
            + v_gen(grid, node, unit_ , s, f, t)
            ) // END sum(gnu_output)
262
;
263
264
265
266
267
268
269
270
271

q_rateOfChangeOfFrequencyTransfer(group, gn(grid, node_), node_fail, sft(s, f, t))
    ${  p_groupPolicy(group, 'defaultFrequency')
        and p_groupPolicy(group, 'ROCOF')
        and gnGroup(grid, node_, group) // only interconnectors where one end is 'inside the group'
        and not gnGroup(grid, node_fail, group) // and the other end is 'outside the group'
        and [ p_gnn(grid, node_, node_fail, 'portion_of_transfer_to_reserve')
              or p_gnn(grid, node_fail, node_, 'portion_of_transfer_to_reserve')
              ]
272
273
        } ..

274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
    // Kinetic/rotational energy in the system
    + p_groupPolicy(group, 'ROCOF')*2
        * [
            + sum(gnu_output(grid, node, unit)${   gnGroup(grid, node, group)
                                                   and gnuft(grid, node, unit, f, t)
                                                   },
                + p_gnu(grid, node, unit, 'inertia')
                    * p_gnu(grid ,node, unit, 'unitSizeMVA')
                    * [
                        + v_online_LP(unit, s, f+df_central(f,t), t)
                            ${uft_onlineLP(unit, f, t)}
                        + v_online_MIP(unit, s, f+df_central(f,t), t)
                            ${uft_onlineMIP(unit, f, t)}
                        + v_gen(grid, node, unit, s, f, t)${not uft_online(unit, f, t)}
                            / p_gnu(grid, node, unit, 'unitSizeGen')
                        ] // * p_gnu
                ) // END sum(gnu_output)
            ] // END * p_groupPolicy
292
293
294

    =G=

295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
    // Demand for kinetic/rotational energy due to a large interconnector that could fail
    + p_groupPolicy(group, 'defaultFrequency')
        * [
            // Loss of import due to potential interconnector failures
            + p_gnn(grid, node_fail, node_, 'portion_of_transfer_to_reserve')
                * v_transferRightward(grid, node_fail, node_, s, f, t)
                * (1 - p_gnn(grid, node_fail, node_, 'transferLoss') )
            + p_gnn(grid, node_, node_fail, 'portion_of_transfer_to_reserve')
                * v_transferLeftward(grid, node_, node_fail, s, f, t)
                * (1 - p_gnn(grid, node_, node_fail, 'transferLoss') )
            // Loss of export due to potential interconnector failures
            + p_gnn(grid, node_fail, node_, 'portion_of_transfer_to_reserve')
                * v_transferLeftward(grid, node_fail, node_, s, f, t)
            + p_gnn(grid, node_, node_fail, 'portion_of_transfer_to_reserve')
                * v_transferRightward(grid, node_, node_fail, s, f, t)
            ] // END * p_groupPolicy
311
;
312

313
* --- N-1 Upward reserve demand due to a possibility that an interconnector that is transferring power to the node fails -------------------------------------------------
314
315
316
// NOTE! Currently, there are multiple identical instances of the reserve balance equation being generated for each forecast branch even when the reserves are committed and identical between the forecasts.
// NOTE! This could be solved by formulating a new "ft_reserves" set to cover only the relevant forecast-time steps, but it would possibly make the reserves even more confusing.

ran li's avatar
ran li committed
317
q_resDemandLargestInfeedTransfer(grid, restypeDirectionNode(restype, 'up', node), node_fail, sft(s, f, t))
318
319
320
    ${  ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')
        and not [ restypeReleasedForRealization(restype)
                  and sft_realized(s, f, t)]
321
322
323
        and [ p_gnn(grid, node, node_fail, 'portion_of_transfer_to_reserve')
              or p_gnn(grid, node_fail, node, 'portion_of_transfer_to_reserve')
              ]
ran li's avatar
ran li committed
324
        and p_nReserves(node, restype, 'LossOfTrans')
325
326
327
        } ..

    // Reserve provision by capable units on this node
ran li's avatar
ran li committed
328
329
    + sum(nuft(node, unit, f, t)${nuRescapable(restype, 'up', node, unit)},
        + v_reserve(restype, 'up', node, unit, s, f+df_reserves(node, restype, f, t), t)
330
331
332
333
334
335
336
337
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                ] // END * v_reserve
        ) // END sum(nuft)

    // Reserve provision from other reserve categories when they can be shared
    + sum((nuft(node, unit, f, t), restype_)${p_nuRes2Res(node, unit, restype_, 'up', restype)},
ran li's avatar
ran li committed
338
339
        + v_reserve(restype_, 'up', node, unit, s, f+df_reserves(node, restype_, f, t), t)
            * p_nuRes2Res(node, unit, restype_, 'up', restype)
340
341
342
343
344
345
346
347
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                    * p_nuReserves(node, unit, restype_, 'reserveReliability')
                ] // END * v_reserve
        ) // END sum(nuft)

    // Reserve provision to this node via transfer links
348
    // SHOULD THE node_fail BE EXCLUDED?
ran li's avatar
ran li committed
349
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, 'up', node_, node)},
350
        + (1 - p_gnn(grid, node_, node, 'transferLoss') )
ran li's avatar
ran li committed
351
            * v_resTransferRightward(restype, 'up', node_, node, s, f+df_reserves(node_, restype, f, t), t)
352
            * [1$(not node_(node_fail)) + p_gnn(grid, node_, node, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
353
        ) // END sum(gn2n_directional)
ran li's avatar
ran li committed
354
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, 'up', node_, node)},
355
        + (1 - p_gnn(grid, node, node_, 'transferLoss') )
ran li's avatar
ran li committed
356
            * v_resTransferLeftward(restype, 'up', node, node_, s, f+df_reserves(node_, restype, f, t), t)
357
            * [1$(not node_(node_fail)) + p_gnn(grid, node, node_, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
358
359
360
361
        ) // END sum(gn2n_directional)

    =G=

362
    // Upward Demand for reserves due to potential interconnector failures
ran li's avatar
ran li committed
363
    + p_gnn(grid, node_fail, node, 'portion_of_transfer_to_reserve')
364
        * v_transferRightward(grid, node_fail, node, s, f, t)
365
    + p_gnn(grid, node, node_fail, 'portion_of_transfer_to_reserve')
ran li's avatar
ran li committed
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
        * v_transferLeftward(grid, node, node_fail, s, f, t)

    // Reserve provisions to another nodes via transfer links
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, 'up', node, node_)},
          // Reserve transfers to other nodes increase the reserve need of the present node
        + v_resTransferRightward(restype, 'up', node, node_, s, f+df_reserves(node, restype, f, t), t)
            * [1$(not node_(node_fail)) + p_gnn(grid, node, node_, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
        ) // END sum(gn2n_directional)
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, 'up', node, node_)},
          // Reserve transfers to other nodes increase the reserve need of the present node
        + v_resTransferLeftward(restype, 'up', node_, node, s, f+df_reserves(node, restype, f, t), t)
            * [1$(not node_(node_fail)) + p_gnn(grid, node_, node, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
        ) // END sum(gn2n_directional)

    // Reserve demand feasibility dummy variables
381
382
*    - vq_resDemand(restype, 'up', node, s, f+df_reserves(node, restype, f, t), t)
*    - vq_resMissing(restype, 'up', node, s, f+df_reserves(node, restype, f, t), t)${ft_reservesFixed(node, restype, f+df_reserves(node, restype, f, t), t)}
ran li's avatar
ran li committed
383
384
385
386
387
388
389
390
391
392
;

* --- N-1 Downward reserve demand due to a possibility that an interconnector that is transferring power to the node fails -------------------------------------------------
// NOTE! Currently, there are multiple identical instances of the reserve balance equation being generated for each forecast branch even when the reserves are committed and identical between the forecasts.
// NOTE! This could be solved by formulating a new "ft_reserves" set to cover only the relevant forecast-time steps, but it would possibly make the reserves even more confusing.

q_resDemandLargestInfeedTransfer2(grid, restypeDirectionNode(restype, 'down', node), node_fail, sft(s, f, t))
    ${  ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')
        and not [ restypeReleasedForRealization(restype)
                  and sft_realized(s, f, t)]
393
394
395
        and [ p_gnn(grid, node, node_fail, 'portion_of_transfer_to_reserve')
              or p_gnn(grid, node_fail, node, 'portion_of_transfer_to_reserve')
              ]
ran li's avatar
ran li committed
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
        and p_nReserves(node, restype, 'LossOfTrans')
        } ..

    // Reserve provision by capable units on this node
    + sum(nuft(node, unit, f, t)${nuRescapable(restype, 'down', node, unit)},
        + v_reserve(restype, 'down', node, unit, s, f+df_reserves(node, restype, f, t), t)
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                ] // END * v_reserve
        ) // END sum(nuft)

    // Reserve provision from other reserve categories when they can be shared
    + sum((nuft(node, unit, f, t), restype_)${p_nuRes2Res(node, unit, restype_, 'down', restype)},
        + v_reserve(restype_, 'down', node, unit, s, f+df_reserves(node, restype_, f, t), t)
            * p_nuRes2Res(node, unit, restype_, 'down', restype)
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                    * p_nuReserves(node, unit, restype_, 'reserveReliability')
                ] // END * v_reserve
        ) // END sum(nuft)

    // Reserve provision to this node via transfer links
    // SHOULD THE node_fail BE EXCLUDED?
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, 'down', node_, node)},
        + (1 - p_gnn(grid, node_, node, 'transferLoss') )
            * v_resTransferRightward(restype, 'down', node_, node, s, f+df_reserves(node_, restype, f, t), t)
            * [1$(not node_(node_fail)) + p_gnn(grid, node_, node, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
        ) // END sum(gn2n_directional)
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, 'down', node_, node)},
        + (1 - p_gnn(grid, node, node_, 'transferLoss') )
            * v_resTransferLeftward(restype, 'down', node, node_, s, f+df_reserves(node_, restype, f, t), t)
            * [1$(not node_(node_fail)) + p_gnn(grid, node, node_, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
        ) // END sum(gn2n_directional)

    =G=

    // Demand for reserves due to potential interconnector failures
    + p_gnn(grid, node_fail, node, 'portion_of_transfer_to_reserve')
ran li's avatar
ran li committed
436
        * v_transferLeftward(grid, node_fail, node, s, f, t)
437
    + p_gnn(grid, node, node_fail, 'portion_of_transfer_to_reserve')
ran li's avatar
ran li committed
438
        * v_transferRightward(grid, node, node_fail, s, f, t)
439
440

    // Reserve provisions to another nodes via transfer links
ran li's avatar
ran li committed
441
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, 'up', node, node_)},
442
          // Reserve transfers to other nodes increase the reserve need of the present node
ran li's avatar
ran li committed
443
        + v_resTransferRightward(restype, 'down', node, node_, s, f+df_reserves(node, restype, f, t), t)
444
445
            * [1$(not node_(node_fail)) + p_gnn(grid, node, node_, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
        ) // END sum(gn2n_directional)
ran li's avatar
ran li committed
446
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, 'down', node, node_)},
447
          // Reserve transfers to other nodes increase the reserve need of the present node
ran li's avatar
ran li committed
448
        + v_resTransferLeftward(restype, 'down', node_, node, s, f+df_reserves(node, restype, f, t), t)
449
450
451
452
            * [1$(not node_(node_fail)) + p_gnn(grid, node_, node, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
        ) // END sum(gn2n_directional)

    // Reserve demand feasibility dummy variables
453
454
*    - vq_resDemand(restype, 'down', node, s, f+df_reserves(node, restype, f, t), t)
*    - vq_resMissing(restype, 'down', node, s, f+df_reserves(node, restype, f, t), t)${ft_reservesFixed(node, restype, f+df_reserves(node, restype, f, t), t)}
455
;
456

ran li's avatar
ran li committed
457
458
459
460
461
462
463
464
* --- N-1 Upward reserve demand due to a possibility that an interconnector that is transferring power to the node fails -------------------------------------------------
// NOTE! Currently, there are multiple identical instances of the reserve balance equation being generated for each forecast branch even when the reserves are committed and identical between the forecasts.
// NOTE! This could be solved by formulating a new "ft_reserves" set to cover only the relevant forecast-time steps, but it would possibly make the reserves even more confusing.

*q_resDemandLargestInfeedTransfer(grid, restypeDirectionNode(restype, up_down, node), node_fail, sft(s, f, t))
*    ${  ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')
*        and not [ restypeReleasedForRealization(restype)
*                  and sft_realized(s, f, t)]
465
*        and (p_gnn(grid, node, node_fail, 'portion_of_transfer_to_reserve') and p_gnn(grid, node_fail, node, 'portion_of_transfer_to_reserve'))
ran li's avatar
ran li committed
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
*        and p_nReserves3D(node, restype, up_down, 'LossOfTrans')
*        } ..
*
*    // Reserve provision by capable units on this node
*    + sum(nuft(node, unit, f, t)${nuRescapable(restype, up_down, node, unit)},
*        + v_reserve(restype, up_down, node, unit, s, f+df_reserves(node, restype, f, t), t)
*            * [ // Account for reliability of reserves
*                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
*                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
*                ] // END * v_reserve
*        ) // END sum(nuft)
*
*    // Reserve provision from other reserve categories when they can be shared
*    + sum((nuft(node, unit, f, t), restype_)${p_nuRes2Res(node, unit, restype_, 'up', restype)},
*        + v_reserve(restype_, up_down, node, unit, s, f+df_reserves(node, restype_, f, t), t)
*            * p_nuRes2Res(node, unit, restype_, up_down, restype)
*            * [ // Account for reliability of reserves
*                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
*                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
*                    * p_nuReserves(node, unit, restype_, 'reserveReliability')
*                ] // END * v_reserve
*        ) // END sum(nuft)
*
*    // Reserve provision to this node via transfer links
*    // SHOULD THE node_fail BE EXCLUDED?
*    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, up_down, node_, node)},
*        + (1 - p_gnn(grid, node_, node, 'transferLoss') )
*            * v_resTransferRightward(restype, up_down, node_, node, s, f+df_reserves(node_, restype, f, t), t)
*            * [1$(not node_(node_fail)) + p_gnn(grid, node_, node, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
*        ) // END sum(gn2n_directional)
*    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, up_down, node_, node)},
*        + (1 - p_gnn(grid, node, node_, 'transferLoss') )
*            * v_resTransferLeftward(restype, up_down, node, node_, s, f+df_reserves(node_, restype, f, t), t)
*            * [1$(not node_(node_fail)) + p_gnn(grid, node, node_, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
*        ) // END sum(gn2n_directional)
*
*    =G=
*
*    // Upward Demand for reserves due to potential interconnector failures
*    [+ p_gnn(grid, node_fail, node, 'portion_of_transfer_to_reserve')
*        * v_transferRightward(grid, node_fail, node, s, f, t)
*    + p_gnn(grid, node, node_fail, 'portion_of_transfer_to_reserve')
*        * v_transferLeftward(grid, node, node_fail, s, f, t)]$(up_down eq 'up')
*    //Downward Demand for reserves due to potential interconnector failures
*    [+ p_gnn(grid, node_fail, node, 'portion_of_transfer_to_reserve')
*        * v_transferLeftward(grid, node_fail, node, s, f, t)
*    + p_gnn(grid, node, node_fail, 'portion_of_transfer_to_reserve')
*        * v_transferRightward(grid, node, node_fail, s, f, t)]$(up_down eq 'down')
*
*    // Reserve provisions to another nodes via transfer links
*    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, up_down, node, node_)},
*          // Reserve transfers to other nodes increase the reserve need of the present node
*        + v_resTransferRightward(restype, up_down, node, node_, s, f+df_reserves(node, restype, f, t), t)
*            * [1$(not node_(node_fail)) + p_gnn(grid, node, node_, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
*        ) // END sum(gn2n_directional)
*    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, up_down, node, node_)},
*          // Reserve transfers to other nodes increase the reserve need of the present node
*        + v_resTransferLeftward(restype, up_down, node_, node, s, f+df_reserves(node, restype, f, t), t)
*            * [1$(not node_(node_fail)) + p_gnn(grid, node_, node, 'portion_of_transfer_to_reserve')$(node_(node_fail))]
*        ) // END sum(gn2n_directional)
*
*    // Reserve demand feasibility dummy variables
*    - vq_resDemand(restype, up_down, node, s, f+df_reserves(node, restype, f, t), t)
*    - vq_resMissing(restype, up_down, node, s, f+df_reserves(node, restype, f, t), t)${ft_reservesFixed(node, restype, f+df_reserves(node, restype, f, t), t)}
*;

532
533
* --- Maximum Downward Capacity -----------------------------------------------

534
q_maxDownward(gnu(grid, node, unit), msft(m, s, f, t))
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
    ${  gnuft(grid, node, unit, f, t)
        and {
            [   ord(t) < tSolveFirst + smax(restype, p_nReserves(node, restype, 'reserve_length')) // Unit is either providing
                and sum(restype, nuRescapable(restype, 'down', node, unit)) // downward reserves
                ]
            // NOTE!!! Could be better to form a gnuft_reserves subset?
            or [ // the unit has an online variable
                uft_online(unit, f, t)
                and [
                    (unit_minLoad(unit) and p_gnu(grid, node, unit, 'unitSizeGen')) // generators with a min. load
                    or p_gnu(grid, node, unit, 'maxCons') // or consuming units with an online variable
                    ]
                ] // END or
            or [ // consuming units with investment possibility
                gnu_input(grid, node, unit)
                and [unit_investLP(unit) or unit_investMIP(unit)]
                ]
        }} ..

554
    // Energy generation/consumption
555
    + v_gen(grid, node, unit, s, f, t)
556
557

    // Considering output constraints (e.g. cV line)
558
559
    + sum(gngnu_constrainedOutputRatio(grid, node, grid_output, node_, unit),
        + p_gnu(grid_output, node_, unit, 'cV')
560
            * v_gen(grid_output, node_, unit, s, f, t)
561
562
563
        ) // END sum(gngnu_constrainedOutputRatio)

    // Downward reserve participation
564
    - sum(nuRescapable(restype, 'down', node, unit)${ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')},
565
        + v_reserve(restype, 'down', node, unit, s, f+df_reserves(node, restype, f, t), t) // (v_reserve can be used only if the unit is capable of providing a particular reserve)
566
567
568
569
570
571
        ) // END sum(nuRescapable)

    =G= // Must be greater than minimum load or maximum consumption  (units with min-load and both generation and consumption are not allowed)

    // Generation units, greater than minload
    + p_gnu(grid, node, unit, 'unitSizeGen')
572
        * sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
573
574
575
576
            + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
            + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
            ) // END sum(effGroup)
        * [ // Online variables should only be generated for units with restrictions
577
578
            + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f+df_central(f,t), t)} // LP online variant
            + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f+df_central(f,t), t)} // MIP online variant
579
580
            ] // END v_online

581
582
583
584
    // Units in run-up phase neet to keep up with the run-up rate
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(unitStarttype(unit, starttype)${uft_startupTrajectory(unit, f, t)},
            sum(runUpCounter(unit, counter)${t_active(t+dt_trajectory(counter))}, // Sum over the run-up intervals
585
586
                + [
                    + v_startup_LP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
587
                        ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
588
                    + v_startup_MIP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
589
                        ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
590
                    ]
591
                    * p_uCounter_runUpMin(unit, counter)
592
593
594
595
596
597
                ) // END sum(runUpCounter)
            ) // END sum(unitStarttype)

    // Units in shutdown phase need to keep up with the shutdown rate
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(shutdownCounter(unit, counter)${t_active(t+dt_trajectory(counter)) and uft_shutdownTrajectory(unit, f, t)}, // Sum over the shutdown intervals
598
599
600
601
602
603
            + [
                + v_shutdown_LP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                    ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                + v_shutdown_MIP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                    ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                ]
604
                * p_uCounter_shutdownMin(unit, counter)
605
            ) // END sum(shutdownCounter)
606

607
608
609
610
611
    // Consuming units, greater than maxCons
    // Available capacity restrictions
    - p_unit(unit, 'availability')
        * [
            // Capacity factors for flow units
612
            + sum(flowUnit(flow, unit),
613
                + ts_cf_(flow, node, f, t, s)
614
615
616
617
618
619
                ) // END sum(flow)
            + 1${not unit_flow(unit)}
            ] // END * p_unit(availability)
        * [
            // Online capacity restriction
            + p_gnu(grid, node, unit, 'maxCons')${not uft_online(unit, f, t)} // Use initial maximum if no online variables
620
621
622
623
            // !!! TEMPORARY SOLUTION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            + [
                + p_gnu(grid, node, unit, 'unitSizeCons')
                + p_gnu(grid, node, unit, 'maxCons')${not p_gnu(grid, node, unit, 'unitSizeCons') > 0}
624
                    / ( p_unit(unit, 'unitCount') + 1${not p_unit(unit, 'unitCount') > 0} )
625
626
                ]
            // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
627
                * [
628
                    // Capacity online
629
630
                    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
                    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}
631
632
633
634
635
636
637
638

                    // Investments to additional non-online capacity
                    + sum(t_invest(t_)${    ord(t_)<=ord(t)
                                            and not uft_online(unit, f, t)
                                            },
                        + v_invest_LP(unit, t_)${unit_investLP(unit)} // NOTE! v_invest_LP also for consuming units is positive
                        + v_invest_MIP(unit, t_)${unit_investMIP(unit)} // NOTE! v_invest_MIP also for consuming units is positive
                        ) // END sum(t_invest)
639
640
                    ] // END * p_gnu(unitSizeCons)
            ] // END * p_unit(availability)
641
;
642
643
644

* --- Maximum Upwards Capacity ------------------------------------------------

645
q_maxUpward(gnu(grid, node, unit), msft(m, s, f, t))
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
    ${  gnuft(grid, node, unit, f, t)
        and {
            [   ord(t) < tSolveFirst + smax(restype, p_nReserves(node, restype, 'reserve_length')) // Unit is either providing
                and sum(restype, nuRescapable(restype, 'up', node, unit)) // upward reserves
                ]
            or [
                uft_online(unit, f, t) // or the unit has an online variable
                and [
                    [unit_minLoad(unit) and p_gnu(grid, node, unit, 'unitSizeCons')] // consuming units with min_load
                    or [p_gnu(grid, node, unit, 'maxGen')]                          // generators with an online variable
                    ]
                ]
            or [
                gnu_output(grid, node, unit) // generators with investment possibility
                and (unit_investLP(unit) or unit_investMIP(unit))
                ]
        }}..

664
    // Energy generation/consumption
665
    + v_gen(grid, node, unit, s, f, t)
666
667
668
669

    // Considering output constraints (e.g. cV line)
    + sum(gngnu_constrainedOutputRatio(grid, node, grid_output, node_, unit),
        + p_gnu(grid_output, node_, unit, 'cV')
670
            * v_gen(grid_output, node_, unit, s, f, t)
671
672
673
        ) // END sum(gngnu_constrainedOutputRatio)

    // Upwards reserve participation
674
    + sum(nuRescapable(restype, 'up', node, unit)${ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')},
675
        + v_reserve(restype, 'up', node, unit, s, f+df_reserves(node, restype, f, t), t)
676
677
678
679
680
        ) // END sum(nuRescapable)

    =L= // must be less than available/online capacity

    // Consuming units
681
    - p_gnu(grid, node, unit, 'unitSizeCons')
682
        * sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
683
684
685
686
            + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
            + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
            ) // END sum(effGroup)
        * [
687
688
            + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)} // Consuming units are restricted by their min. load (consuming is negative)
            + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)} // Consuming units are restricted by their min. load (consuming is negative)
689
690
691
692
693
694
695
            ] // END * p_gnu(unitSizeCons)

    // Generation units
    // Available capacity restrictions
    + p_unit(unit, 'availability') // Generation units are restricted by their (available) capacity
        * [
            // Capacity factor for flow units
696
            + sum(flowUnit(flow, unit),
697
                + ts_cf_(flow, node, f, t, s)
698
699
700
701
702
703
704
705
                ) // END sum(flow)
            + 1${not unit_flow(unit)}
            ] // END * p_unit(availability)
        * [
            // Online capacity restriction
            + p_gnu(grid, node, unit, 'maxGen')${not uft_online(unit, f, t)} // Use initial maxGen if no online variables
            + p_gnu(grid, node, unit, 'unitSizeGen')
                * [
706
                    // Capacity online
707
708
                    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f ,t)}
                    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}
709
710
711
712
713
714
715
716

                    // Investments to non-online capacity
                    + sum(t_invest(t_)${    ord(t_)<=ord(t)
                                            and not uft_online(unit, f ,t)
                                            },
                        + v_invest_LP(unit, t_)${unit_investLP(unit)}
                        + v_invest_MIP(unit, t_)${unit_investMIP(unit)}
                        ) // END sum(t_invest)
717
718
                    ] // END * p_gnu(unitSizeGen)
            ] // END * p_unit(availability)
719

720
721
722
723
    // Units in run-up phase neet to keep up with the run-up rate
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(unitStarttype(unit, starttype)${uft_startupTrajectory(unit, f, t)},
            sum(runUpCounter(unit, counter)${t_active(t+dt_trajectory(counter))}, // Sum over the run-up intervals
724
725
                + [
                    + v_startup_LP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
726
                        ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
727
                    + v_startup_MIP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
728
                        ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
729
                    ]
730
                    * p_uCounter_runUpMax(unit, counter)
731
732
733
734
735
736
                ) // END sum(runUpCounter)
            ) // END sum(unitStarttype)

    // Units in shutdown phase need to keep up with the shutdown rate
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(shutdownCounter(unit, counter)${t_active(t+dt_trajectory(counter)) and uft_shutdownTrajectory(unit, f, t)}, // Sum over the shutdown intervals
737
738
739
740
741
742
            + [
                + v_shutdown_LP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                    ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                + v_shutdown_MIP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                    ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                ]
743
                * p_uCounter_shutdownMax(unit, counter)
744
            ) // END sum(shutdownCounter)
745
;
746

747
748
* --- Reserve Provision of Units with Investments -----------------------------

749
750
751
752
q_reserveProvision(nuRescapable(restypeDirectionNode(restype, up_down, node), unit), sft(s, f, t))
    ${  ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')
        and nuft(node, unit, f, t)
        and (unit_investLP(unit) or unit_investMIP(unit))
753
754
        and not sum(restypeDirectionGridNodeGroup(restype, up_down, grid, node, group),
                    ft_reservesFixed(group, restype, f+df_reservesGroup(group, restype, f, t), t))
755
756
        } ..

757
    + v_reserve(restype, up_down, node, unit, s, f+df_reserves(node, restype, f, t), t)
758
759
760
761
762
763
764
765
766
767
768
769
770

    =L=

    + p_nuReserves(node, unit, restype, up_down)
        * [
            + sum(grid, p_gnu(grid, node, unit, 'maxGen') + p_gnu(grid, node, unit, 'maxCons') )  // Reserve sets and variables are currently lacking the grid dimension...
            + sum(t_invest(t_)${ ord(t_)<=ord(t) },
                + v_invest_LP(unit, t_)${unit_investLP(unit)}
                    * sum(grid, p_gnu(grid, node, unit, 'unitSizeTot')) // Reserve sets and variables are currently lacking the grid dimension...
                + v_invest_MIP(unit, t_)${unit_investMIP(unit)}
                    * sum(grid, p_gnu(grid, node, unit, 'unitSizeTot')) // Reserve sets and variables are currently lacking the grid dimension...
                ) // END sum(t_)
            ]
771
772
773
774
        * p_unit(unit, 'availability') // Taking into account availability...
        * [
            // ... and capacity factor for flow units
            + sum(flowUnit(flow, unit),
775
                + ts_cf_(flow, node, f, t, s)
776
777
                ) // END sum(flow)
            + 1${not unit_flow(unit)}
778
779
780
            ] // How to consider reserveReliability in the case of investments when we typically only have "realized" time steps?
;

781
782
* --- Unit Startup and Shutdown -----------------------------------------------

783
784
785
786
q_startshut(ms(m, s), uft_online(unit, f, t))
    ${  msft(m, s, f, t)
        }..

787
    // Units currently online
788
789
    + v_online_LP (unit, s, f+df_central(f,t), t)${uft_onlineLP (unit, f, t)}
    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}
790
791

    // Units previously online
792
    // The same units
793
    - v_online_LP (unit, s+ds(s,t), f+df(f,t+dt(t)), t+dt(t))${ uft_onlineLP_withPrevious(unit, f+df(f,t+dt(t)), t+dt(t))
794
                                                             and not uft_aggregator_first(unit, f, t) } // This reaches to tFirstSolve when dt = -1
795
    - v_online_MIP(unit, s+ds(s,t), f+df(f,t+dt(t)), t+dt(t))${ uft_onlineMIP_withPrevious(unit, f+df(f,t+dt(t)), t+dt(t))
796
797
798
799
                                                             and not uft_aggregator_first(unit, f, t) }

    // Aggregated units just before they are turned into aggregator units
    - sum(unit_${unitAggregator_unit(unit, unit_)},
800
801
        + v_online_LP (unit_, s, f+df(f,t+dt(t)), t+dt(t))${uft_onlineLP_withPrevious(unit_, f+df(f,t+dt(t)), t+dt(t))}
        + v_online_MIP(unit_, s, f+df(f,t+dt(t)), t+dt(t))${uft_onlineMIP_withPrevious(unit_, f+df(f,t+dt(t)), t+dt(t))}
802
        )${uft_aggregator_first(unit, f, t)} // END sum(unit_)
803

804
805
    =E=

806
    // Unit startup and shutdown
807

808
    // Add startup of units dt_toStartup before the current t (no start-ups for aggregator units before they become active)
809
    + sum(unitStarttype(unit, starttype),
810
        + v_startup_LP(unit, starttype, s, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t))
811
            ${ uft_onlineLP_withPrevious(unit, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) }
812
        + v_startup_MIP(unit, starttype, s, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t))
813
            ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) }
814
        )${not [unit_aggregator(unit) and ord(t) + dt_toStartup(unit, t) <= tSolveFirst + p_unit(unit, 'lastStepNotAggregated')]} // END sum(starttype)
815

816
817
818
819
    // NOTE! According to 3d_setVariableLimits,
    // cannot start a unit if the time when the unit would become online is outside
    // the horizon when the unit has an online variable
    // --> no need to add start-ups of aggregated units to aggregator units
820

821
    // Shutdown of units at time t
822
823
824
825
    - v_shutdown_LP(unit, s, f, t)
        ${ uft_onlineLP(unit, f, t) }
    - v_shutdown_MIP(unit, s, f, t)
        ${ uft_onlineMIP(unit, f, t) }
826
;
827

828
*--- Startup Type -------------------------------------------------------------
829
// !!! NOTE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
830
831
832
// This formulation doesn't work as intended when unitCount > 1, as one recent
// shutdown allows for multiple hot/warm startups on subsequent time steps.
// Pending changes.
833

834
835
836
837
q_startuptype(ms(m, s), starttypeConstrained(starttype), uft_online(unit, f, t))
    ${  msft(m, s, f, t)
        and unitStarttype(unit, starttype)
        } ..
838
839

    // Startup type
840
841
    + v_startup_LP(unit, starttype, s, f, t)${ uft_onlineLP(unit, f, t) }
    + v_startup_MIP(unit, starttype, s, f, t)${ uft_onlineMIP(unit, f, t) }
842
843
844
845

    =L=

    // Subunit shutdowns within special startup timeframe
846
847
848
849
850
851
852
    + sum(unitCounter(unit, counter)${  dt_starttypeUnitCounter(starttype, unit, counter)
                                        and t_active(t+(dt_starttypeUnitCounter(starttype, unit, counter)+1))
                                        },
        + v_shutdown_LP(unit, s, f+df(f,t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)), t+(dt_starttypeUnitCounter(starttype, unit, counter)+1))
            ${ uft_onlineLP_withPrevious(unit, f+df(f,t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)), t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)) }
        + v_shutdown_MIP(unit, s, f+df(f,t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)), t+(dt_starttypeUnitCounter(starttype, unit, counter)+1))
            ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)), t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)) }
853
854
855
        ) // END sum(counter)

    // NOTE: for aggregator units, shutdowns for aggregated units are not considered
856
;
857

858

859
860
*--- Online Limits with Startup Type Constraints and Investments --------------

861
862
863
864
865
866
867
868
869
q_onlineLimit(ms(m, s), uft_online(unit, f, t))
    ${  msft(m, s, f, t)
        and {
            p_unit(unit, 'minShutdownHours')
            or p_u_runUpTimeIntervals(unit)
            or unit_investLP(unit)
            or unit_investMIP(unit)
        }} ..

870
    // Online variables
871
872
    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f ,t)}
873
874
875
876
877
878

    =L=

    // Number of existing units
    + p_unit(unit, 'unitCount')

879
    // Number of units unable to become online due to restrictions
880
881
882
883
884
885
886
    - sum(unitCounter(unit, counter)${  dt_downtimeUnitCounter(unit, counter)
                                        and t_active(t+(dt_downtimeUnitCounter(unit, counter) + 1))
                                        },
        + v_shutdown_LP(unit, s, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1))
            ${ uft_onlineLP_withPrevious(unit, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1)) }
        + v_shutdown_MIP(unit, s, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1))
            ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1)) }
887
888
889
        ) // END sum(counter)

    // Number of units unable to become online due to restrictions (aggregated units in the past horizon or if they have an online variable)
890
    - sum(unitAggregator_unit(unit, unit_),
891
892
893
894
895
896
897
        + sum(unitCounter(unit, counter)${  dt_downtimeUnitCounter(unit, counter)
                                            and t_active(t+(dt_downtimeUnitCounter(unit, counter) + 1))
                                            },
            + v_shutdown_LP(unit_, s, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1))
                ${ uft_onlineLP_withPrevious(unit_, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1)) }
            + v_shutdown_MIP(unit_, s, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1))
                ${ uft_onlineMIP_withPrevious(unit_, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1)) }
898
899
            ) // END sum(counter)
        )${unit_aggregator(unit)} // END sum(unit_)
900
901
902

    // Investments into units
    + sum(t_invest(t_)${ord(t_)<=ord(t)},
903
904
        + v_invest_LP(unit, t_)${unit_investLP(unit)}
        + v_invest_MIP(unit, t_)${unit_investMIP(unit)}
905
906
907
        ) // END sum(t_invest)
;

908
909
910
911
*--- Both q_offlineAfterShutdown and q_onlineOnStartup work when there is only one unit.
*    These equations prohibit single units turning on and off at the same time step.
*    Unfortunately there seems to be no way to prohibit this when unit count is > 1.
*    (it shouldn't be worthwhile anyway if there is a startup cost, but it can fall within the solution gap).
912
913
914
915
q_onlineOnStartUp(s_active(s), uft_online(unit, f, t))
    ${  sft(s, f, t)
        and sum(starttype, unitStarttype(unit, starttype))
        }..
916
917

    // Units currently online
918
919
    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}
920
921
922
923

    =G=

    + sum(unitStarttype(unit, starttype),
924
        + v_startup_LP(unit, starttype, s, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) //dt_toStartup displaces the time step to the one where the unit would be started up in order to reach online at t
925
            ${ uft_onlineLP_withPrevious(unit, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) }
926
        + v_startup_MIP(unit, starttype, s, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) //dt_toStartup displaces the time step to the one where the unit would be started up in order to reach online at t
927
            ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) }
928
929
930
      ) // END sum(starttype)
;

931
932
933
934
q_offlineAfterShutdown(s_active(s), uft_online(unit, f, t))
    ${  sft(s, f, t)
        and sum(starttype, unitStarttype(unit, starttype))
        }..
935

936
937
938
939
940
    // Number of existing units
    + p_unit(unit, 'unitCount')

    // Investments into units
    + sum(t_invest(t_)${ord(t_)<=ord(t)},
941
942
        + v_invest_LP(unit, t_)${unit_investLP(unit)}
        + v_invest_MIP(unit, t_)${unit_investMIP(unit)}
943
944
        ) // END sum(t_invest)

945
    // Units currently online
946
947
    - v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
    - v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}
948
949
950

    =G=

951
952
953
954
    + v_shutdown_LP(unit, s, f, t)
        ${ uft_onlineLP(unit, f, t) }
    + v_shutdown_MIP(unit, s, f, t)
        ${ uft_onlineMIP(unit, f, t) }
955
956
;

957
958
*--- Minimum Unit Uptime ------------------------------------------------------

959
960
961
962
q_onlineMinUptime(ms(m, s), uft_online(unit, f, t))
    ${  msft(m, s, f, t)
        and  p_unit(unit, 'minOperationHours')
        } ..
963
964

    // Units currently online
965
966
    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}
967
968
969
970

    =G=

    // Units that have minimum operation time requirements active
971
972
973
    + sum(unitCounter(unit, counter)${  dt_uptimeUnitCounter(unit, counter)
                                        and t_active(t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) // Don't sum over counters that don't point to an active time step
                                        },
974
        + sum(unitStarttype(unit, starttype),
975
            + v_startup_LP(unit, starttype, s, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1))
976
                ${ uft_onlineLP_withPrevious(unit, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) }
977
            + v_startup_MIP(unit, starttype, s, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1))
978
                ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) }
979
            ) // END sum(starttype)
980
981
982
        ) // END sum(counter)

    // Units that have minimum operation time requirements active (aggregated units in the past horizon or if they have an online variable)
Topi Rasku's avatar
Topi Rasku committed
983
    + sum(unitAggregator_unit(unit, unit_),
984
985
986
        + sum(unitCounter(unit, counter)${  dt_uptimeUnitCounter(unit, counter)
                                            and t_active(t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) // Don't sum over counters that don't point to an active time step
                                            },
987
            + sum(unitStarttype(unit, starttype),
988
                + v_startup_LP(unit, starttype, s, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1))
989
                    ${ uft_onlineLP_withPrevious(unit, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) }
990
                + v_startup_MIP(unit, starttype, s, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1))
991
                    ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) }
992
993
994
                ) // END sum(starttype)
            ) // END sum(counter)
        )${unit_aggregator(unit)} // END sum(unit_)
995
996
;

997
998
* --- Cyclic Boundary Conditions for Online State -----------------------------

999
1000
q_onlineCyclic(uss_bound(unit, s_, s), m)
    ${  ms(m, s_)