3c_inputsLoop.gms 27.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
$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
17

18
* =============================================================================
19
* --- Update the Forecast Data ------------------------------------------------
20
* =============================================================================
21

22
put log 'ord tSolve: ';
23
put log ord(tSolve):0:0 /;
24
putclose log;
25

26
27
// Determine the necessary horizon for updating data
option clear = tmp;
28
29
tmp = max(  mSettings(mSolve, 't_forecastLengthUnchanging'),
            mSettings(mSolve, 't_forecastLengthDecreasesFrom')
30
31
32
33
34
35
36
            );

// Find time steps until the forecast horizon
option clear = tt_forecast;
tt_forecast(t_current(t))
    ${ ord(t) <= tSolveFirst + tmp }
    = yes;
37

38
if (ord(tSolve) = tForecastNext(mSolve) - mSettings(mSolve, 't_forecastJump'), // tForecastNext updated already in periodicLoop!
39

40
41
42
    // Update ts_unit
    if (mTimeseries_loop_read(mSolve, 'ts_unit'),
        put_utility 'gdxin' / '%input_dir%/ts_unit/' tSolve.tl:0 '.gdx';
43
        execute_load ts_unit_update=ts_unit;
44
        ts_unit(unit_timeseries(unit), param_unit, f_solve(f), tt_forecast(t)) // Only update if time series enabled for the unit
45
            ${not mf_realization(mSolve, f) // Realization not updated
46
              and (mSettings(mSolve, 'onlyExistingForecasts')
47
                   -> ts_unit_update(unit, param_unit, f, t)) // Update only existing values (zeroes need to be EPS)
48
                }
49
            = ts_unit_update(unit, param_unit, f, t);
50
    ); // END if('ts_unit')
51
52
53
54
$ontext
* !!! NOTE !!!
* These probably shouldn't be read at all, as p_effUnit and p_effGroupUnit are
* not input data, but calculated based on p_unit
55
    // Update ts_effUnit
56
    if (mTimeseries_loop_read(mSolve, 'ts_effUnit'),
57
        put_utility 'gdxin' / '%input_dir%/ts_effUnit/' tSolve.tl:0 '.gdx';
58
        execute_load ts_effUnit_update=ts_effUnit;
59
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f_solve(f), tt_forecast(t)) // Only update if time series enabled for the unit
60
            ${  not mf_realization(mSolve, f) // Realization not updated
61
*               and ts_effUnit_update(effSelector, unit, effSelector, param_eff, f, t) // Update only existing values (zeroes need to be EPS)
62
                }
63
            = ts_effUnit_update(effSelector, unit, effSelector, param_eff, f, t);
64
    ); // END if('ts_effUnit')
65
66

    // Update ts_effGroupUnit
67
    if (mTimeseries_loop_read(mSolve, 'ts_effGroupUnit'),
68
        put_utility 'gdxin' / '%input_dir%/ts_effGroupUnit/' tSolve.tl:0 '.gdx';
69
        execute_load ts_effGroupUnit_update=ts_effGroupUnit;
70
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f_solve(f), tt_forecast(t)) // Only update if time series enabled for the unit
71
            ${  not mf_realization(mSolve, f) // Realization not updated
72
*               and ts_effGroupUnit_update(effSelector, unit, param_eff, f, t) // Update only existing values (zeroes need to be EPS)
73
                }
74
            = ts_effGroupUnit_update(effSelector, unit, param_eff, f, t);
75
76
    ); // END if('ts_effGroupUnit')
$offtext
77
    // Update ts_influx
78
    if (mTimeseries_loop_read(mSolve, 'ts_influx'),
79
        put_utility 'gdxin' / '%input_dir%/ts_influx/' tSolve.tl:0 '.gdx';
80
81
        execute_load ts_influx_update=ts_influx;
        ts_influx(gn(grid, node), f_solve(f), tt_forecast(t))
82
            ${  not mf_realization(mSolve, f) // Realization not updated
83
                and (mSettings(mSolve, 'onlyExistingForecasts')
84
                     -> ts_influx_update(grid, node, f, t)) // Update only existing values (zeroes need to be EPS)
85
                }
86
87
            = ts_influx_update(grid, node, f, t);
    ); // END if('ts_influx')
