2b_equations.gms 41.2 KB
Newer Older
1
2
equations
    q_obj "Objective function"
3
    q_balance(grid, node, mType, f, t) "Energy demand must be satisfied at each node"
4
    q_resDemand(restype, up_down, node, f, t) "Procurement for each reserve type is greater than demand"
5
    q_resTransfer(grid, node, node, f, t) "Transfer of energy and capacity reservations are less than the transfer capacity"
6
    q_maxDownward(grid, node, unit, f, t) "Downward commitments will not undercut power plant minimum load constraints or maximum elec. consumption"
7
    q_maxUpward(grid, node, unit, f, t) "Upward commitments will not exceed maximum available capacity or consumed power"
Juha Kiviluoma's avatar
Juha Kiviluoma committed
8
    q_startup(unit, f, t) "Capacity started up is greater than the difference of online cap. now and in the previous time step"
9
10
    q_genRamp(grid, node, mType, unit, f, t) "Record the ramps of units with ramp restricitions or costs"
    q_genRampChange(grid, node, mType, unit, f, t) "Record the ramp rates of units with ramping costs"
11
    q_conversionDirectInputOutput(effSelector, unit, f, t) "Direct conversion of inputs to outputs (no piece-wise linear part-load efficiencies)"
Juha Kiviluoma's avatar
Juha Kiviluoma committed
12
13
14
15
16
    q_conversionSOS2InputIntermediate(effSelector, unit, f, t)   "Intermediate output when using SOS2 variable based part-load piece-wise linearization"
    q_conversionSOS2Constraint(effSelector, unit, f, t)          "Sum of v_sos2 has to equal v_online"
    q_conversionSOS2IntermediateOutput(effSelector, unit, f, t)  "Output is forced equal with v_sos2 output"
    q_outputRatioFixed(grid, node, grid, node, unit, f, t) "Force fixed ratio between two energy outputs into different energy grids"
    q_outputRatioConstrained(grid, node, grid, node, unit, f, t) "Constrained ratio between two grids of energy output; e.g. electricity generation is greater than cV times unit_heat generation in extraction plants"
17
    q_stateSlack(grid, node, slack, f, t) "Slack variable greater than the difference between v_state and the slack boundary"
18
19
    q_stateUpwardLimit(grid, node, mType, f, t) "Limit the commitments of a node with a state variable to the available headrooms"
    q_stateDownwardLimit(grid, node, mType, f, t) "Limit the commitments of a node with a state variable to the available headrooms"
20
    q_boundState(grid, node, node, mType, f, t) "Node state variables bounded by other nodes"
21
    q_boundStateMaxDiff(grid, node, node, mType, f, t) "Node state variables bounded by other nodes (maximum state difference)"
22
    q_boundCyclic(grid, node, mType, f, t, t_) "Cyclic bound for the first and the last state"
23
    q_bidirectionalTransfer(grid, node, node, f, t) "Possible common transfer capacity constraint for interconnected transfer variables"
24
25
26
;


27
$setlocal def_penalty 1e9
28
29
30
31
Scalars
    PENALTY "Default equation violation penalty" / %def_penalty% /
;
Parameters
32
    PENALTY_BALANCE(grid) "Penalty on violating energy balance eq. (/MWh)"
33
    PENALTY_RES(restype, up_down) "Penalty on violating a reserve (/MW)"
34
;
35
PENALTY_BALANCE(grid) = %def_penalty%;
36
PENALTY_RES(restype, up_down) =  1e-3*%def_penalty%;
37
38
39

* -----------------------------------------------------------------------------
q_obj ..
Juha Kiviluoma's avatar
Juha Kiviluoma committed
40
  + v_obj * 1000000
41
42
  =E=
  + sum(msft(m, s, f, t),
43
        p_sft_Probability(s,f,t) *
44
        (
45
46
         // Variable O&M costs
         + sum(gnu_output(grid, node, unit),  // Calculated only for output energy
Juha Kiviluoma's avatar
Juha Kiviluoma committed
47
                p_unit(unit, 'omCosts') *
48
49
                $$ifi not '%rampSched%' == 'yes' p_stepLength(m, f, t) *
                $$ifi '%rampSched%' == 'yes' (p_stepLength(m, f, t) + p_stepLength(m, f, t+1))/2 *
50
                     v_gen(grid, node, unit, f, t)$nuft(node, unit, f, t)
51
           )
52
         // Fuel and emission costs
Juha Kiviluoma's avatar
Juha Kiviluoma committed
53
54
55
56
57
58
59
60
61
62
63
         + sum((node, unit_fuel, fuel)$(nu(node, unit_fuel) and uFuel(unit_fuel, 'main', fuel)),
              + p_stepLength(m, f, t)
              * v_fuelUse(fuel, unit_fuel, f, t)$nuft(node, unit_fuel, f, t)
                  * ( + sum{tFuel$[ord(tFuel) <= ord(t)],
                            ts_fuelPriceChangenode(fuel, node, tFuel) }  // Fuel costs, sum initial fuel price plus all subsequent changes to the fuelprice
                      + sum{emission,         // Emission taxes
                            p_fuelEmission(fuel, emission) / 1e3
                              * sum(grid$gnu_output(grid, node, unit_fuel), p_gnPolicy(grid, node, 'emissionTax', emission))  // Sum emission costs from different output energy types
                        }
                     )
           )
64
         // Start-up costs
65
         + sum(uft_online(unit, f, t),
66
             + {
67
68
                 + v_startup(unit, f, t) // Cost of starting up
               } / p_unit(unit, 'unitCount')
69
70
             * {
                  // Startup variable costs
71
                 + p_unit(unit, 'startCost')
72
                 * p_unit(unit, 'outputCapacityTotal')
73
                  // Start-up fuel and emission costs
Juha Kiviluoma's avatar
Juha Kiviluoma committed
74
                 + sum(uFuel(unit_fuel, 'startup', fuel),
75
                     + p_unit(unit, 'startFuelCons')
76
                     * p_unit(unit, 'outputCapacityTotal')
77
                     * sum(gnu_output(grid, node, unit),
78
                           // Fuel costs for start-up fuel use
Juha Kiviluoma's avatar
Juha Kiviluoma committed
79
80
81
82
83
84
85
                         + ( + sum{tFuel$[ord(tFuel) <= ord(t)],
                                   ts_fuelPriceChangenode(fuel, node, tFuel) }
                               // Emission taxes of startup fuel use
                             + sum(emission,
                                p_fuelEmission(fuel, emission) / 1e3
                                  * p_gnPolicy(grid, node, 'emissionTax', emission)  // Sum emission costs from different output energy types
                               )
86
87
                           ) / p_gnu(grid, node, unit, 'maxGen')  // Calculate these in relation to maximum output ratios between multiple outputs
                       ) * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'maxGen'))  // see line above
88
89
90
                   )
               }
           )
