3c_inputsLoop.gms 28.2 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 ts_unit_update(unit, param_unit, f, t) // Update only existing values (zeroes need to be EPS)
47
                }
48
            = ts_unit_update(unit, param_unit, f, t);
49
    ); // END if('ts_unit')
50
51
52
53
$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
54
    // Update ts_effUnit
55
    if (mTimeseries_loop_read(mSolve, 'ts_effUnit'),
56
        put_utility 'gdxin' / '%input_dir%/ts_effUnit/' tSolve.tl:0 '.gdx';
57
        execute_load ts_effUnit_update=ts_effUnit;
58
        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
59
            ${  not mf_realization(mSolve, f) // Realization not updated
60
*               and ts_effUnit_update(effSelector, unit, effSelector, param_eff, f, t) // Update only existing values (zeroes need to be EPS)
61
                }
62
            = ts_effUnit_update(effSelector, unit, effSelector, param_eff, f, t);
63
    ); // END if('ts_effUnit')
64
65

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

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

    // Update ts_reserveDemand
99
    if (mTimeseries_loop_read(mSolve, 'ts_reserveDemand'),
100
        put_utility 'gdxin' / '%input_dir%/ts_reserveDemand/' tSolve.tl:0 '.gdx';
101
102
        execute_load ts_reserveDemand_update=ts_reserveDemand;
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f_solve(f), tt_forecast(t))
103
104
105
            ${  not mf_realization(mSolve, f) // Realization not updated
*                and ts_reserveDemand_update(restype, up_down, node, f, t) // Update only existing values (zeroes need to be EPS)
                }
106
107
            = ts_reserveDemand_update(restype, up_down, node, f, t);
    ); // END if('ts_reserveDemand')
108
109

    // Update ts_node
110
    if (mTimeseries_loop_read(mSolve, 'ts_node'),
111
        put_utility 'gdxin' / '%input_dir%/ts_node/' tSolve.tl:0 '.gdx';
112
113
        execute_load ts_node_update=ts_node;
        ts_node(gn(grid, node), param_gnBoundaryTypes, f_solve(f), tt_forecast(t))
114
115
116
            ${  not mf_realization(mSolve, f) // Realization not updated
*                and ts_node_update(grid, node, param_gnBoundaryTypes, f ,t) // Update only existing values (zeroes need to be EPS)
                }
117
118
119
120
121
            = 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?
122
123

    // Update ts_fuelPriceChange
124
    if (mTimeseries_loop_read(mSolve, 'ts_fuelPriceChange'),
125
        put_utility 'gdxin' / '%input_dir%/ts_fuelPriceChange/' tSolve.tl:0 '.gdx';
126
127
        execute_load ts_fuelPriceChange_update=ts_fuelPriceChange;
        ts_fuelPriceChange(fuel, tt_forecast(t))
128
*            ${ ts_fuelPriceChange_update(fuel, t) } // Update only existing values (zeroes need to be EPS)
129
130
            = ts_fuelPriceChange_update(fuel, t);
    ); // END if('ts_fuelPriceChange')
131
132

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

141
); // END if(tForecastNext)
142

143
* =============================================================================
144
* --- Optional forecast improvement -------------------------------------------
145
146
* =============================================================================

147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
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
174
175
176
177
        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
178
        // ts_effUnit
179
180
        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);
181
        // ts_effGroupUnit
182
183
        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);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
$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
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f, tt(t))
            = ts_reserveDemand(restype, up_down, node, f, t) - ts_reserveDemand(restype, up_down, node, f+ddf(f), t);
        // 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
203
        ts_unit(unit_timeseries(unit), param_unit, f, tt(t)) // Only update for units with time series enabled
204
            = [ + (ord(t) - tSolveFirst)
205
                    * ts_unit(unit, param_unit, f, t)
206
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
207
                    * ts_unit(unit, param_unit, f+ddf_(f), t)
208
                ] / mSettings(mSolve, 't_improveForecast');
209
210
$ontext
* Should these be handled here at all? See above note
211
        // ts_effUnit
212
        ts_effUnit(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f, tt(t)) // Only update for units with time series enabled
213
            = [ + (ord(t) - tSolveFirst)
214
                    * ts_effUnit(effSelector, unit, effSelector, param_eff, f, t)
215
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
216
                    * ts_effUnit(effSelector, unit, effSelector, param_eff, f+ddf_(f), t)
217
218
                ] / mSettings(mSolve, 't_improveForecast');
        // ts_effGroupUnit