88
89

    // Update ts_cf
90
    if (mTimeseries_loop_read(mSolve, 'ts_cf'),
91
        put_utility 'gdxin' / '%input_dir%/ts_cf/' tSolve.tl:0 '.gdx';
92
93
        execute_load ts_cf_update=ts_cf;
        ts_cf(flowNode(flow, node), f_solve(f), tt_forecast(t))
94
            ${  not mf_realization(mSolve, f) // Realization not updated
95
                and (mSettings(mSolve, 'onlyExistingForecasts')
96
                     -> ts_cf_update(flow, node, f, t)) // Update only existing values (zeroes need to be EPS)
97
                }
98
99
            = ts_cf_update(flow, node, f, t);
    ); // END if('ts_cf')
100
101

    // Update ts_reserveDemand
102
    if (mTimeseries_loop_read(mSolve, 'ts_reserveDemand'),
103
        put_utility 'gdxin' / '%input_dir%/ts_reserveDemand/' tSolve.tl:0 '.gdx';
104
        execute_load ts_reserveDemand_update=ts_reserveDemand;
105
        ts_reserveDemand(restypeDirectionGroup(restype, up_down, group), f_solve(f), tt_forecast(t))
106
            ${  not mf_realization(mSolve, f) // Realization not updated
107
                and (mSettings(mSolve, 'onlyExistingForecasts')
108
                     -> ts_reserveDemand_update(restype, up_down, group, f, t)) // Update only existing values (zeroes need to be EPS)
109
                }
110
            = ts_reserveDemand_update(restype, up_down, group, f, t);
111
    ); // END if('ts_reserveDemand')
112
113

    // Update ts_node
114
    if (mTimeseries_loop_read(mSolve, 'ts_node'),
115
        put_utility 'gdxin' / '%input_dir%/ts_node/' tSolve.tl:0 '.gdx';
116
117
        execute_load ts_node_update=ts_node;
        ts_node(gn(grid, node), param_gnBoundaryTypes, f_solve(f), tt_forecast(t))
118
            ${  not mf_realization(mSolve, f) // Realization not updated
119
                and (mSettings(mSolve, 'onlyExistingForecasts')
120
                     -> ts_node_update(grid, node, param_gnBoundaryTypes, f ,t)) // Update only existing values (zeroes need to be EPS)
121
                }
122
123
124
125
126
            = ts_node_update(grid, node, param_gnBoundaryTypes, f, t);
    ); // END if('ts_node')

* --- NO FORECAST DIMENSION, SHOULD THESE BE HANDLED SEPARATELY? --------------
// Currently, only updated until the forecast horizon, but is this correct?
127

128
129
    // Update ts_priceChange
    if (mTimeseries_loop_read(mSolve, 'ts_priceChange'),
130
131
        put log '!!! Abort: mTimeseries_loop_read(mSolve, ts_priceChange) currently not working!' /;
        abort "mTimeseries_loop_read(mSolve, ts_priceChange) currently not working!";
132
133
134
135
136
137
        put_utility 'gdxin' / '%input_dir%/ts_priceChange/' tSolve.tl:0 '.gdx';
        execute_load ts_priceChange_update=ts_priceChange;
        ts_priceChange(node, tt_forecast(t))
*            ${ ts_priceChange_update(fuel, t) } // Update only existing values (zeroes need to be EPS)
            = ts_priceChange_update(node, t);
    ); // END if('ts_priceChange')
138
139

    // Update ts_unavailability
140
    if (mTimeseries_loop_read(mSolve, 'ts_unavailability'),
141
        put_utility 'gdxin' / '%input_dir%/ts_unavailability/' tSolve.tl:0 '.gdx';
142
143
        execute_load ts_unavailability_update=ts_unavailability;
        ts_unavailability(unit, tt_forecast(t))
144
*            ${ ts_unavailability_update(unit, t) } // Update only existing values (zeroes need to be EPS)
145
146
            = ts_unavailability_update(unit, t);
    ); // END if('ts_unavailability')
147

148
); // END if(tForecastNext)
149

150
* =============================================================================
151
* --- Optional forecast improvement -------------------------------------------
152
153
* =============================================================================