91
92
93
94
95
96
97
98
         // Ramping costs
         + sum(gnuft_ramp(grid, node, unit, f, t)${ p_gnu(grid, node, unit, 'rampUpCost') OR p_gnu(grid, node, unit, 'rampDownCost') },
            + (p_gnu(grid, node, unit, 'maxGen') + p_gnu(grid, node, unit, 'maxCons')) // NOTE! Doens't work correctly if a gnu has both! Is that even possible, though?
            * ( // Changes in ramp rates
                + p_gnu(grid, node, unit, 'rampUpCost') * v_genRampChange(grid, node, unit, 'up', f, t)
                + p_gnu(grid, node, unit, 'rampDownCost') * v_genRampChange(grid, node, unit, 'down', f, t)
              )
           )
99
        )  // END * p_sft_probability(s,f,t)
100
    ) // END sum over msft(m, s, f, t)
101
    // Value of energy storage change
102
103
  - sum(mftLastSteps(m, f, t)$active('storageValue'),
         sum(gn_state(grid, node),
Juha Kiviluoma's avatar
Juha Kiviluoma committed
104
              p_storageValue(grid, node, t) *
105
106
107
108
109
110
111
112
                  sum(s$p_sft_probability(s,f,t), p_sft_probability(s, f, t) * v_state(grid, node, f, t))
         )
    )
  + sum(mftStart(m, f, t)$active('storageValue'),
         sum(gn_state(grid, node),
              p_storageValue(grid, node, t) *
                  sum(s$p_sft_probability(s,f,t), p_sft_probability(s, f, t) * v_state(grid, node, f, t))
         )
Juha Kiviluoma's avatar
Juha Kiviluoma committed
113
    )
114
115
116
  - sum([s, m, uft_online(unit, f, t)]$mftStart(m, f, t), p_sft_probability(s, f, t) * 0.5 * v_online(unit, f, t))     // minus value of avoiding startup costs before
  - sum((s, uft_online(unit, f, t))$(ord(t) = p_uft_online_last(unit, f, t)), p_sft_probability(s, f, t) * 0.5 * v_online(unit, f, t)) // or after the model solve

Juha Kiviluoma's avatar
Juha Kiviluoma committed
117
    // Dummy variables
118
  + sum(msft(m, s, f, t), p_sft_probability(s, f, t) * (
Juha Kiviluoma's avatar
Juha Kiviluoma committed
119
        sum(inc_dec,
120
            sum( gn(grid, node), vq_gen(inc_dec, grid, node, f, t) * (p_stepLength(m, f, t) + p_stepLength(m, f+pf(f,t), t+pt(t))${not p_stepLength(m, f, t)}) * PENALTY_BALANCE(grid) )
121
        )
122
123
        + sum(restypeDirectionNode(restype, up_down, node),
              vq_resDemand(restype, up_down, node, f, t)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
124
            * p_stepLength(m, f, t)
125
            * PENALTY_RES(restype, up_down)
126
127
          )
      )
Juha Kiviluoma's avatar
Juha Kiviluoma committed
128
129
    )
    // Node state slack variable penalties
130
  + sum(gn_stateSlack(grid, node),
131
      + sum(msft(m, s, f, t),
132
          + sum(upwardSlack$p_gnBoundaryPropertiesForStates(grid, node, upwardSlack, 'slackCost'),
133
              + p_sft_probability(s, f, t) * ( p_gnBoundaryPropertiesForStates(grid, node,   upwardSlack, 'slackCost') * v_stateSlack(grid, node,   upwardSlack, f, t) )
134
            )
135
          + sum(downwardSlack$p_gnBoundaryPropertiesForStates(grid, node, downwardSlack, 'slackCost'),
136
              + p_sft_probability(s, f, t) * ( p_gnBoundaryPropertiesForStates(grid, node, downwardSlack, 'slackCost') * v_stateSlack(grid, node, downwardSlack, f, t) )
137
            )
Juha Kiviluoma's avatar
Juha Kiviluoma committed
138
139
        )
    )
140
;
141

142
* -----------------------------------------------------------------------------
143
q_balance(gn(grid, node), m, ft_dynamic(f, t))${p_stepLength(m, f+pf(f,t), t+pt(t)) and not p_gn(grid, node, 'boundAll')} .. // Energy/power balance dynamics solved using implicit Euler discretization
144
  // The left side of the equation is the change in the state (will be zero if the node doesn't have a state)
145
  + 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)