219
        ts_effGroupUnit(effSelector, unit_timeseries(unit), param_eff, f, tt(t)) // Only update for units with time series enabled
220
            = [ + (ord(t) - tSolveFirst)
221
                    * ts_effGroupUnit(effSelector, unit, param_eff, f, t)
222
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
223
                    * ts_effGroupUnit(effSelector, unit, param_eff, f+ddf_(f), t)
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
                ] / 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
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f, tt(t))
            = [ + (ord(t) - tSolveFirst)
                    * ts_reserveDemand(restype, up_down, node, f, t)
                + (tSolveFirst - ord(t) + mSettings(mSolve, 't_improveForecast'))
                    * ts_reserveDemand(restype, up_down, node, f+ddf_(f), t)
                ] / 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
260
261
262
263
        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
264
        // ts_effUnit
265
266
        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);
267
        // ts_effGroupUnit
268
269
        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);
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
$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
        ts_reserveDemand(restypeDirectionNode(restype, up_down, node), f, tt(t))
            = max(ts_reserveDemand(restype, up_down, node, f, t) + ts_reserveDemand(restype, up_down, node, f+ddf(f), t), 0); // Ensure that reserve demand forecasts remains positive
        // 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)
286

287
288
289
290
291
292
293
294
295
296
* =============================================================================
* --- 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);
297
298
299
    // Make a temporary clone of tt_interval(t)
    option clear = tt_;
    tt_(tt_interval) = yes;
