3c_inputsLoop.gms 24.8 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

    // Update ts_fuelPriceChange
129
    if (mTimeseries_loop_read(mSolve, 'ts_fuelPriceChange'),
130
        put_utility 'gdxin' / '%input_dir%/ts_fuelPriceChange/' tSolve.tl:0 '.gdx';
131
132
        execute_load ts_fuelPriceChange_update=ts_fuelPriceChange;
        ts_fuelPriceChange(fuel, tt_forecast(t))
133
*            ${ ts_fuelPriceChange_update(fuel, t) } // Update only existing values (zeroes need to be EPS)
134
135
            = ts_fuelPriceChange_update(fuel, t);
    ); // END if('ts_fuelPriceChange')
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
381
382
                )
                $(sameas(param_gnBoundaryTypes, 'upwardLimit') or upwardSlack(param_gnBoundaryTypes));
    // Fuel price time series
    ts_fuelPrice_(fuel, tt_interval(t))
        ${ p_fuelPrice(fuel, 'useTimeSeries') }
383
384
        = sum(tt_aggregate(t, t_),
            + ts_fuelPrice(fuel, t_+dt_circular(t_))
385
386
387
            )
            / mInterval(mSolve, 'stepsPerInterval', counter);

388
389
); // END loop(counter)

390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
* --- 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)

Erkka Rinne's avatar
Erkka Rinne committed
414
415
416
417
* =============================================================================
* --- Input data processing ---------------------------------------------------
* =============================================================================

Erkka Rinne's avatar
Erkka Rinne committed
418
$ifthen.scenarios defined scenario
Erkka Rinne's avatar
Erkka Rinne committed
419
* --- Scenario reduction ------------------------------------------------------
420
if(active(mSolve, 'scenred') and mSettings('schedule', 'scenarios') > 1,
Erkka Rinne's avatar
Erkka Rinne committed
421
422
    $$include 'inc/scenred.gms'
);
Erkka Rinne's avatar
Erkka Rinne committed
423
$endif.scenarios
Erkka Rinne's avatar
Erkka Rinne committed
424

425
426
427
428
429
* --- Update probabilities ----------------------------------------------------
Option clear = p_msft_probability;
p_msft_probability(msft(mSolve, s, f, t))
    = p_mfProbability(mSolve, f)
        / sum(f_${ft(f_, t)},
430
431
              p_mfProbability(mSolve, f_)) * p_msProbability(mSolve, s)
              * p_msWeight(mSolve, s);
432
433


434
435
* --- Calculate sample displacements ------------------------------------------
Options clear = ds, clear = ds_state;
Erkka Rinne's avatar
Erkka Rinne committed
436
437
loop((mst_start(mSolve, s, t), ss(s, s_)),
    ds(s, t) = -(ord(s) - ord(s_));
438
439
440
441
442
443
    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
444
445
* --- Smooting of stochastic scenarios ----------------------------------------
$ontext
446
447
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
448
449
450
451
452
453

[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

454
455
456
* Check that we have values for the autocorrelations
$ifthen.autocorr defined p_autocorrelation

457
// Do smoothing
458
459
460
461
462
463
464
465
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
    );
466
); // END if('scenarios')
467
$endif.autocorr