146
      * ( + v_state(grid, node, f+cf(f,t), t)                           // The difference between current
147
          - v_state(grid, node, f+pf(f,t), t+pt(t))                     // ... and previous state of the node
148
149
150
151
152
153
        )
  =E=
  // The right side of the equation contains all the changes converted to energy terms
  + (
      + (
          // Self discharge out of the model boundaries
154
          - p_gn(grid, node, 'selfDischargeLoss')$gn_state(grid, node) * (
155
              + v_state(grid, node, f+cf(f,t), t)                                               // The current state of the node
156
              $$ifi '%rampSched%' == 'yes' + v_state(grid, node, f+pf(f,t), t+pt(t))    // and possibly averaging with the previous state of the node
157
158
159
160
            )
          // Energy diffusion from this node to neighbouring nodes
          - sum(to_node$(gnn_state(grid, node, to_node)),
              + p_gnn(grid, node, to_node, 'diffCoeff') * (
161
                  + v_state(grid, node, f+cf(f,t), t)
162
163
164
165
166
167
                  $$ifi '%rampSched%' == 'yes' + v_state(grid, node, f+pf(f,t), t+pt(t))
                )
            )
          // Energy diffusion from neighbouring nodes to this node
          + sum(from_node$(gnn_state(grid, from_node, node)),
              + p_gnn(grid, from_node, node, 'diffCoeff') * (
168
                  + v_state(grid, from_node, f+cf(f,t), t)                                             // Incoming diffusion based on the state of the neighbouring node
169
170
171
172
173
174
175
                  $$ifi '%rampSched%' == 'yes' + v_state(grid, from_node, f+pf(f,t), t+pt(t))  // Ramp schedule averaging, NOTE! State and other terms use different indeces for non-ramp-schedule!
                )
            )
          // Controlled energy transfer from other nodes to this one
          + sum(from_node$(gn2n(grid, from_node, node)),
              + (1 - p_gnn(grid, from_node, node, 'transferLoss')) * (    // Include transfer losses
                  + v_transfer(grid, from_node, node, f+pf(f,t), t+pt(t))
176
                  $$ifi '%rampSched%' == 'yes' + v_transfer(grid, from_node, node, f+cf(f,t), t)    // Ramp schedule averaging, NOTE! State and other terms use different indeces for non-ramp-schedule!
177
178
179
180
181
                )
            )
          // Controlled energy transfer to other nodes from this one
          - sum(to_node$(gn2n(grid, node, to_node)),
              + v_transfer(grid, node, to_node, f+pf(f,t), t+pt(t))   // Transfer losses accounted for in the previous term
182
              $$ifi '%rampSched%' == 'yes' + v_transfer(grid, node, to_node, f+cf(f,t), t)    // Ramp schedule averaging
183
184
            )
          // Interactions between the node and its units
185
186
          + sum(gnuft(grid, node, unit, f+pf(f,t), t+pt(t)),
              + v_gen(grid, node, unit, f+pf(f,t), t+pt(t)) // Unit energy generation and consumption
187
              $$ifi '%rampSched%' == 'yes' + v_gen(grid, node, unit, f+cf(f,t), t)
188
189
190
            )
          // Spilling energy out of the endogenous grids in the model
          - v_spill(grid, node, f+pf(f,t), t+pt(t))$node_spill(node)
191
          $$ifi '%rampSched%' == 'yes' - v_spill(grid, node, f+cf(f,t), t)$node_spill(node)
192
          // Power inflow and outflow timeseries to/from the node
193
          + ts_influx_(grid, node, f+pf(f,t), t+pt(t))   // Incoming (positive) and outgoing (negative) absolute value time series
194
          $$ifi '%rampSched%' == 'yes' + ts_influx_(grid, node, f+cf(f,t), t)
195
196
          // Dummy generation variables, for feasibility purposes
          + vq_gen('increase', grid, node, f+pf(f,t), t+pt(t)) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
197
          $$ifi '%rampSched%' == 'yes' + vq_gen('increase', grid, node, f+cf(f,t), t)
198
          - vq_gen('decrease', grid, node, f+pf(f,t), t+pt(t)) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
199
          $$ifi '%rampSched%' == 'yes' - vq_gen('decrease', grid, node, f+cf(f,t), t)
200
201
202
        ) * p_stepLength(m, f+pf(f,t), t+pt(t))   // Multiply by time step to get energy terms
    )
  $$ifi '%rampSched%' == 'yes' / 2    // Averaging all the terms on the right side of the equation over the timestep here.
203
204
;
* -----------------------------------------------------------------------------
205
q_resDemand(restypeDirectionNode(restype, up_down, node), ft(f, t))${ord(t) le tSolveFirst + sum[mf(m, f), mSettings(m, 't_reserveLength')] and [not (restype('tertiary') and ord(t) le tSolveFirst + tDispatchCurrent)]} ..
206
  + sum(nuft(node, unit, f, t)$nuRescapable(restype, up_down, node, unit),   // Reserve capable units on this node
207
        v_reserve(restype, up_down, node, unit, f, t) * p_nuReserves(node, unit, restype, 'reserveContribution')
208
    )
209
  + sum(gnu_input(grid, node, unit)${gnuft(grid, node, unit, f, t) AND nuRescapable(restype, up_down, node, unit)},
210
        v_reserve(restype, up_down, node, unit, f, t) * p_nuReserves(node, unit, restype, 'reserveContribution') // Reserve capable units with input from this node
211
    )
212
  + sum(gn2n(grid, from_node, node)$restypeDirectionNode(restype, up_down, from_node),
213
        (1 - p_gnn(grid, from_node, node, 'transferLoss')
214
        ) * v_resTransfer(restype, up_down, from_node, node, f, t)             // Reserves from another node - reduces the need for reserves in the node
215
216
    )
  =G=
217
218
  + ts_reserveDemand_(restype, up_down, node, f, t)$p_nReserves(node, restype, 'use_time_series')
  + p_nReserves(node, restype, up_down)${not p_nReserves(node, restype, 'use_time_series')}
219
220
221
  - vq_resDemand(restype, up_down, node, f, t)
  + sum(gn2n(grid, node, to_node)$restypeDirectionNode(restype, up_down, to_node),   // If trasferring reserves to another node, increase your own reserves by same amount
        v_resTransfer(restype, up_down, node, to_node, f, t)
222
223
224
    )