300
301
302
303

    // If stepsPerInterval equals one, simply use all the steps within the block
    if(mInterval(mSolve, 'stepsPerInterval', counter) = 1,

304
305
* --- Select time series data matching the intervals, for stepsPerInterval = 1, this is trivial.

306
307
308
        loop(ft(f_solve(f), tt_interval(t)),
            ts_unit_(unit_timeseries(unit), param_unit, f, t) // Only if time series enabled for the unit
                = ts_unit(unit, param_unit, f, t+dt_circular(t));
309
310
311
312
313
314
315
$ontext
* Should these be handled here at all? See above comment
            ts_effUnit_(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f_solve, t)
                = ts_effUnit(effSelector, unit, effSelector, param_eff, f_solve, t+dt_circular(t));
            ts_effGroupUnit_(effSelector, unit_timeseries(unit), param_eff, f_solve, t)
                = ts_effGroupUnit(effSelector, unit, param_eff, f_solve, t+dt_circular(t));
$offtext
316
317
318
319
320
321
            ts_cf_(flowNode(flow, node), f, t, s)$sft(s, f, t)
                = ts_cf(flow, node, f_solve, t + (dt_scenarioOffset(flow, node, 'ts_cf', s)
                                                  + dt_circular(t)$(not gn_scenarios(flow, node, 'ts_cf'))));
            ts_influx_(gn(grid, node), f, t, s)$sft(s, f, t)
                = ts_influx(grid, node, f, t + (dt_scenarioOffset(grid, node, 'ts_influx', s)
                                                + dt_circular(t)$(not gn_scenarios(grid, node, 'ts_influx'))));
322
323
324
325
326
            // Reserve demand relevant only up until reserve_length
            ts_reserveDemand_(restypeDirectionNode(restype, up_down, node), f_solve, t)
              ${ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')}
                = ts_reserveDemand(restype, up_down, node, f_solve, t+dt_circular(t));
            ts_node_(gn_state(grid, node), param_gnBoundaryTypes, f_solve, t, s)
327
328
329
330
              $p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')

                = ts_node(grid, node, param_gnBoundaryTypes, f_solve, t + (dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                                                                           + dt_circular(t)$(not gn_scenarios(grid, node, 'ts_node'))));
331
332
333
334
335
            // Fuel price time series
            ts_fuelPrice_(fuel, t)
                = ts_fuelPrice(fuel, t+dt_circular(t));
        ); // END loop(ft)

336
337
* --- If stepsPerInterval exceeds 1 (stepsPerInterval < 1 not defined) --------

338
339
340
341
    elseif mInterval(mSolve, 'stepsPerInterval', counter) > 1,

        // Select and average time series data matching the intervals, for stepsPerInterval > 1
        // Loop over the t:s of the interval
342
        loop(ft(f_solve(f), tt_interval(t)),
343
344
            // Select t:s within the interval
            Option clear = tt;
345
346
347
            tt(tt_(t_))
                ${  ord(t_) >= ord(t)
                    and ord(t_) < ord(t) + mInterval(mSolve, 'stepsPerInterval', counter)
348
349
                 }
                = yes;
350
351
352

            ts_unit_(unit_timeseries(unit), param_unit, f, t)
                = sum(tt(t_), ts_unit(unit, param_unit, f, t_+dt_circular(t_)))
353
354
355
356
357
358
359
360
361
362
                    / mInterval(mSolve, 'stepsPerInterval', counter);
$ontext
* Should these be handled here at all? See above comment
            ts_effUnit_(effGroupSelectorUnit(effSelector, unit_timeseries(unit), effSelector), param_eff, f_solve, t)
                = sum(tt(t_), ts_effUnit(effSelector, unit, effSelector, param_eff, f_solve, t_+dt_circular(t_))))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            ts_effGroupUnit_(effSelector, unit_timeseries(unit), param_eff, f_solve, t)
                = sum(tt(t_), ts_effGroupUnit(effSelector, unit, param_eff, f_solve, t_+dt_circular(t_))))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
$offtext
363
364
365
            ts_influx_(gn(grid, node), f_solve, t, s)$sft(s, f, t)
                = sum(tt(t_), ts_influx(grid, node, f_solve, t_ + (dt_scenarioOffset(grid, node, 'ts_influx', s)
                                                                   + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_influx')))))
366
                    / mInterval(mSolve, 'stepsPerInterval', counter);
367
368
369
            ts_cf_(flowNode(flow, node), f_solve, t, s)$sft(s, f, t)
                = sum(tt(t_), ts_cf(flow, node, f_solve, t_ + (dt_scenarioOffset(flow, node, 'ts_cf', s)
                                                               + dt_circular(t_)$(not gn_scenarios(flow, node, 'ts_cf')))))
370
371
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            // Reserves relevant only until reserve_length
372
            ts_reserveDemand_(restypeDirectionNode(restype, up_down, node), f, t)
373
374
375
              ${ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')  }
                = sum(tt(t_), ts_reserveDemand(restype, up_down, node, f_solve, t_+dt_circular(t_)))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
376
            ts_node_(gn_state(grid, node), param_gnBoundaryTypes, f, t, s)
377
              ${p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')
378
                and sft(s, f, t)}
379
                   // Take average if not a limit type
380
381
                = (sum(tt(t_), ts_node(grid, node, param_gnBoundaryTypes, f_solve, t_ + (dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                                                                                         + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node')))))
382
                    / mInterval(mSolve, 'stepsPerInterval', counter))$(not (sameas(param_gnBoundaryTypes, 'upwardLimit') or sameas(param_gnBoundaryTypes, 'downwardLimit')))
383
                  // Maximum lower limit
384
385
                  + smax(tt(t_), ts_node(grid, node, param_gnBoundaryTypes, f_solve, t_ + (dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                                                                                           + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node')))))
386
387
                      $sameas(param_gnBoundaryTypes, 'downwardLimit')
                  // Minimum upper limit
388
389
                  + smin(tt(t_), ts_node(grid, node, param_gnBoundaryTypes, f_solve, t_ + (dt_scenarioOffset(grid, node, param_gnBoundaryTypes, s)
                                                                                           + dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_node')))))
390
391
392
393
394
395
396
397
398
399
                       $sameas(param_gnBoundaryTypes, 'upwardLimit');
            // Fuel price time series
            ts_fuelPrice_(fuel, t)
                = sum(tt(t_), ts_fuelPrice(fuel, t_+dt_circular(t_)))
                    / mInterval(mSolve, 'stepsPerInterval', counter);
            ); // END loop(ft)

    ); // END if(stepsPerInterval)
); // END loop(counter)

400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
* --- 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
424
425
426
427
428
429
430
431
432
* =============================================================================
* --- Input data processing ---------------------------------------------------
* =============================================================================

* --- Scenario reduction ------------------------------------------------------
if(active(mSolve, 'scenred'),
    $$include 'inc/scenred.gms'
);

433
434
435
436
437
438
439
440
* --- Update probabilities ----------------------------------------------------
Option clear = p_msft_probability;
p_msft_probability(msft(mSolve, s, f, t))
    = p_mfProbability(mSolve, f)
        / sum(f_${ft(f_, t)},
              p_mfProbability(mSolve, f_)) * p_msProbability(mSolve, s);


441
442
443
444
445
446
447
448
449
450
* --- Calculate sample displacements ------------------------------------------
Options clear = ds, clear = ds_state;
loop(ss(s, s_),
    ds(s, t)$(ord(t) = tSolveFirst + msStart(mSolve, s)) = -(ord(s) - ord(s_));
    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
451
452
453
454
455
456
457
458
459
460
461
* --- Smooting of stochastic scenarios ----------------------------------------
$ontext
First calculate standard deviation for over all samples, then smoothen the scenarios
following the methodology presented in [1, p. 443]. This avoids a discontinuity
`jump' after the initial sample.

[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

462
463
* Only do smoothing if scenarios are used
if(mSettings(msolve, 'scenarios'),
464

Erkka Rinne's avatar
Erkka Rinne committed
465
466
* Influx
loop(gn(grid, node)$p_autocorrelation(grid, node, 'ts_influx'),
467
468
469
    ts_influx_mean(grid, node, ft(f, t))$(mf_central(mSolve, f) and sum(sft(s, f, t), 1))
        = sum(sft(s, f, t), ts_influx_(grid, node, f, t, s))
                / sum(sft(s, f, t), 1);
Erkka Rinne's avatar
Erkka Rinne committed
470

471
472
    ts_influx_std(grid, node, ft(f, t))$(mf_central(mSolve, f) and sum(sft(s, f, t), 1))
        = sqrt(sum(sft(s, f, t), sqr(ts_influx_(grid, node, f, t, s)
Erkka Rinne's avatar
Erkka Rinne committed
473
                                         - ts_influx_mean(grid, node, f, t)))
474
                / sum(sft(s, f, t), 1)
Erkka Rinne's avatar
Erkka Rinne committed
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
          );

    // Do smoothing
    loop(mst_end(ms_initial(mSolve, s_), t_),
        ts_influx_(grid, node, ft(f, t), s)$(ts_influx_std(grid, node, f, t_+dt_circular(t_))
                                             and sft(s, f, t)
                                             and not ms_initial(mSolve, s))
            = min(p_tsMaxValue(node, 'ts_influx'), max(p_tsMinValue(node, 'ts_influx'),
              ts_influx_(grid, node, f, t, s)
              + (ts_influx_(grid, node, f, t_, s_)
                 - ts_influx_(grid, node, f, t_, s))
                * (ts_influx_std(grid, node, f, t+dt_circular(t))
                    / ts_influx_std(grid, node, f, t_+dt_circular(t_)))
                * power(p_autocorrelation(grid, node, 'ts_influx'), abs(ord(t) - ord(t_)))
              ));
    );
);

* CF
loop(flowNode(flow, node)$p_autocorrelation(flow, node, 'ts_cf'),
495
496
497
    ts_cf_mean(flow, node, ft(f, t))$(mf_central(mSolve, f) and sum(sft(s, f, t), 1))
        = sum(sft(s, f, t), ts_cf_(flow, node, f, t, s))
                / sum(sft(s, f, t), 1);
Erkka Rinne's avatar
Erkka Rinne committed
498

499
500
    ts_cf_std(flow, node, ft(f, t))$(mf_central(mSolve, f) and sum(sft(s, f, t), 1))
        = sqrt(sum(sft(s, f, t), sqr(ts_cf_(flow, node, f, t, s)
Erkka Rinne's avatar
Erkka Rinne committed
501
                                         - ts_cf_mean(flow, node, f, t)))
502
                / sum(sft(s, f, t), 1)
Erkka Rinne's avatar
Erkka Rinne committed
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
          );

    // Do smoothing
    loop(mst_end(ms_initial(mSolve, s_), t_),
        ts_cf_(flow, node, ft(f, t), s)$(ts_cf_std(flow, node, f, t_+dt_circular(t_))
                                         and sft(s, f, t)
                                         and not ms_initial(mSolve, s))
            = min(p_tsMaxValue(node, 'ts_cf'), max(p_tsMinValue(node, 'ts_cf'),
              ts_cf_(flow, node, f, t, s)
              + (ts_cf_(flow, node, f, t_, s_)
                 - ts_cf_(flow, node, f, t_, s))
                * (ts_cf_std(flow, node, f, t+dt_circular(t))
                    / ts_cf_std(flow, node, f, t_+dt_circular(t_)))
                * power(p_autocorrelation(flow, node, 'ts_cf'), abs(ord(t) - ord(t_)))
              ));
    );
);
520

521
); // end if(mSettings(msolve, 'scenarios'),