154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
if(mSettings(mSolve, 't_improveForecast'),
* Linear improvement of the central forecast towards the realized forecast,
* while preserving the difference between the central forecast and the
* remaining forecasts

    // Determine the set of improved time steps
    option clear = tt;
    tt(tt_forecast(t))
        ${ ord(t) <= tSolveFirst + mSettings(mSolve, 't_improveForecast') }
        = yes;

    // Temporary forecast displacement to reach the central forecast
    option clear = ddf;
    ddf(f_solve(f))
        ${ not mf_central(mSolve, f) }
        = sum(mf_central(mSolve, f_), ord(f_) - ord(f));

    // Temporary forecast displacement to reach the realized forecast
    option clear = ddf_;
    ddf_(f_solve(f))
        ${ not mf_realization(mSolve, f) }
        = sum(mf_realization(mSolve, f_), ord(f_) - ord(f));

* --- Calculate the other forecasts relative to the central one ---------------

    loop(f_solve(f)${ not mf_realization(mSolve, f) and not mf_central(mSolve, f) },
        // ts_unit
181
182
183
184
        ts_unit(unit_timeseries(unit), param_unit, f, tt(t))// Only update for units with time series enabled
            = ts_unit(unit, param_unit, f, t) - ts_unit(unit, param_unit, f+ddf(f), t);
$ontext
* Should these be handled here at all? See above note
185
        // ts_effUnit
186
187
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f, tt(t)) // Only update for units with time series enabled
            = ts_effUnit(effSelector, unit, effSelector, param_eff, f, t) - ts_effUnit(effSelector, unit, effSelector, param_eff, f+ddf(f), t);
188
        // ts_effGroupUnit
189
190
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f, tt(t)) // Only update for units with time series enabled
            = ts_effGroupUnit(effSelector, unit, param_eff, f, t) - ts_effGroupUnit(effSelector, unit, param_eff, f+ddf(f), t);
191
192
193
194
195
196
197
198
$offtext
        // ts_influx
        ts_influx(gn(grid, node), f, tt(t))
            = ts_influx(grid, node, f, t) - ts_influx(grid, node, f+ddf(f), t);
        // ts_cf
        ts_cf(flowNode(flow, node), f, tt(t))
            = ts_cf(flow, node, f, t) - ts_cf(flow, node, f+ddf(f), t);
        // ts_reserveDemand
199
200
        ts_reserveDemand(restypeDirectionGroup(restype, up_down, group), f, tt(t))
            = ts_reserveDemand(restype, up_down, group, f, t) - ts_reserveDemand(restype, up_down, group, f+ddf(f), t);
201
202
203
204
205
206
207
208
209
        // ts_node
        ts_node(gn(grid, node), param_gnBoundaryTypes, f, tt(t))
            = ts_node(grid, node, param_gnBoundaryTypes, f, t) - ts_node(grid, node, param_gnBoundaryTypes, f+ddf(f), t);
    ); // END loop(f_solve)

* --- Linear improvement of the central forecast ------------------------------

    loop(mf_central(mSolve, f),
        // ts_unit
210
        ts_unit(unit_timeseries(unit), param_unit, f, tt(t)) // Only update for units with time series enabled
211
            = [ + (ord(t) - tSolveFirst)
212
                    * ts_unit(unit, param_unit, f, t)
213
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
214
                    * ts_unit(unit, param_unit, f+ddf_(f), t)
215
                ] / mSettings(mSolve, 't_improveForecast');
216
217
$ontext
* Should these be handled here at all? See above note
218
        // ts_effUnit
219
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f, tt(t)) // Only update for units with time series enabled
220
            = [ + (ord(t) - tSolveFirst)
221
                    * ts_effUnit(effSelector, unit, effSelector, param_eff, f, t)
222
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
223
                    * ts_effUnit(effSelector, unit, effSelector, param_eff, f+ddf_(f), t)
224
225
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_effGroupUnit
226
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f, tt(t)) // Only update for units with time series enabled
227
            = [ + (ord(t) - tSolveFirst)
228
                    * ts_effGroupUnit(effSelector, unit, param_eff, f, t)
229
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
230
                    * ts_effGroupUnit(effSelector, unit, param_eff, f+ddf_(f), t)
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
                ] / mSettings(mSolve, 't_improveForecast');