;
* -----------------------------------------------------------------------------
225
226
q_resTransfer(gn2n(grid, from_node, to_node), ft(f, t))${ sum(restypeDirection(restype, up_down), restypeDirectionNode(restype, up_down, from_node))
                                                            OR sum(restypeDirection(restype, up_down), restypeDirectionNode(restype, up_down, to_node))
227
228
                                                            } ..
  + v_transfer(grid, from_node, to_node, f, t)
229
230
  + sum(restypeDirection(restype, up_down)$(restypeDirectionNode(restype, up_down, from_node) and restypeDirectionNode(restype, up_down, to_node)),
        + v_resTransfer(restype, up_down, from_node, to_node, f, t)
231
232
233
234
235
    )
  =L=
  + p_gnn(grid, from_node, to_node, 'transferCap')
;
* -----------------------------------------------------------------------------
236
237
238
239
240
241
242
243
244
245
q_maxDownward(gnuft(grid, node, unit, f, t))${ [     ord(t) < tSolveFirst + settings('t_reserveLength')               // Unit is either providing
                                               and sum(restype, nuRescapable(restype, 'down', node, unit))            // downward reserves
                                             ] or
                                             [ uft_online(unit, f, t) and                                             // or the unit has an online varaible
                                                 (
                                                      [unit_minLoad(unit) and p_gnu(grid, node, unit, 'maxGen')]      // generators with a min. load
                                                   or [p_gnu(grid, node, unit, 'maxCons')]                        // or consuming units with an online variable
                                                 )
                                             ]
                                           }..
246
  + v_gen(grid, node, unit, f, t)                                                                                    // energy generation/consumption
Juha Kiviluoma's avatar
Juha Kiviluoma committed
247
248
  + sum( gngnu_constrainedOutputRatio(grid, node, grid_, node_, unit),
        p_gnu(grid_, node_, unit, 'cV') * v_gen(grid_, node_, unit, f, t) )                              // considering output constraints (e.g. cV line)
Topi Rasku's avatar
Topi Rasku committed
249
  - sum(nuRescapable(restype, 'down', node, unit)$[unit_elec(unit) and ord(t) < tSolveFirst + settings('t_reserveLength')],               // minus downward reserve participation
250
        v_reserve(restype, 'down', node, unit, f, t)                                                              // (v_reserve can be used only if the unit is capable of providing a particular reserve)
251
    )
Juha Kiviluoma's avatar
Juha Kiviluoma committed
252
  =G=                                                                                                                // must be greater than minimum load or maximum consumption  (units with min-load and both generation and consumption are not allowed)
253
  + v_online(unit, f, t)${uft_online(unit, f, t) and p_gnu(grid, node, unit, 'maxGen')} // Online variables should only be generated for units with restrictions
254
    / p_unit(unit, 'unitCount')
255
256
257
258
    * p_gnu(grid, node, unit, 'maxGen')
    * sum(effGroup, // Uses the minimum 'lb' for the current efficiency approximation
        + (p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)} + ts_effGroupUnit(effGroup, unit, 'lb', f, t))
      )
259
260
261
262
  + v_gen.lo(grid, node, unit, f, t)$p_gnu(grid, node, unit, 'maxCons')           // notice: v_gen.lo for consuming units is negative
  - v_online(unit, f, t)${uft_online(unit, f, t) and p_gnu(grid, node, unit, 'maxCons')} // Online variables should only be generated for units with restrictions
    / p_unit(unit, 'unitCount')
    * p_gnu(grid, node, unit, 'maxGen')
263
264
;
* -----------------------------------------------------------------------------
265
266
267
268
269
270
271
272
273
274
275
276
277
q_maxUpward(gnuft(grid, node, unit, f, t))${ [     ord(t) < tSolveFirst + settings('t_reserveLength')               // Unit is either providing
                                               and sum(restype, nuRescapable(restype, 'up', node, unit))         // upward reserves
                                             ] or
                                             [ uft_online(unit, f, t) and                                           // or the unit has an online varaible
                                                 (
                                                      [unit_minLoad(unit) and p_gnu(grid, node, unit, 'maxCons')]  // consuming units with min_load
                                                   or [p_gnu(grid, node, unit, 'maxGen')]                          // generators with an online variable
                                                 )
                                             ]
                                           }..
  + v_gen(grid, node, unit, f, t)                                                                                   // energy generation/consumption
  + sum( gngnu_constrainedOutputRatio(grid, node, grid_output, node_, unit),
         p_gnu(grid_output, node_, unit, 'cV') * v_gen(grid_output, node_, unit, f, t) )                            // considering output constraints (e.g. cV line)
Topi Rasku's avatar
Topi Rasku committed
278
  + sum(nuRescapable(restype, 'up', node, unit)$[unit_elec(unit) and ord(t) < tSolveFirst + settings('t_reserveLength')],                                                    // plus upward reserve participation
279
        v_reserve(restype, 'up', node, unit, f, t)                                                                  // (v_reserve can be used only if the unit can provide a particular reserve)
280
    )
281
282
  =L=                                                                                                               // must be less than available/online capacity
  - v_online(unit, f, t)${uft_online(unit, f, t) and p_gnu(grid, node, unit, 'maxCons')}       // Consuming units are restricted by their min. load (consuming is negative)
283
284
285
286
287
    / p_unit(unit, 'unitCount')
    * p_gnu(grid, node, unit, 'maxCons')
    * sum(effGroup, // Uses the minimum 'lb' for the current efficiency approximation
        + (p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)} + ts_effGroupUnit(effGroup, unit, 'lb', f, t))
      )
288
289
290
291
  + v_gen.up(grid, node, unit, f, t)${not uft_online(unit, f, t) and p_gnu(grid, node, unit, 'maxGen')}            // Generation units are restricted by their (online) capacity
  + v_online(unit, f, t)${uft_online(unit, f, t) and p_gnu(grid, node, unit, 'maxGen')}       // Consuming units are restricted by their min. load (consuming is negative)
    / p_unit(unit, 'unitCount')
    * p_gnu(grid, node, unit, 'maxGen')
