3c_inputsLoop.gms 26.9 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
130
131
132
133
134
135
    // Update ts_priceChange
    if (mTimeseries_loop_read(mSolve, 'ts_priceChange'),
        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')
136
137

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

146
); // END if(tForecastNext)
147

148
* =============================================================================
149
* --- Optional forecast improvement -------------------------------------------
150
151
* =============================================================================

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
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
179
180
181
182
        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
183
        // ts_effUnit
184
185
        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);
186
        // ts_effGroupUnit
187
188
        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);
189
190
191
192
193
194
195
196
$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
197
198
        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);
199
200
201
202
203
204
205
206
207
        // 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
208
        ts_unit(unit_timeseries(unit), param_unit, f, tt(t)) // Only update for units with time series enabled
209
            = [ + (ord(t) - tSolveFirst)
210
                    * ts_unit(unit, param_unit, f, t)
211
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
212
                    * ts_unit(unit, param_unit, f+ddf_(f), t)
213
                ] / mSettings(mSolve, 't_improveForecast');
214
215
$ontext
* Should these be handled here at all? See above note
216
        // ts_effUnit
217
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f, tt(t)) // Only update for units with time series enabled
218
            = [ + (ord(t) - tSolveFirst)
219
                    * ts_effUnit(effSelector, unit, effSelector, param_eff, f, t)
220
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
221
                    * ts_effUnit(effSelector, unit, effSelector, param_eff, f+ddf_(f), t)
222
223
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_effGroupUnit
224
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f, tt(t)) // Only update for units with time series enabled
225
            = [ + (ord(t) - tSolveFirst)
226
                    * ts_effGroupUnit(effSelector, unit, param_eff, f, t)
227
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
228
                    * ts_effGroupUnit(effSelector, unit, param_eff, f+ddf_(f), t)
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
                ] / 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
246
        ts_reserveDemand(restypeDirectionGroup(restype, up_down, group), f, tt(t))
247
            = [ + (ord(t) - tSolveFirst)
248
                    * ts_reserveDemand(restype, up_down, group, f, t)
249
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
250
                    * ts_reserveDemand(restype, up_down, group, f+ddf_(f), t)
251
252
253
254
255
256
257
258
259
260
261
262
263
264
                ] / 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
265
266
267
268
        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
269
        // ts_effUnit
270
271
        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);
272
        // ts_effGroupUnit
273
274
        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);
275
276
277
278
279
280
281
282
$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
283
284
        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
285
286
287
288
289
290
        // 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)
291

292
293
294
295
296
297
298
299
300
301
302
* =============================================================================
* --- 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);

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

381
382
383
    // Commodity price time series
    ts_vomCost_(gnu(grid, node, unit), tt_interval(t))
        = + p_gnu(grid, node, unit, 'vomCosts')
384
          // input commodity cost
385
          + sum(gnu_input(grid, commodity, unit)$un_commodity_in(unit,commodity), 
386
387
388
389
390
391
              + p_price(commodity, 'price')$p_price(commodity, 'useConstant')
              + sum(tt_aggregate(t, t_)$p_price(commodity, 'useTimeSeries'),
                  + ts_price(node, t_+dt_circular(t_))
                )
                / mInterval(mSolve, 'stepsPerInterval', counter)
            )
392
          // output commodity cost
393
          - sum(gnu_output(grid, commodity, unit)$un_commodity_out(unit,commodity),  
394
395
396
397
398
399
400
              + p_price(commodity, 'price')$p_price(commodity, 'useConstant')
              + sum(tt_aggregate(t, t_)$p_price(commodity, 'useTimeSeries'),
                  + ts_price(node, t_+dt_circular(t_))
                )
                / mInterval(mSolve, 'stepsPerInterval', counter)
            )
          // emission cost
401
402
403
404
          + sum(emission$p_unitEmissionCost(unit, node, emission), // Emission taxes
              + p_unitEmissionCost(unit, node, emission)
            ); // END sum(emission)

405
    p_unStartup(unit, node, starttype)$p_uStartupfuel(unit, node, 'fixedFuelFraction')
406
      =
407
408
        + p_uStartup(unit, starttype, 'consumption')
            * p_uStartupfuel(unit, node, 'fixedFuelFraction');
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426

    // Calculating startup cost time series
    ts_startupCost_(unit, starttype, tt_interval(t))
      =
        + p_uStartup(unit, starttype, 'cost')
        // Start-up fuel and emission costs
        + sum(commodity$p_uStartupfuel(unit, commodity, 'fixedFuelFraction'),
            + p_uStartup(unit, starttype, 'consumption')
              * p_uStartupfuel(unit, commodity, 'fixedFuelFraction')
              * [
                  + p_price(commodity, 'price')$p_price(commodity, 'useConstant')
                  + sum(tt_aggregate(t, t_)$p_price(commodity, 'useTimeseries'),
                      + ts_price(commodity, t_+dt_circular(t_))
                    )
                    / mInterval(mSolve, 'stepsPerInterval', counter)
                ] // END * p_uStartup
          ) // END sum(commodity)
        + sum((node, emission)$p_unitEmissionCost(unit, node, emission),
427
            + p_unStartup(unit, node, starttype)
428
              * p_unitEmissionCost(unit, node, emission)
429
          ) // END sum(node, emission)
430
431
); // END loop(counter)

432

433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
* --- 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)

457

Erkka Rinne's avatar
Erkka Rinne committed
458
459
460
461
* =============================================================================
* --- Input data processing ---------------------------------------------------
* =============================================================================

Erkka Rinne's avatar
Erkka Rinne committed
462
$ifthen.scenarios defined scenario
Erkka Rinne's avatar
Erkka Rinne committed
463
* --- Scenario reduction ------------------------------------------------------
464
if(active(mSolve, 'scenred') and mSettings('schedule', 'scenarios') > 1,
Erkka Rinne's avatar
Erkka Rinne committed
465
466
    $$include 'inc/scenred.gms'
);
Erkka Rinne's avatar
Erkka Rinne committed
467
$endif.scenarios
Erkka Rinne's avatar
Erkka Rinne committed
468

469
470
471
472
* --- Update probabilities ----------------------------------------------------
Option clear = p_msft_probability;
p_msft_probability(msft(mSolve, s, f, t))
    = p_mfProbability(mSolve, f)
473
474
475
476
        / sum(f_$ft(f_, t),
              p_mfProbability(mSolve, f_)
          ) * p_msProbability(mSolve, s)
            * p_msWeight(mSolve, s);
477
478


479
480
* --- Calculate sample displacements ------------------------------------------
Options clear = ds, clear = ds_state;
Erkka Rinne's avatar
Erkka Rinne committed
481
482
loop((mst_start(mSolve, s, t), ss(s, s_)),
    ds(s, t) = -(ord(s) - ord(s_));
483
484
485
486
487
488
    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
489
490
* --- Smooting of stochastic scenarios ----------------------------------------
$ontext
491
492
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
493
494
495
496
497
498

[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

499
500
501
* Check that we have values for the autocorrelations
$ifthen.autocorr defined p_autocorrelation

502
// Do smoothing
503
504
505
506
507
508
509
510
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
    );
511
); // END if('scenarios')
512
$endif.autocorr