$offtext
        // ts_influx
        ts_influx(gn(grid, node), f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_influx(grid, node, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_influx(grid, node, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_cf
        ts_cf(flowNode(flow, node), f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_cf(flow, node, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_cf(flow, node, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_reserveDemand
248
        ts_reserveDemand(restypeDirectionGroup(restype, up_down, group), f, tt(t))
249
            = [ + (ord(t) - tSolveFirst)
250
                    * ts_reserveDemand(restype, up_down, group, f, t)
251
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
252
                    * ts_reserveDemand(restype, up_down, group, f+ddf_(f), t)
253
254
255
256
257
258
259
260
261
262
263
264
265
266
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_node
        ts_node(gn(grid, node), param_gnBoundaryTypes, f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_node(grid, node, param_gnBoundaryTypes, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_node(grid, node, param_gnBoundaryTypes, f+ddf_(f), t)
                ] / mSettings(mSolve, 't_improveForecast');
    ); // END loop(mf_central)

* --- Recalculate the other forecasts based on the improved central forecast --

    loop(f_solve(f)${ not mf_realization(mSolve, f) and not mf_central(mSolve, f) },
        // ts_unit
267
268
269
270
        ts_unit(unit_timeseries(unit), param_unit, f, tt(t)) // Only update for units with time series enabled
            = ts_unit(unit, param_unit, f, t) + ts_unit(unit, param_unit, f+ddf(f), t);
$ontext
* Should these be handled here at all? See above note
271
        // ts_effUnit
272
273
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f, tt(t)) // Only update for units with time series enabled
            = ts_effUnit(effSelector, unit, effSelector, param_eff, f, t) + ts_effUnit(effSelector, unit, effSelector, param_eff, f+ddf(f), t);
274
        // ts_effGroupUnit
275
276
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f, tt(t)) // Only update for units with time series enabled
            = ts_effGroupUnit(effSelector, unit, param_eff, f, t) + ts_effGroupUnit(effSelector, unit, param_eff, f+ddf(f), t);
277
278
279
280
281
282
283
284
$offtext
        // ts_influx
        ts_influx(gn(grid, node), f, tt(t))
            = ts_influx(grid, node, f, t) + ts_influx(grid, node, f+ddf(f), t);
        // ts_cf
        ts_cf(flowNode(flow, node), f, tt(t))
            = max(min(ts_cf(flow, node, f, t) + ts_cf(flow, node, f+ddf(f), t), 1), 0); // Ensure that capacity factor forecasts remain between 0-1
        // ts_reserveDemand
285
286
        ts_reserveDemand(restypeDirectionGroup(restype, up_down, group), f, tt(t))
            = max(ts_reserveDemand(restype, up_down, group, f, t) + ts_reserveDemand(restype, up_down, group, f+ddf(f), t), 0); // Ensure that reserve demand forecasts remains positive
287
288
289
290
291
292
        // ts_node
        ts_node(gn(grid, node), param_gnBoundaryTypes, f, tt(t))
            = ts_node(grid, node, param_gnBoundaryTypes, f, t) + ts_node(grid, node, param_gnBoundaryTypes, f+ddf(f), t);
    ); // END loop(f_solve)

); // END if(t_improveForecast)
293

294
295
296
297
298
299
300
301
302
303
304
* =============================================================================
* --- Aggregate time series data for the time intervals -----------------------
* =============================================================================

// Loop over the defined blocks of intervals
loop(cc(counter),

    // Retrieve interval block time steps
    option clear = tt_interval;
    tt_interval(t) = tt_block(counter, t);

305
306
    // Select and average time series data matching the intervals
    ts_unit_(unit_timeseries(unit), param_unit, ft(f, tt_interval(t)))
307
308
        = sum(tt_aggregate(t, t_),
            ts_unit(unit, param_unit, f, t_+dt_circular(t_))
309
310
            )
            / mInterval(mSolve, 'stepsPerInterval', counter);
311
312
$ontext
* Should these be handled here at all? See above comment
313
    ts_effUnit_(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, ft(f, tt_interval(t)))
314
315
        = sum(tt_aggregate(t, t_),
            ts_effUnit(effSelector, unit, effSelector, param_eff, f, t_+dt_circular(t_))
316
317
318
            )
            / mInterval(mSolve, 'stepsPerInterval', counter);
    ts_effGroupUnit_(effSelector, unit_timeseries(unit), param_eff, ft(f, tt_interval(t)))
319
320
        = sum(tt_aggregate(t, t_),
            ts_effGroupUnit(effSelector, unit, param_eff, f, t_+dt_circular(t_))
321
322
            )
            / mInterval(mSolve, 'stepsPerInterval', counter);
323
$offtext
324
    ts_influx_(gn(grid, node), sft(s, f, tt_interval(t)))
325
        = sum(tt_aggregate(t, t_),
326
            ts_influx(grid, node,
327
328
329
330
                f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_influx'))
                     + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_influx')),
                t_+ (+ dt_scenarioOffset(grid, node, 'ts_influx', s)
                     + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_influx'))))
331
332
            )
            / mInterval(mSolve, 'stepsPerInterval', counter);
333
    ts_cf_(flowNode(flow, node), sft(s, f, tt_interval(t)))
334
        = sum(tt_aggregate(t, t_),
335
            ts_cf(flow, node,
336
337
338
339
                f + (  df_realization(f, t)$(not gn_forecasts(flow, node, 'ts_cf'))
                     + df_scenario(f, t)$gn_scenarios(flow, node, 'ts_cf')),
                t_+ (  dt_scenarioOffset(flow, node, 'ts_cf', s)
                     + dt_circular(t_)$(not gn_scenarios(flow, node, 'ts_cf'))))
340
341
342
            )
            / mInterval(mSolve, 'stepsPerInterval', counter);
    // Reserves relevant only until reserve_length
343
344
    ts_reserveDemand_(restypeDirectionGroup(restype, up_down, group), ft(f, tt_interval(t)))
      ${ord(t) <= tSolveFirst + p_groupReserves(group, restype, 'reserve_length')  }
345
        = sum(tt_aggregate(t, t_),
346
            ts_reserveDemand(restype, up_down, group,
347
348
                f + (  df_realization(f, t)${not sum(gnGroup(grid, node, group), gn_forecasts(restype, node, 'ts_reserveDemand'))}
                     + df_scenario(f, t)${sum(gnGroup(grid, node, group), gn_scenarios(restype, node, 'ts_reserveDemand'))} ),                t_+ dt_circular(t_))
349
350
            )
            / mInterval(mSolve, 'stepsPerInterval', counter);
351
    ts_node_(gn_state(grid, node), param_gnBoundaryTypes, sft(s, f, tt_interval(t)))
352
353
      $p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')
           // Take average if not a limit type
354
        = (sum(tt_aggregate(t, t_),
355
                ts_node(grid, node, param_gnBoundaryTypes,
356
357
                    f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_node'))
                         + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
358
359
                    t_+ (   + dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                            + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node'))))
360
361
362
363
364
            )
            / mInterval(mSolve, 'stepsPerInterval', counter))$( not (sameas(param_gnBoundaryTypes, 'upwardLimit')
                                                                or sameas(param_gnBoundaryTypes, 'downwardLimit')
                                                                or slack(param_gnBoundaryTypes)))
          // Maximum lower limit
365
          + smax(tt_aggregate(t, t_),
366
                ts_node(grid, node, param_gnBoundaryTypes,
367
368
                    f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_node'))
                         + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
369
370
                    t_+ (   + dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                            + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node'))))
371
372
373
                )
                $(sameas(param_gnBoundaryTypes, 'downwardLimit') or downwardSlack(param_gnBoundaryTypes))
          // Minimum upper limit
374
          + smin(tt_aggregate(t, t_),
375
                ts_node(grid, node, param_gnBoundaryTypes,
376
377
                    f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_node'))
                         + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_node')),
378
379
                    t_+ (   + dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                            + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node'))))
380
381
382
                )
                $(sameas(param_gnBoundaryTypes, 'upwardLimit') or upwardSlack(param_gnBoundaryTypes));

383
    // Node price time series
384
385
    ts_vomCost_(gnu(grid, node, unit), tt_interval(t))
        = + p_gnu(grid, node, unit, 'vomCosts')
386
387
388
389
390
391
392
393
394
395
396
397
          // input node cost
          + p_price(node, 'price')${un_commodity_in(unit, node) and p_price(node, 'useConstant')}
          + sum(tt_aggregate(t, t_)${un_commodity_in(unit, node) and p_price(node, 'useTimeSeries')},
              + ts_price(node, t_+dt_circular(t_))
              )
              / mInterval(mSolve, 'stepsPerInterval', counter)
          // output node cost (if price > 0 --> ts_vomCost_ < 0, i.e. considered as revenue)
          - p_price(node, 'price')${un_commodity_out(unit, node) and p_price(node, 'useConstant')}
          - sum(tt_aggregate(t, t_)${un_commodity_out(unit, node) and p_price(node, 'useTimeSeries')},
              + ts_price(node, t_+dt_circular(t_))
              )
              / mInterval(mSolve, 'stepsPerInterval', counter)
398
          // emission cost
399
400
401
402
          + sum(emission$p_unitEmissionCost(unit, node, emission), // Emission taxes
              + p_unitEmissionCost(unit, node, emission)
            ); // END sum(emission)

403
    p_unStartup(unit, node, starttype)$p_uStartupfuel(unit, node, 'fixedFuelFraction')
404
      =
405
406
        + p_uStartup(unit, starttype, 'consumption')
            * p_uStartupfuel(unit, node, 'fixedFuelFraction');
407
408
409
410

    // Calculating startup cost time series
    ts_startupCost_(unit, starttype, tt_interval(t))
      =
411
        + p_uStartup(unit, starttype, 'cost') // CUR/start-up
412
        // Start-up fuel and emission costs
413
414
        + sum(nu(node,unit)$p_unStartup(unit, node, starttype),
            + p_unStartup(unit, node, starttype) // MWh/start-up
415
              * [
416
417
418
                  + p_price(node, 'price')$p_price(node, 'useConstant') // CUR/MWh
                  + sum(tt_aggregate(t, t_)$p_price(node, 'useTimeseries'),
                      + ts_price(node, t_+dt_circular(t_)) // CUR/MWh
419
420
421
                    )
                    / mInterval(mSolve, 'stepsPerInterval', counter)
                ] // END * p_uStartup
422
423
424
425
          ) // END sum(node)
        + sum((nu(node, unit), emission)$p_unitEmissionCost(unit, node, emission),
            + p_unStartup(unit, node, starttype) // MWh/start-up
              * p_unitEmissionCost(unit, node, emission) // CUR/MWh
Topi Rasku's avatar
Topi Rasku committed
426
          ); // END sum(nu, emission)
427
428

    // `storageValue`
429
    ts_storageValue_(gn_state(grid, node), sft(s, f, tt_interval(t)))${ p_gn(grid, node, 'storageValueUseTimeSeries') }
430
431
432
433
434
435
436
437
438
        = sum(tt_aggregate(t, t_),
            ts_storageValue(grid, node,
                f + (  df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_storageValue'))
                     + df_scenario(f, t)$gn_scenarios(grid, node, 'ts_storageValue')),
                t_+ (+ dt_scenarioOffset(grid, node, 'ts_storageValue', s)
                     + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_storageValue'))))
            )
            / mInterval(mSolve, 'stepsPerInterval', counter);

439
440
); // END loop(counter)

441

442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
* --- Process unit time series data -------------------------------------------

// Calculate time series form parameters for units using direct input output conversion without online variable
// Always constant 'lb', 'rb', and 'section', so need only to define 'slope'.
loop(effGroupSelectorUnit(effDirectOff, unit_timeseries(unit), effDirectOff_),
    ts_effUnit(effDirectOff, unit, effDirectOff_, 'slope', ft(f, t))
        ${ sum(eff, ts_unit(unit, eff, f, t)) } // NOTE!!! Averages the slope over all available data.
        = sum(eff${ts_unit(unit, eff, f, t)}, 1 / ts_unit(unit, eff, f, t))
            / sum(eff${ts_unit(unit, eff, f, t)}, 1);
); // END loop(effGroupSelectorUnit)

// NOTE! Using the same methodology for the directOn and lambda approximations in time series form might require looping over ft(f,t) to find the min and max 'eff' and 'rb'
// Alternatively, one might require that the 'rb' is defined in a similar structure, so that the max 'rb' is located in the same index for all ft(f,t)

// Calculate unit wide parameters for each efficiency group
loop(effLevelGroupUnit(effLevel, effGroup, unit)${  mSettingsEff(mSolve, effLevel)
                                                    and p_unit(unit, 'useTimeseries')
                                                    },
    ts_effGroupUnit(effGroup, unit, 'lb', ft(f, t))${   sum(effSelector, ts_effUnit(effGroup, unit, effSelector, 'lb', f, t))}
        = smin(effSelector${effGroupSelectorUnit(effGroup, unit, effSelector)}, ts_effUnit(effGroup, unit, effSelector, 'lb', f, t));
    ts_effGroupUnit(effGroup, unit, 'slope', ft(f, t))${sum(effSelector, ts_effUnit(effGroup, unit, effSelector, 'slope', f, t))}
        = smin(effSelector$effGroupSelectorUnit(effGroup, unit, effSelector), ts_effUnit(effGroup, unit, effSelector, 'slope', f, t)); // Uses maximum efficiency for the group
); // END loop(effLevelGroupUnit)

466

Erkka Rinne's avatar
Erkka Rinne committed
467
468
469
470
* =============================================================================
* --- Input data processing ---------------------------------------------------
* =============================================================================

Erkka Rinne's avatar
Erkka Rinne committed
471
$ifthen.scenarios defined scenario
Erkka Rinne's avatar
Erkka Rinne committed
472
* --- Scenario reduction ------------------------------------------------------
473
if(active(mSolve, 'scenred') and mSettings('schedule', 'scenarios') > 1,
Erkka Rinne's avatar
Erkka Rinne committed
474
475
    $$include 'inc/scenred.gms'
);
Erkka Rinne's avatar
Erkka Rinne committed
476
$endif.scenarios
Erkka Rinne's avatar
Erkka Rinne committed
477

478
479
480
481
* --- Update probabilities ----------------------------------------------------
Option clear = p_msft_probability;
p_msft_probability(msft(mSolve, s, f, t))
    = p_mfProbability(mSolve, f)
482
483
484
485
        / sum(f_$ft(f_, t),
              p_mfProbability(mSolve, f_)
          ) * p_msProbability(mSolve, s)
            * p_msWeight(mSolve, s);
486
487


488
489
* --- Calculate sample displacements ------------------------------------------
Options clear = ds, clear = ds_state;
Erkka Rinne's avatar
Erkka Rinne committed
490
491
loop((mst_start(mSolve, s, t), ss(s, s_)),
    ds(s, t) = -(ord(s) - ord(s_));
492
493
494
495
496
497
    ds_state(gn_state(grid, node), s, t)
      ${not sum(s__, gnss_bound(grid, node, s__, s))
        and not sum(s__, gnss_bound(grid, node, s, s__))} = ds(s, t);
);


Erkka Rinne's avatar
Erkka Rinne committed
498
499
* --- Smooting of stochastic scenarios ----------------------------------------
$ontext
500
501
Smoothen the scenarios following the methodology presented in [1, p. 443].
This avoids a discontinuity `jump' after the initial sample.
Erkka Rinne's avatar
Erkka Rinne committed
502
503
504
505
506
507

[1] A. Helseth, B. Mo, A. Lote Henden, and G. Warland, "Detailed long-term hydro-
    thermal scheduling for expansion planning in the Nordic power system," IET Gener.
    Transm. Distrib., vol. 12, no. 2, pp. 441 - 447, 2018.
$offtext

508
509
510
* Check that we have values for the autocorrelations
$ifthen.autocorr defined p_autocorrelation

511
// Do smoothing
512
513
514
515
516
517
518
519
if(mSettings(mSolve, 'scenarios'),  // Only do smooting if using long-term scenarios
    // Select the initial sample, first `t` not in it and `f` of the last `t` in it
    loop((ms_initial(mSolve, s_), t_, ft(f_, t__))
        $[ord(t_) = msEnd(mSolve, s_) and mst_end(mSolve, s_, t__)
          and (mf_realization(mSolve, f_) xor mf_central(mSolve, f_))],
        $$batinclude 'inc/smoothing.gms' ts_influx
        $$batinclude 'inc/smoothing.gms' ts_cf
    );
520
); // END if('scenarios')
521
$endif.autocorr