292
293
;
* -----------------------------------------------------------------------------
294
q_startup(uft_online(unit, ft_dynamic(f, t))) ..
295
  + v_online(unit, f+cf(f,t), t)
Topi Rasku's avatar
Topi Rasku committed
296
297
298
299
  =E=
  + v_online(unit, f+pf(f,t), t+pt(t)) // This reaches to tFirstSolve when pt = -1
  + v_startup(unit, f+pf(f,t), t+pt(t))
  - v_shutdown(unit, f+pf(f,t), t+pt(t))
300
;
301
* -----------------------------------------------------------------------------
302
q_genRamp(gn(grid, node), m, unit, ft_dynamic(f, t))${gnuft_ramp(grid, node, unit, f, t)} ..
303
    + v_genRamp(grid, node, unit, f, t+pt(t))
304
305
306
307
    * ( p_gnu(grid, node, unit, 'maxGen') + p_gnu(grid, node, unit, 'maxCons') )  * 60 / 100 // Unit conversion from [p.u./min] to [MW/h]
    * p_stepLength(m, f+pf(f,t), t+pt(t))
    =E=
    // Change in generation over the time step
308
    + v_gen(grid, node, unit, f+cf(f,t), t)
309
310
311
    - v_gen(grid, node, unit, f+pf(f,t), t+pt(t))
    // Correction term to account for online variables and min loads
    - (
312
        + v_online(unit, f+cf(f,t), t)${uft_online(unit, f, t)}
313
314
315
316
        - v_online(unit, f+pf(f,t), t+pt(t))${uft_online(unit, f+pf(f,t), t+pt(t))}
      )
        / p_unit(unit, 'unitCount')
        * ( p_gnu(grid, node, unit, 'maxGen') - p_gnu(grid, node, unit, 'maxCons') )
317
        * sum(suft(effGroup, unit, f+cf(f,t), t), p_effGroupUnit(effGroup, unit, 'lb')) // Newly started units are assumed to start to their minload.
318
319
;
* -----------------------------------------------------------------------------
320
q_genRampChange(gn(grid, node), m, unit, ft_dynamic(f, t))${ gnuft_ramp(grid, node, unit, f, t) AND [ p_gnu(grid, node, unit, 'rampUpCost') OR p_gnu(grid, node, unit, 'rampDownCost') ]} ..
321
322
    + v_genRampChange(grid, node, unit, 'up', f, t+pt(t))
    - v_genRampChange(grid, node, unit, 'down', f, t+pt(t))
323
    =E=
324
325
    + v_genRamp(grid, node, unit, f, t+pt(t))
    - v_genRamp(grid, node, unit, f, t+pt(t))
326
327
;
* -----------------------------------------------------------------------------
328
q_conversionDirectInputOutput(suft(effDirect, unit, f, t)) ..
Juha Kiviluoma's avatar
Juha Kiviluoma committed
329
330
331
332
333
334
  - sum(gnu_input(grid, node, unit),
      + v_gen(grid, node, unit, f, t)
    )
  + sum(uFuel(unit, 'main', fuel),
      + v_fuelUse(fuel, unit, f, t)
    )
335
  =E=
Juha Kiviluoma's avatar
Juha Kiviluoma committed
336
337
  + sum(gnu_output(grid, node, unit),
      + v_gen(grid, node, unit, f, t)
338
          * (p_effUnit(effDirect, unit, effDirect, 'slope')${not ts_effUnit(effDirect, unit, effDirect, 'slope', f, t)} + ts_effUnit(effDirect, unit, effDirect, 'slope', f, t))
339
340
    )
  + v_online(unit, f, t)${uft_online(unit, f, t)}
341
342
343
      / p_unit(unit, 'unitCount')
      * sum( gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'maxGen') )
      * (p_effGroupUnit(effDirect, unit, 'section')${not ts_effUnit(effDirect, unit, effDirect, 'section', f, t)} + ts_effUnit(effDirect, unit, effDirect, 'section', f, t))
344
;
Juha Kiviluoma's avatar
Juha Kiviluoma committed
345
346
347
348
349
350
351
352
353
* -----------------------------------------------------------------------------
q_conversionSOS2InputIntermediate(suft(effGroup, unit, f, t))$effLambda(effGroup) ..
  - sum(gnu_input(grid, node, unit),
      + v_gen(grid, node, unit, f, t)
    )
  + sum(uFuel(unit, 'main', fuel),
      + v_fuelUse(fuel, unit, f, t)
    )
  =E=
354
355
356
357
358
359
360
  + ( + sum(effSelector$effGroupSelectorUnit(effGroup, unit, effSelector),
          + v_sos2(unit, f, t, effSelector)
              * ( p_effUnit(effGroup, unit, effSelector, 'op')${not ts_effUnit(effGroup, unit, effSelector, 'op', f, t)} + ts_effUnit(effGroup, unit, effSelector, 'op', f, t))
              * ( p_effUnit(effGroup, unit, effSelector, 'slope')${not ts_effUnit(effGroup, unit, effSelector, 'slope', f, t)} + ts_effUnit(effGroup, unit, effSelector, 'slope', f, t) )
        )
      + v_online(unit, f, t)
          * p_effGroupUnit(effGroup, unit, 'section')
Juha Kiviluoma's avatar
Juha Kiviluoma committed
361
    )
362
363
      / p_unit(unit, 'unitCount')
      * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'maxGen'))
364
;
365
* -----------------------------------------------------------------------------
Juha Kiviluoma's avatar
Juha Kiviluoma committed
366
q_conversionSOS2Constraint(suft(effGroup, unit, f, t))$effLambda(effGroup) ..
367
  + sum(effSelector$effGroupSelectorUnit(effGroup, unit, effSelector),
Juha Kiviluoma's avatar
Juha Kiviluoma committed
368
369
      + v_sos2(unit, f, t, effSelector)
    )
370
  =E=
371
  + v_online(unit, f, t)${uft_online(unit, f, t)}
372
*  + 1${not uft_online(unit, f, t)} // Should not be required, as effLambda implies online variables
373
;
374
* -----------------------------------------------------------------------------
Juha Kiviluoma's avatar
Juha Kiviluoma committed
375
q_conversionSOS2IntermediateOutput(suft(effGroup, unit, f, t))$effLambda(effGroup) ..
376
  + sum(effSelector$effGroupSelectorUnit(effGroup, unit, effSelector),
Juha Kiviluoma's avatar
Juha Kiviluoma committed
377
      + v_sos2(unit, f, t, effSelector)
378
      * (p_effUnit(effGroup, unit, effSelector, 'op')${not ts_effUnit(effGroup, unit, effSelector, 'op', f, t)} + ts_effUnit(effGroup, unit, effSelector, 'op', f, t))
Juha Kiviluoma's avatar
Juha Kiviluoma committed
379
    )
380
  / p_unit(unit, 'unitCount')
Juha Kiviluoma's avatar
Juha Kiviluoma committed
381
  * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'maxGen'))
382
  =E=
Juha Kiviluoma's avatar
Juha Kiviluoma committed
383
384
385
  + sum(gnu_output(grid, node, unit),
      + v_gen(grid, node, unit, f, t)
    )
386
;
387
* -----------------------------------------------------------------------------
388
q_outputRatioFixed(gngnu_fixedOutputRatio(grid, node, grid_, node_, unit), ft(f, t))${uft(unit, f, t)} ..
389
  + v_gen(grid, node, unit, f, t)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
390
391
      / p_gnu(grid, node, unit, 'cB')
  =E=
392
  + v_gen(grid_, node_, unit, f, t)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
393
394
395
      / p_gnu(grid_, node_, unit, 'cB')
;
* -----------------------------------------------------------------------------
396
397
q_outputRatioConstrained(gngnu_constrainedOutputRatio(grid, node, grid_, node_, unit), ft(f, t))${uft(unit, f, t)} ..
  + v_gen(grid, node, unit, f, t)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
398
      / p_gnu(grid, node, unit, 'cB')
399
  =G=
400
  + v_gen(grid_, node_, unit, f, t)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
401
      / p_gnu(grid_, node_, unit, 'cB')
402
;
403
* -----------------------------------------------------------------------------
404
q_stateSlack(gn_stateSlack(grid, node), slack, ft_dynamic(f, t))$p_gnBoundaryPropertiesForStates(grid, node, slack, 'slackCost') ..
405
  + v_stateSlack(grid, node, slack, f, t)
406
  =G=
407
408
409
410
411
  + p_slackDirection(slack) * (
      + v_state(grid, node, f, t)
      - p_gnBoundaryPropertiesForStates(grid, node, slack, 'constant')$p_gnBoundaryPropertiesForStates(grid, node, slack, 'useConstant')
      - ts_nodeState(grid, node, slack, f, t)$p_gnBoundaryPropertiesForStates(grid, node, slack, 'useTimeSeries')
    )
412
;
Juha Kiviluoma's avatar
Cleanup    
Juha Kiviluoma committed
413
* -----------------------------------------------------------------------------
414
415
q_stateUpwardLimit(gn_state(grid, node), m, ft_dynamic(f, t))$(    sum(gn2gnu(grid, node, grid_, node_output, unit)$(sum(restype, nuRescapable(restype, 'down', node_output, unit))), 1)  // nodes that have units with endogenous output with possible reserve provision
                                                        or sum(gn2gnu(grid_, node_input, grid, node, unit) $(sum(restype, nuRescapable(restype, 'down', node_input , unit))), 1)  // or nodes that have units with endogenous input with possible reserve provision
416
                                                      ) ..
417
418
419
420
421
422
423
424
425
426
427
428
  ( // Utilizable headroom in the state variable
      + p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'useConstant')   * p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'constant')
      + p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'useTimeSeries') * ts_nodeState(grid, node, 'upwardLimit', f, t)
      - v_state(grid, node, f, t)
  )
      * ( // Accounting for the energyStoredPerUnitOfState ...
          + p_gn(grid, node, 'energyStoredPerUnitOfState')
          + p_stepLength(m , f+pf(f,t), t+pt(t))
              * ( // ... and the change in energy losses from the node
                  + p_gn(grid, node, 'selfDischargeLoss')
                  + sum(to_node, p_gnn(grid, node, to_node, 'diffCoeff'))
                )
429
        )
430
431
432
  =G=
  + p_stepLength(m, f+pf(f,t), t+pt(t))
      * ( // Reserve provision from units that have output to this node
433
          + sum(gn2gnu(grid_, node_input, grid, node, unit)${uft(unit, f+pf(f,t), t+pt(t))},
Topi Rasku's avatar
Topi Rasku committed
434
              + sum(restype$[nuRescapable(restype, 'down', node_input, unit) and unit_elec(unit) and ord(t) < tSolveFirst + settings('t_reserveLength')], // Downward reserves from units that output energy to the node
435
                  + v_reserve(restype, 'down', node_input, unit, f+pf(f,t), t+pt(t))
436
                      / sum(effGroup${suft(effGroup, unit, f+pf(f,t), t+pt(t))}, (p_effGroupUnit(effGroup, unit, 'slope')${not ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t))} + ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t)))) // Efficiency approximated using maximum slope of effGroup?
437
438
439
                )
            )
          // Reserve provision from units that take input from this node
440
          + sum(gn2gnu(grid, node, grid_, node_output, unit)${uft(unit, f+pf(f,t), t+pt(t))},
Topi Rasku's avatar
Topi Rasku committed
441
              + sum(restype$[nuRescapable(restype, 'down', node_output, unit) and unit_elec(unit) and ord(t) < tSolveFirst + settings('t_reserveLength')], // Downward reserves from units that use the node as energy input
442
                  + v_reserve(restype, 'down', node_output, unit, f+pf(f,t), t+pt(t))
443
                      * sum(effGroup${suft(effGroup, unit, f+pf(f,t), t+pt(t))}, (p_effGroupUnit(effGroup, unit, 'slope')${not ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t))} + ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t)))) // Efficiency approximated using maximum slope of effGroup?
444
                )
445
            )
446
447
      // Here we could have a term for using the energy in the node to offer reserves as well as imports and exports of reserves, but as long as reserves are only
      // considered in power grids that do not have state variables, these terms are not needed. Earlier commit (29.11.2016) contains a draft of those terms.
448
        )
449
;
Juha Kiviluoma's avatar
Cleanup    
Juha Kiviluoma committed
450
* -----------------------------------------------------------------------------
451
452
q_stateDownwardLimit(gn_state(grid, node), m, ft_dynamic(f, t))$(    sum(gn2gnu(grid, node, grid_, node_output, unit)$(sum(restype, nuRescapable(restype, 'up', node_output, unit))), 1)  // nodes that have units with endogenous output with possible reserve provision
                                                          or sum(gn2gnu(grid_, node_input, grid, node, unit) $(sum(restype, nuRescapable(restype, 'up', node_input , unit))), 1)  // or nodes that have units with endogenous input with possible reserve provision
453
                                                        ) ..
454
455
456
457
458
459
460
461
462
463
464
465
  ( // Utilizable headroom in the state variable
      + v_state(grid, node, f, t)
      - p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'useConstant')   * p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'constant')
      - p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'useTimeSeries') * ts_nodeState(grid, node, 'downwardLimit', f, t)
  )
      * ( // Accounting for the energyStoredPerUnitOfState ...
          + p_gn(grid, node, 'energyStoredPerUnitOfState')
          + p_stepLength(m , f+pf(f,t), t+pt(t))
              * ( // ... and the change in energy losses from the node
                  + p_gn(grid, node, 'selfDischargeLoss')
                  + sum(to_node, p_gnn(grid, node, to_node, 'diffCoeff'))
                )
466
        )
467
468
469
  =G=
  + p_stepLength(m, f+pf(f,t), t+pt(t))
      * ( // Reserve provision from units that have output to this node
470
          + sum(gn2gnu(grid_, node_input, grid, node, unit)${uft(unit, f+pf(f,t), t+pt(t))},
Topi Rasku's avatar
Topi Rasku committed
471
              + sum(restype$[nuRescapable(restype, 'up', node_input, unit) and unit_elec(unit) and ord(t) < tSolveFirst + settings('t_reserveLength')], // Upward reserves from units that output energy to the node
472
                  + v_reserve(restype, 'up', node_input, unit, f+pf(f,t), t+pt(t))
473
                      / sum(effGroup${suft(effGroup, unit, f+pf(f,t), t+pt(t))}, (p_effGroupUnit(effGroup, unit, 'slope')${not ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t))} + ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t)))) // Efficiency approximated using maximum slope of effGroup?
474
475
476
                )
            )
          // Reserve provision from units that take input from this node
477
          + sum(gn2gnu(grid, node, grid_, node_output, unit)${uft(unit, f+pf(f,t), t+pt(t))},
Topi Rasku's avatar
Topi Rasku committed
478
              + sum(restype$[nuRescapable(restype, 'up', node_output, unit) and unit_elec(unit) and ord(t) < tSolveFirst + settings('t_reserveLength')], // Upward reserves from units that use the node as energy input
479
                  + v_reserve(restype, 'up', node_output, unit, f+pf(f,t), t+pt(t))
480
                      * sum(effGroup${suft(effGroup, unit, f+pf(f,t), t+pt(t))}, (p_effGroupUnit(effGroup, unit, 'slope')${not ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t))} + ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t)))) // Efficiency approximated using maximum slope of effGroup?
481
                )
482
483
484
            )
      // Here we could have a term for using the energy in the node to offer reserves as well as imports and exports of reserves, but as long as reserves are only
      // considered in power grids that do not have state variables, these terms are not needed. Earlier commit (29.11.2016) contains a draft of those terms.
485
        )
486
487
;
* -----------------------------------------------------------------------------
488
q_boundState(gnn_boundState(grid, node, node_), m, ft_dynamic(f, t)) ..
489
    + v_state(grid, node, f, t)   // The state of the first node sets the upper limit of the second
490
    + ( // Downward reserve provided by units in node
491
492
*        - sum(nuRescapable(restype, 'down', node, unit)${nuft(unit, f+pf(f,t), t+pt(t))},
*            + v_reserve(restype, 'down', node, unit, f+pf(f,t), t+pt(t))
493
494
*          )
        // Upwards reserve provided by input units
Topi Rasku's avatar
Topi Rasku committed
495
        - sum(nuRescapable(restype, 'up', node_input, unit)${sum(grid_, gn2gnu(grid_, node_input, grid, node, unit)) AND uft(unit, f+pf(f,t), t+pt(t)) and ord(t) < tSolveFirst + settings('t_reserveLength')},
496
            + v_reserve(restype, 'up', node_input, unit, f+pf(f,t), t+pt(t))
497
                / sum(effGroup${suft(effGroup, unit, f+pf(f,t), t+pt(t))}, (p_effGroupUnit(effGroup, unit, 'slope')${not ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t))} + ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t)))) // Efficiency approximated using maximum slope of effGroup?
498
          )
499
        // Upwards reserve providewd by output units
Topi Rasku's avatar
Topi Rasku committed
500
        - sum(nuRescapable(restype, 'up', node_output, unit)${sum(grid_, gn2gnu(grid, node, grid_, node_output, unit)) AND uft(unit, f+pf(f,t), t+pt(t)) and ord(t) < tSolveFirst + settings('t_reserveLength')},
501
            + v_reserve(restype, 'up', node_output, unit, f+pf(f,t), t+pt(t))
502
                / sum(effGroup${suft(effGroup, unit, f+pf(f,t), t+pt(t))}, (p_effGroupUnit(effGroup, unit, 'slope')${not ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t))} + ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t)))) // Efficiency approximated using maximum slope of effGroup?
503
          )
504
505
        // Here we could have a term for using the energy in the node to offer reserves as well as imports and exports of reserves, but as long as reserves are only
        // considered in power grids that do not have state variables, these terms are not needed. Earlier commit (16.2.2017) contains a draft of those terms.
506
      )
507
        / (p_gn(grid, node, 'energyStoredPerUnitOfState') )                                     // Divide by the stored energy in the node per unit of v_state to obtain same unit as the v_state
508
509
510
511
        * p_stepLength(m, f+pf(f,t), t+pt(t))                                                   // Multiply with time step to obtain change in state over the step
    =G=
    + v_state(grid, node_, f, t)
    + p_gnn(grid, node, node_, 'boundStateOffset')                                              // Affected by the offset parameter
512
    + (  // Possible reserve by this node
513
514
*        + sum(nuRescapable(restype, 'up', node_, unit)${nuft(node, unit, f+pf(f,t), t+pt(t))},
*            + v_reserve(restype, 'up', node_, unit, f+pf(f,t), t+pt(t))
515
516
*          )
        // Possible reserve by input node
Topi Rasku's avatar
Topi Rasku committed
517
        + sum(nuRescapable(restype, 'down', node_input, unit)${sum(grid_, gn2gnu(grid_, node_input, grid, node_, unit)) AND uft(unit, f+pf(f,t), t+pt(t)) and ord(t) < tSolveFirst + settings('t_reserveLength')},
518
            + v_reserve(restype, 'down', node_input, unit, f+pf(f,t), t+pt(t))               // NOTE! If elec-elec conversion, this might result in weird reserve requirements!
519
                / sum(effGroup${suft(effGroup, unit, f+pf(f,t), t+pt(t))}, (p_effGroupUnit(effGroup, unit, 'slope')${not ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t))} + ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t)))) // Efficiency approximated using maximum slope of effGroup?
520
          )
521
        // Possible reserve by output node
Topi Rasku's avatar
Topi Rasku committed
522
        + sum(nuRescapable(restype, 'down', node_output, unit)${sum(grid_, gn2gnu(grid, node_, grid_, node_output, unit)) AND uft(unit, f+pf(f,t), t+pt(t)) and ord(t) < tSolveFirst + settings('t_reserveLength')},
523
            + v_reserve(restype, 'down', node_output, unit, f+pf(f,t), t+pt(t))               // NOTE! If elec-elec conversion, this might result in weird reserve requirements!
524
                / sum(effGroup${suft(effGroup, unit, f+pf(f,t), t+pt(t))}, (p_effGroupUnit(effGroup, unit, 'slope')${not ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t))} + ts_effGroupUnit(effGroup, unit, 'slope', f+pf(f,t), t+pt(t)))) // Efficiency approximated using maximum slope of effGroup?
525
          )
526
527
        // Here we could have a term for using the energy in the node to offer reserves as well as imports and exports of reserves, but as long as reserves are only
        // considered in power grids that do not have state variables, these terms are not needed. Earlier commit (16.2.2017) contains a draft of those terms.
528
      )
529
        / (p_gn(grid, node_, 'energyStoredPerUnitOfState') )                                    // Divide by the stored energy in the node per unit of v_state to obtain same unit as the v_state
530
        * p_stepLength(m, f+pf(f,t), t+pt(t))                                                   // Multiply with time step to obtain change in state over the step
531
;
532
* -----------------------------------------------------------------------------
533
534
535
536
537
538
539
q_boundStateMaxDiff(gnn_state(grid, node, node_), m, ft_dynamic(f, t))$p_gnn(grid, node, node_, 'boundStateMaxDiff') ..
    + v_state(grid, node, f, t)   // The state of the first node sets the lower limit of the second
    =L=
    + v_state(grid, node_, f, t)
    + p_gnn(grid, node, node_, 'boundStateMaxDiff')                                              // Affected by the maximum difference parameter
;
* -----------------------------------------------------------------------------
540
q_boundCyclic(gn_state(grid, node), mft(m, f, t), t_)${  p_gn(grid, node, 'boundCyclic')         // Bind variables if parameter found
541
                                                        AND tSolveFirst = mSettings(m, 't_start') // For the very first model solve only
542
                                                        AND mftStart(m, f+pf(f,t), t)                   // Use only the starting time step of the model solve
543
544
                                                        AND mftLastSteps(m, f, t_)              // Use only the ending time step of the model solve
                                                        }..
545
    + v_state(grid, node, f+pf(f,t), t)
546
547
548
    =E=
    + v_state(grid, node, f, t_)
;
549
550
551
552
* -----------------------------------------------------------------------------
q_bidirectionalTransfer(gn2n_bidirectional(grid, node, node_), ft(f, t))${p_gnn(grid, node, node_, 'transferCapBidirectional')} ..
    + v_transfer(grid, node, node_, f, t) // Transfers in one direction
    + v_transfer(grid, node_, node, f, t) // Transfers in the other direction
553
554
555
    + sum(restypeDirection(restype, up_down)${restypeDirectionNode(restype, up_down, node) AND restypeDirectionNode(restype, up_down, node_)},
        + v_resTransfer(restype, up_down, node, node_, f, t) // Reserve transfers in one direction
        + v_resTransfer(restype, up_down, node_, node, f, t) // Reserve transfers in the other direction
556
557
558
559
      )
    =L=
    p_gnn(grid, node, node_, 'transferCapBidirectional')
;
560