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

18
* =============================================================================
19
* --- Generate model rules from basic patterns defined in the model definition files
20
* =============================================================================
21
22
// NOTE! Correctly defining multiple models still needs to be implemented!
// Pending changes?
23

24
// Initialize various sets
25
26
Option clear = t_full;
Option clear = f_solve;
27
28
Option clear = tmp;

29
30
31
32
33
// Abort model run if more than one model type is defined - unsupported at the moment
if(sum(m$mType(m), 1) > 1,
    abort "Backbone does not currently support more than one model type - you have defined more than one m";
);

34
// Loop over m
35
36
loop(m,

37
38
* --- Time Steps within Model Horizon -----------------------------------------

39
    // Determine the full set of timesteps to be considered by the defined simulation
40
    t_full(t)${ ord(t) >= mSettings(m, 't_start')
41
42
                and ord(t) <= mSettings(m, 't_end') + mSettings(m, 't_horizon')
                }
43
        = yes;
44

45
* --- Samples and Forecasts ---------------------------------------------------
46
47
48
49
50
51
$ontext
    // Check that forecast length is feasible
    if(mSettings(m, 't_forecastLength') > mSettings(m, 't_horizon'),
        abort "t_forecastLength should be less than or equal to t_horizon";
    );
$offtext
52
53
54

    // Set the time for the next available forecast.
    tForecastNext(m) = mSettings(m, 't_forecastStart');
55

56
57
58
    // Update number of samples
    if(mSettings(m, 'scenarios') = 1, mSettings(m, 'samples') = 2);

59
60
61
62
63
    // Select samples for the model
    if (not sum(s, ms(m, s)),  // unless they have been provided as input
        ms(m, s)$(ord(s) <= mSettings(m, 'samples')) = yes;
    );

64
65
66
67
68
69
70
71
72
    // Set active and previous samples
    loop(ms(m, s),
        s_active(s) = yes;
        loop(s_ $ms(m, s_),
            // Set previous samples for samples
            ss(s, s_)$(msStart(m, s) = msEnd(m, s_)) = yes;
        );
    );

73
74
75
76
77
    // Select forecasts in use for the models
    if (not sum(f, mf(m, f)),  // unless they have been provided as input
        mf(m, f)$(ord(f) <= 1 + mSettings(m, 'forecasts')) = yes;  // realization needs one f, therefore 1 + number of forecasts
    );

78
    // Select the forecasts included in the modes to be solved
79
    f_solve(f)${mf(m,f) and p_mfProbability(m, f)}
80
81
        = yes;

82
83
    // Select combinations of models, samples and forecasts to be solved
    msf(m, s, f_solve(f))$(ms(m, s) and mf(m, f)) = yes;
84
85
86
    if(mSettings(m, 'scenarios'),
        msf(ms_central(m, s), f_solve(f))$mf_realization(m, f) = no;
    );
87

88
89
    // Check the modelSolves for preset patterns for model solve timings
    // If not found, then use mSettings to set the model solve timings
90
    if(sum(modelSolves(m, t_full(t)), 1) = 0,
91
        t_skip_counter = 0;
92
        loop(t_full(t)${ ord(t) = mSettings(m, 't_start') + mSettings(m, 't_jump') * t_skip_counter
93
                        and ord(t) <= mSettings(m, 't_end')
94
                        },
95
            modelSolves(m, t) = yes;
96
97
98
99
100

            // Forecast index displacement between realized and forecasted timesteps for the initial values
            df(f, t)${  mf(m, f)
                        and ord(t) = mSettings(m, 't_start')
                        }
101
            = sum(mf_realization(m, f_), ord(f_) - ord(f));
102
103

            // Increase the t_skip counter
104
105
106
            t_skip_counter = t_skip_counter + 1;
        );
    );
107

108
* --- Intervals and Time Series -----------------------------------------------
Juha Kiviluoma's avatar
Juha Kiviluoma committed
109

110
    // Check whether the defined intervals are feasible
Juha Kiviluoma's avatar
Juha Kiviluoma committed
111
    continueLoop = 1;
112
    loop(counter${ continueLoop },
113
        if(not mInterval(m, 'lastStepInIntervalBlock', counter),
Juha Kiviluoma's avatar
Juha Kiviluoma committed
114
            continueLoop = 0;
115
        elseif mod(mInterval(m, 'lastStepInIntervalBlock', counter) - mInterval(m, 'lastStepInIntervalBlock', counter-1), mInterval(m, 'stepsPerInterval', counter)),
116
117
            put log "!!! Error occurred on interval block ", counter.tl:0 /;
            put log "!!! Abort: stepsPerInterval is not evenly divisible within the interval"
118
            abort "stepsPerInterval is not evenly divisible within the interval", m, continueLoop;
Juha Kiviluoma's avatar
Juha Kiviluoma committed
119
120
121
122
123
        else
            continueLoop = continueLoop + 1;
        );
    );

124
125
126
127
    // Determine maximum data length, if not provided in the model definition file.
    if(mSettings(m, 'dataLength'),
        tmp = max(mSettings(m, 'dataLength') + 1, tmp); // 'dataLength' increased by one to account for t000000 in ord(t)
    else
128
        put log '!!! Warning: mSettings(m, dataLength) is not defined! Calculating dataLength based on ts_influx and ts_node.' /;
129
        // Calculate the length of the time series data (based on realized forecast)
130
131
132
        option clear = tt; // Find the time steps with input time series data (ts_influx and ts_node)
        loop(gn(grid, node),
            loop(mf_realization(m, f), // Select only the realized forecast
133
                tt(t)${ts_influx(grid, node, f, t)} = yes;
134
                loop(param_gnBoundaryTypes,
135
                    tt(t)${ts_node(grid, node, param_gnBoundaryTypes, f, t)} = yes;
136
137
138
139
                ); // END loop(param_gnBoundaryTypes)
            ); // END loop(mf_realization)
        ); // END loop(gn)
        tmp = smax(tt(t), ord(t)); // Find the maximum ord(t) defined in time series data.
140
    ); // END if(mSettings(dataLength))
141

142
); // END loop(m)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
143

144
145
146
* --- Calculate Time Series Length and Circular Time Displacement -------------

// Circular displacement of time index for data loop
147
148
149
dt_circular(t_full(t))${ ord(t) > tmp }
    = - (tmp - 1) // (tmp - 1) used in order to not circulate initial values at t000000
        * floor(ord(t) / (tmp));
150

151
152
153
154
* =============================================================================
* --- Initialize Unit Efficiency Approximations -------------------------------
* =============================================================================

155
loop(m,
156
157

* --- Unit Aggregation --------------------------------------------------------
Topi Rasku's avatar
Topi Rasku committed
158

159
160
    unitAggregator_unit(unit, unit_)$sum(effLevel$(mSettingsEff(m, effLevel)), unitUnitEffLevel(unit, unit_, effLevel)) = yes;

Topi Rasku's avatar
Topi Rasku committed
161
    // Define unit aggregation sets
162
163
164
165
166
167
168
169
    unit_aggregator(unit)${ sum(unit_, unitAggregator_unit(unit, unit_)) }
        = yes; // Set of aggregator units
    unit_aggregated(unit)${ sum(unit_, unitAggregator_unit(unit_, unit)) }
        = yes; // Set of aggregated units
    unit_noAggregate(unit) = yes; // Set of units that are not aggregated into any aggregate, or are not aggregates themselves
    unit_noAggregate(unit)$unit_aggregated(unit) = no;
    unit_noAggregate(unit)${ sum((unit_, effLevel), unitUnitEffLevel(unit, unit_, effLevel)) } = no;

Topi Rasku's avatar
Topi Rasku committed
170
171
    // Process data for unit aggregations
    // Aggregate maxGen as the sum of aggregated maxGen
172
173
174
175
    p_gnu(grid, node, unit_aggregator(unit), 'maxGen')
        = sum(unit_$unitAggregator_unit(unit, unit_),
            + p_gnu(grid, node, unit_, 'maxGen')
            );
Topi Rasku's avatar
Topi Rasku committed
176
    // Aggregate maxCons as the sum of aggregated maxCons
177
178
179
180
181
    p_gnu(grid, node, unit_aggregator(unit), 'maxCons')
        = sum(unit_$unitAggregator_unit(unit, unit_),
            + p_gnu(grid, node, unit_, 'maxCons')
            );

Topi Rasku's avatar
Topi Rasku committed
182
183
* --- Calculate 'lastStepNotAggregated' for aggregated units and aggregator units ---

184
185
186
187
188
189
190
191
    loop(effLevel$mSettingsEff(m, effLevel),
        loop(effLevel_${mSettingsEff(m, effLevel_) and ord(effLevel_) < ord(effLevel)},
            p_unit(unit_aggregated(unit), 'lastStepNotAggregated')${ sum(unit_,unitUnitEffLevel(unit_, unit, effLevel)) }
                = mSettingsEff(m, effLevel_);
            p_unit(unit_aggregator(unit), 'lastStepNotAggregated')${ sum(unit_,unitUnitEffLevel(unit, unit_, effLevel)) }
                = mSettingsEff(m, effLevel_);
        );
    );
192
193
);

Topi Rasku's avatar
Topi Rasku committed
194
* --- Ensure that efficiency levels extend to the end of the model horizon and do not go beyond ---
195
196
197

loop(m,
    // First check how many efficiency levels there are and cut levels going beyond the t_horizon
198
    tmp = 0;
199
    loop(effLevel$mSettingsEff(m, effLevel),
200
201
202
203
204
        continueLoop = ord(effLevel);
        // Check if the level extends to the end of the t_horizon
        if (mSettingsEff(m, effLevel) = mSettings(m, 't_horizon'),
            tmp = 1;
        );
205
        if (mSettingsEff(m, effLevel) > mSettings(m, 't_horizon'),
206
207
208
209
210
211
212
213
214
215
            // Cut the first level going beyond the t_horizon (if the previous levels did not extend to the t_horizon)
            if (tmp = 0,
                mSettingsEff(m, effLevel) = mSettings(m, 't_horizon');
                tmp = 1;
                put log '!!! Set mSettingsEff(', m.tl:0, ', ', effLevel.tl:0, ') to ', mSettings(m, 't_horizon'):0:0 /;
            // Remove other levels going beyond the t_horizon
            else
                mSettingsEff(m, effLevel) = no;
                put log '!!! Removed mSettingsEff(', m.tl:0, ', ', effLevel.tl:0, ')' /;
            );
216
217
        );
    );
218
219
220
221
    // Ensure that that the last active level extends to the end of the t_horizon
    if ( tmp = 0,
        mSettingsEff(m, effLevel)${ord(effLevel) = continueLoop} = mSettings(m, 't_horizon');
        put log '!!! Set mSettingsEff(', m.tl:0, ', level', continueLoop, ') to ', mSettings(m, 't_horizon'):0:0 /;
222
    );
223
    // Remove effLevels with same end time step (keep the last one)
224
    loop(effLevel$mSettingsEff(m, effLevel),
225
226
        loop(effLevel_${mSettingsEff(m, effLevel_) and ord(effLevel) <> ord(effLevel_)},
            if (mSettingsEff(m, effLevel_) = mSettingsEff(m, effLevel),
227
228
                mSettingsEff(m, effLevel) = no;
                put log '!!! Removed mSettingsEff(', m.tl:0, ', ', effLevel.tl:0, ')' /;
229
            );
230
231
232
233
        );
    );
);

234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
* --- Units with online variables in the first active effLevel  ---------------

loop(m,
    continueLoop = 0;
    loop(effLevel$mSettingsEff(m, effLevel),
        continueLoop = continueLoop + 1;
        if (continueLoop = 1,
            unit_online(unit)${ sum(effSelector$effOnline(effSelector), effLevelGroupUnit(effLevel, effSelector, unit)) }
                = yes;
            unit_online_LP(unit)${ sum(effSelector, effLevelGroupUnit(effLevel, 'directOnLP', unit)) }
                = yes;
            unit_online_MIP(unit) = unit_online(unit) - unit_online_LP(unit);
        );
    );
);

250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
* --- Parse through effLevelGroupUnit and convert selected effSelectors into sets representing those selections

// Loop over effLevelGroupUnit
loop(effLevelGroupUnit(effLevel, effSelector, unit)${sum(m, mSettingsEff(m, effLevel))},

    // effSelector using DirectOff
    effGroup(effDirectOff(effSelector)) = yes;
    effGroupSelector(effDirectOff(effSelector), effSelector) = yes;
    effGroupSelectorUnit(effDirectOff(effSelector), unit, effSelector) = yes;

    // effSelector using DirectOn
    effGroup(effDirectOn(effSelector)) = yes;
    effGroupSelector(effDirectOn(effSelector), effSelector) = yes;
    effGroupSelectorUnit(effDirectOn(effSelector), unit, effSelector) = yes;

265
266
267
268
269
    // effSelector using IncHR
    effGroup(effIncHR(effSelector)) = yes;
    effGroupSelector(effIncHR(effSelector), effSelector) = yes;
    effGroupSelectorUnit(effIncHR(effSelector), unit, effSelector) = yes;

270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
    // effSelector using Lambda
    effGroup(effLambda(effSelector)) = yes;
    loop(effLambda_${ord(effLambda_) <= ord(effSelector)},
        effGroupSelector(effLambda(effSelector), effLambda_) = yes;
        effGroupSelectorUnit(effLambda(effSelector), unit, effLambda_) = yes;
        ); // END loop(effLambda_)
    ); // END loop(effLevelGroupUnit)

* --- Loop over effGroupSelectorUnit to generate efficiency approximation parameters for units

loop(effGroupSelectorUnit(effSelector, unit, effSelector_),

    // Determine the last operating point in use for the unit
    Option clear = tmp_count_op;
    loop(op${   p_unit(unit, op)    },
        tmp_count_op = ord(op);
    ); // END loop(op)

    // Parameters for direct conversion units without online variables
    if(effDirectOff(effSelector),
        p_effUnit(effSelector, unit, effSelector, 'lb') = 0; // No min load for the DirectOff approximation
291
        p_effUnit(effSelector, unit, effSelector, 'op') = smax(op, p_unit(unit, op)); // Maximum operating point
292
293
294
295
296
297
298
299
300
301
302
303
304
        p_effUnit(effSelector, unit, effSelector, 'slope') = 1 / smax(eff${p_unit(unit, eff)}, p_unit(unit, eff)); // Uses maximum found (nonzero) efficiency.
        p_effUnit(effSelector, unit, effSelector, 'section') = 0; // No section for the DirectOff approximation
    ); // END if(effDirectOff)

    // Parameters for direct conversion units with online variables
    if(effDirectOn(effSelector),
        p_effUnit(effSelector, unit, effSelector_, 'lb') = p_unit(unit, 'op00'); // op00 contains the minimum load of the unit
        p_effUnit(effSelector, unit, effSelector_, 'op') = smax(op, p_unit(unit, op)); // Maximum load determined by the largest 'op' parameter found in data
        loop(op__$(ord(op__) = tmp_count_op), // Find the maximum defined 'op'.
            loop(eff__${ord(eff__) = ord(op__)}, // ...  and the corresponding 'eff'.

                // If the minimum operating point is at zero, then the section and slope are calculated with the assumption that the efficiency curve crosses at opFirstCross
                if(p_unit(unit, 'op00') = 0,
305

306
                    // Heat rate at the cross between real efficiency curve and approximated efficiency curve
307
308
                    // !!! NOTE !!! It is advised not to define opFirstCross as any of the op points to avoid accidental division by zero!
                    heat_rate = 1 / [
309
310
311
312
313
314
315
316
317
318
319
320
321
                                    + p_unit(unit, 'eff00')
                                        * [ p_unit(unit, op__) - p_unit(unit, 'opFirstCross') ]
                                        / [ p_unit(unit, op__) - p_unit(unit, 'op00') ]
                                    + p_unit(unit, eff__)
                                        * [ p_unit(unit, 'opFirstCross') - p_unit(unit, 'op00') ]
                                        / [ p_unit(unit, op__) - p_unit(unit, 'op00') ]
                                    ];

                    // Unless section has been defined, it is calculated based on the opFirstCross
                    p_effGroupUnit(effSelector, unit, 'section') = p_unit(unit, 'section');
                    p_effGroupUnit(effSelector, unit, 'section')${ not p_effGroupUnit(effSelector, unit, 'section') }
                        = p_unit(unit, 'opFirstCross')
                            * ( heat_rate - 1 / p_unit(unit, eff__) )
322
                            / ( p_unit(unit, op__) - p_unit(unit, 'op00') );
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
                    p_effUnit(effSelector, unit, effSelector_, 'slope')
                        = 1 / p_unit(unit, eff__)
                            - p_effGroupUnit(effSelector, unit, 'section') / p_unit(unit, op__);

                // If the minimum operating point is above zero, then the approximate efficiency curve crosses the real efficiency curve at minimum and maximum.
                else
                    // Calculating the slope based on the first nonzero and the last defined data points.
                    p_effUnit(effSelector, unit, effSelector_, 'slope')
                        = (p_unit(unit, op__) / p_unit(unit, eff__) - p_unit(unit, 'op00') / p_unit(unit, 'eff00'))
                            / (p_unit(unit, op__) - p_unit(unit, 'op00'));

                    // Calculating the section based on the slope and the last defined point.
                    p_effGroupUnit(effSelector, unit, 'section')
                        = ( 1 / p_unit(unit, eff__) - p_effUnit(effSelector, unit, effSelector_, 'slope') )
                            * p_unit(unit, op__);
                ); // END if(p_unit)
            ); // END loop(eff__)
        ); // END loop(op__)
    ); // END if(effDirectOn)

    // Calculate lambdas
    if(effLambda(effSelector),
        p_effUnit(effSelector, unit, effSelector_, 'lb') = p_unit(unit, 'op00'); // op00 contains the min load of the unit

        // Calculate the relative location of the operating point in the lambdas
        tmp_op = p_unit(unit, 'op00')
                    + (ord(effSelector_)-1) / (ord(effSelector) - 1)
                        * (smax(op, p_unit(unit, op)) - p_unit(unit, 'op00'));
        p_effUnit(effSelector, unit, effSelector_, 'op') = tmp_op; // Copy the operating point to the p_effUnit

        // tmp_op falls between two p_unit defined operating points or then it is equal to one of them
        loop((op_, op__)${  (   [tmp_op > p_unit(unit, op_) and tmp_op < p_unit(unit, op__) and ord(op_) = ord(op__) - 1]
                                or [p_unit(unit, op_) = tmp_op and ord(op_) = ord(op__)]
                                )
                            and ord(op__) <= tmp_count_op
                            },
            // Find the corresponding efficiencies
            loop((eff_, eff__)${    ord(op_) = ord(eff_)
                                    and ord(op__) = ord(eff__)
                                    },
                // Calculate the distance between the operating points (zero if the points are the same)
                tmp_dist = p_unit(unit, op__) - p_unit(unit, op_);

                // If the operating points are not the same
                if (tmp_dist,
                    // Heat rate is a weighted average of the heat rates at the p_unit operating points
                    heat_rate = 1 / [
                                    + p_unit(unit, eff_) * [ p_unit(unit, op__) - tmp_op ] / tmp_dist
                                    + p_unit(unit, eff__) * [ tmp_op - p_unit(unit, op_) ] / tmp_dist
                                    ];

                // If the operating point is the same, the the heat rate can be used directly
                else
                    heat_rate = 1 / p_unit(unit, eff_);
                ); // END if(tmp_dist)

                // Special considerations for the first lambda
                if (ord(effSelector_) = 1,
                    // If the min. load of the unit is not zero or the section has been pre-defined, then section is copied directly from the unit properties
                    if(p_unit(unit, 'op00') or p_unit(unit, 'section'),
                        p_effGroupUnit(effSelector, unit, 'section') = p_unit(unit, 'section');

                    // Calculate section based on the opFirstCross, which has been calculated into p_effUnit(effLambda, unit, effLambda_, 'op')
                    else
                        p_effGroupUnit(effSelector, unit, 'section')
388
389
                            = p_unit(unit, 'opFirstCross')
                                * ( heat_rate - 1 / p_unit(unit, 'eff01') )
390
391
392
393
394
                                / ( p_unit(unit, 'op01') - tmp_op );
                    ); // END if(p_unit)
                ); // END if(ord(effSelector))

                // Calculate the slope
395
396
                p_effUnit(effSelector, unit, effSelector_, 'slope')
                    = heat_rate - p_effGroupUnit(effSelector, unit, 'section') / [tmp_op + 1${not tmp_op}];
397
398
399
            ); // END loop(eff_,eff__)
        ); // END loop(op_,op__)
    ); // END if(effLambda)
400
401
402
403
404

// Parameters for incremental heat rates
    if(effIncHR(effSelector),
        p_effUnit(effSelector, unit, effSelector, 'lb') = p_unit(unit, 'hrop00'); // hrop00 contains the minimum load of the unit
        p_effUnit(effSelector, unit, effSelector, 'op') = smax(hrop, p_unit(unit, hrop)); // Maximum operating point
405
406
        p_effUnit(effSelector, unit, effSelector, 'slope') = 1 / smax(eff${p_unit(unit, eff)}, p_unit(unit, eff)); // Uses maximum found (nonzero) efficiency.
        p_effUnit(effSelector, unit, effSelector, 'section') = p_unit(unit, 'hrsection'); // pre-defined
407
408
409
410
411
412
413
414
415
416
417
418
419

        // Whether to use q_conversionIncHR_help1 and q_conversionIncHR_help2 or not
        loop(m,
            loop(hr${p_unit(unit, hr)},
                if (mSettings(m, 'incHRAdditionalConstraints') = 0,
                    if (p_unit(unit, hr) < p_unit(unit, hr-1),
                        unit_incHRAdditionalConstraints(unit) = yes;
                    ); // END if(hr)
                else
                    unit_incHRAdditionalConstraints(unit) = yes;
                ); // END if(incHRAdditionalConstraints)
            ); // END loop(hr)
        ); // END loop(m)
420
    ); // END if(effIncHR)
421
); // END loop(effGroupSelectorUnit)
422
423

// Calculate unit wide parameters for each efficiency group
424
425
426
loop(effLevelGroupUnit(effLevel, effGroup, unit)${sum(m, mSettingsEff(m, effLevel))},
    p_effGroupUnit(effGroup, unit, 'op') = smax(effGroupSelectorUnit(effGroup, unit, effSelector), p_effUnit(effGroup, unit, effSelector, 'op'));
    p_effGroupUnit(effGroup, unit, 'lb') = smin(effGroupSelectorUnit(effGroup, unit, effSelector), p_effUnit(effGroup, unit, effSelector, 'lb'));
427
    p_effGroupUnit(effGroup, unit, 'slope') = smin(effGroupSelectorUnit(effGroup, unit, effSelector), p_effUnit(effGroup, unit, effSelector, 'slope'));
428
); // END loop(effLevelGroupUnit)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
429
430


431
432
433
434
435
* =============================================================================
* --- Initialize Unit Startup and Shutdown Counters ---------------------------
* =============================================================================

* --- Unit Start-up Generation Levels -----------------------------------------
436

437
438
loop(m,
    loop(unit$(p_unit(unit, 'rampSpeedToMinLoad') and p_unit(unit,'op00')),
439

440
        // Calculate time intervals needed for the run-up phase
441
        tmp = [ p_unit(unit,'op00') / (p_unit(unit, 'rampSpeedToMinLoad') * 60) ] / mSettings(m, 'stepLengthInHours');
442
        p_u_runUpTimeIntervals(unit) = tmp;
443
444
445
446
447
448
449
450
        p_u_runUpTimeIntervalsCeil(unit) = ceil(p_u_runUpTimeIntervals(unit));
        runUpCounter(unit, counter) // Store the required number of run-up intervals for each unit
            ${ ord(counter) <= p_u_runUpTimeIntervalsCeil(unit) }
            = yes;
        dt_trajectory(counter)
            ${ runUpCounter(unit, counter) }
            = - ord(counter) + 1; // Runup starts immediately at v_startup

451
452
        // Calculate minimum output during the run-up phase; partial intervals calculated using weighted averaging with min load
        p_uCounter_runUpMin(runUpCounter(unit, counter))
453
454
455
456
457
458
459
            = + p_unit(unit, 'rampSpeedToMinLoad')
                * ( + min(ord(counter), p_u_runUpTimeIntervals(unit)) // Location on ramp
                    - 0.5 * min(p_u_runUpTimeIntervals(unit) - ord(counter) + 1, 1) // Average ramp section
                    )
                * min(p_u_runUpTimeIntervals(unit) - ord(counter) + 1, 1) // Portion of time interval spent ramping
                * mSettings(m, 'stepLengthInHours') // Ramp length in hours
                * 60 // unit conversion from [p.u./min] to [p.u./h]
460
461
462
463
464
465
466
467
468
469
              + p_unit(unit, 'op00')${ not runUpCounter(unit, counter+1) } // Time potentially spent at min load during the last run-up interval
                * ( p_u_runUpTimeIntervalsCeil(unit) - p_u_runUpTimeIntervals(unit) );

        // Maximum output on the last run-up interval can be higher, otherwise the same as minimum.
        p_uCounter_runUpMax(runUpCounter(unit, counter))
            = p_uCounter_runUpMin(unit, counter);
        p_uCounter_runUpMax(runUpCounter(unit, counter))${ not runUpCounter(unit, counter+1) }
            = p_uCounter_runUpMax(unit, counter)
                + ( 1 - p_uCounter_runUpMax(unit, counter) )
                    * ( p_u_runUpTimeIntervalsCeil(unit) - p_u_runUpTimeIntervals(unit) );
470

471
472
473
474
475
476
477
        // Minimum ramp speed in the last interval for the run-up to min. load (p.u./min)
        p_u_minRampSpeedInLastRunUpInterval(unit)
            = p_unit(unit, 'rampSpeedToMinLoad')
                * ( p_u_runUpTimeIntervals(unit) * (p_u_runUpTimeIntervalsCeil(unit) - 0.5 * p_u_runUpTimeIntervals(unit))
                    - 0.5 * p_u_runUpTimeIntervalsCeil(unit) * p_u_runUpTimeIntervalsCeil(unit) + 1
                    );

478
    ); // END loop(unit)
479
480
); // END loop(m)

481
482
* --- Unit Shutdown Generation Levels -----------------------------------------

483
484
485
loop(m,
    loop(unit$(p_unit(unit, 'rampSpeedFromMinLoad') and p_unit(unit,'op00')),
        // Calculate time intervals needed for the shutdown phase
486
        tmp = [ p_unit(unit,'op00') / (p_unit(unit, 'rampSpeedFromMinLoad') * 60) ] / mSettings(m, 'stepLengthInHours');
487
        p_u_shutdownTimeIntervals(unit) = tmp;
488
489
        p_u_shutdownTimeIntervalsCeil(unit) = ceil(p_u_shutdownTimeIntervals(unit));
        shutdownCounter(unit, counter) // Store the required number of shutdown intervals for each unit
490
            ${ ord(counter) <= p_u_shutDownTimeIntervalsCeil(unit)}
491
492
493
494
495
            = yes;
        dt_trajectory(counter)
            ${ shutdownCounter(unit, counter) }
            = - ord(counter) + 1; // Shutdown starts immediately at v_shutdown

496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
        // Calculate minimum output during the shutdown phase; partial intervals calculated using weighted average with zero load
        p_uCounter_shutdownMin(shutdownCounter(unit, counter))
            = + p_unit(unit, 'rampSpeedFromMinLoad')
                * ( min(p_u_shutdownTimeIntervalsCeil(unit) - ord(counter) + 1, p_u_shutdownTimeIntervals(unit)) // Location on ramp
                    - 0.5 * min(p_u_shutdownTimeIntervals(unit) - p_u_shutdownTimeIntervalsCeil(unit) + ord(counter), 1) // Average ramp section
                    )
                * min(p_u_shutdownTimeIntervals(unit) - p_u_shutdownTimeIntervalsCeil(unit) + ord(counter), 1) // Portion of time interval spent ramping
                * mSettings(m, 'stepLengthInHours') // Ramp length in hours
                * 60 // unit conversion from [p.u./min] to [p.u./h]
              + p_unit(unit, 'op00')${ not shutdownCounter(unit, counter-1) } // Time potentially spent at min load on the first shutdown interval
                * ( p_u_shutdownTimeIntervalsCeil(unit) - p_u_shutdownTimeIntervals(unit) );

        // Maximum output on the first shutdown interval can be higher, otherwise the same as minimum.
        p_uCounter_shutdownMax(shutdownCounter(unit, counter))
            = p_uCounter_shutdownMin(unit, counter);
        p_uCounter_shutdownMax(shutdownCounter(unit, counter))${ not shutdownCounter(unit, counter-1) }
            = p_uCounter_shutdownMax(unit, counter)
                + ( 1 - p_uCounter_shutdownMax(unit, counter) )
                    * ( p_u_shutdownTimeIntervalsCeil(unit) - p_u_shutdownTimeIntervals(unit) );
515

516
517
518
519
520
521
522
        // Minimum ramp speed in the first interval for the shutdown from min. load (p.u./min)
        p_u_minRampSpeedInFirstShutdownInterval(unit)
            = p_unit(unit, 'rampSpeedFromMinLoad')
                * ( p_u_shutdownTimeIntervals(unit) * (p_u_shutdownTimeIntervalsCeil(unit) - 0.5 * p_u_shutdownTimeIntervals(unit))
                    - 0.5 * p_u_shutdownTimeIntervalsCeil(unit) * p_u_shutdownTimeIntervalsCeil(unit) + 1
                    );

523
    ); // END loop(unit)
524
525
); // END loop(m)

526
* --- Unit Starttype, Uptime and Downtime Counters ----------------------------
Topi Rasku's avatar
Topi Rasku committed
527

528
529
530
loop(m,
    // Loop over units with online approximations in the model
    loop(effLevelGroupUnit(effLevel, effOnline(effGroup), unit)${mSettingsEff(m, effLevel)},
531
        // Loop over the constrained start types
532
        loop(starttypeConstrained(starttype),
533
            // Find the time step displacements needed to define the start-up time frame
534
            Option clear = cc;
535
536
            cc(counter)${   ord(counter) <= p_uNonoperational(unit, starttype, 'max') / mSettings(m, 'stepLengthInHours')
                            and ord(counter) > p_uNonoperational(unit, starttype, 'min') / mSettings(m, 'stepLengthInHours')
537
538
                            }
                = yes;
Topi Rasku's avatar
Topi Rasku committed
539
            unitCounter(unit, cc(counter)) = yes;
540
541
542
            dt_starttypeUnitCounter(starttype, unit, cc(counter)) = - ord(counter);
        ); // END loop(starttypeConstrained)

543
        // Find the time step displacements needed to define the downtime requirements (include run-up phase and shutdown phase)
544
        Option clear = cc;
545
        cc(counter)${   ord(counter) <= ceil(p_unit(unit, 'minShutdownHours') / mSettings(m, 'stepLengthInHours'))
546
547
                                        + ceil(p_u_runUpTimeIntervals(unit)) // NOTE! Check this
                                        + ceil(p_u_shutdownTimeIntervals(unit)) // NOTE! Check this
548
                        }
549
            = yes;
Topi Rasku's avatar
Topi Rasku committed
550
        unitCounter(unit, cc(counter)) = yes;
551
552
553
554
        dt_downtimeUnitCounter(unit, cc(counter)) = - ord(counter);

        // Find the time step displacements needed to define the uptime requirements
        Option clear = cc;
555
        cc(counter)${ ord(counter) <= ceil(p_unit(unit, 'minOperationHours') / mSettings(m, 'stepLengthInHours'))}
556
            = yes;
Topi Rasku's avatar
Topi Rasku committed
557
        unitCounter(unit, cc(counter)) = yes;
558
559
560
561
        dt_uptimeUnitCounter(unit, cc(counter)) = - ord(counter);
    ); // END loop(effLevelGroupUnit)
); // END loop(m)

562
563
564
565
// Initialize tmp_dt based on the first model interval
loop(m,
    tmp_dt = -mInterval(m, 'stepsPerInterval', 'c000');
); // END loop(m)
566
// Estimate the maximum amount of history required for the model (very rough estimate atm, just sums all possible delays together)
567
loop(unit,
568
    tmp_dt = min(   tmp_dt, // dt operators have negative values, thus use min instead of max
569
570
571
                    smin((starttype, unitCounter(unit, counter)), dt_starttypeUnitCounter(starttype, unit, counter))
                    + smin(unitCounter(unit, counter), dt_downtimeUnitCounter(unit, counter))
                    + smin(unitCounter(unit, counter), dt_uptimeUnitCounter(unit, counter))
572
573
                    - p_u_runUpTimeIntervalsCeil(unit) // NOTE! p_u_runUpTimeIntervalsCeil is positive, whereas all dt operators are negative
                    - p_u_shutdownTimeIntervalsCeil(unit) // NOTE! p_u_shutdownTimeIntervalsCeil is positive, whereas all dt operators are negative
Topi Rasku's avatar
Topi Rasku committed
574
575
576
                    );
); // END loop(starttype, unitCounter)

577
* =============================================================================
578
579
580
581
* --- Disable reserves according to model definition --------------------------
* =============================================================================

loop(m,
582
583
584
585
586
587
588
589
590
591
592

    // Disable group reserve requirements
    restypeDirectionGroup(restype, up_down, group)
        ${  not mSettingsReservesInUse(m, restype, up_down)
            }
        = no;
    restypeDirectionGridNodeGroup(restype, up_down, grid, node, group)
        ${  not mSettingsReservesInUse(m, restype, up_down)
            }
        = no;

593
    // Disable node reserve requirements
594
    restypeDirectionGridNode(restype, up_down, grid, node)
595
596
597
598
599
        ${  not mSettingsReservesInUse(m, restype, up_down)
            }
        = no;

    // Disable node-node reserve connections
600
    restypeDirectionGridNodeNode(restype, up_down, grid, node, node_)
601
602
603
604
605
        ${  not mSettingsReservesInUse(m, restype, up_down)
            }
      = no;

    // Disable reserve provision capability from units
606
    gnuRescapable(restype, up_down, grid, node, unit)
607
608
609
610
611
612
        ${  not mSettingsReservesInUse(m, restype, up_down)
            }
      = no;
); // END loop(m)

* =============================================================================
613
614
615
* --- Various Initial Values and Calculations ---------------------------------
* =============================================================================

616
* --- Calculating fuel price time series --------------------------------------
617

618
loop(fuel${ p_fuelPrice(fuel, 'useTimeSeries') },
619
620
    // Determine the time steps where the prices change
    Option clear = tt;
621
    tt(t)${ ts_fuelPriceChange(fuel ,t) }
622
        = yes;
623
    ts_fuelPrice(fuel, t_full(t)) = sum(tt(t_)${ ord(t_) <= ord(t) }, ts_fuelPriceChange(fuel, t_));
624
); // END loop(fuel)
625

626
* --- Slack Direction ---------------------------------------------------------
627
628
629
630
631

// Upwards slack is positive, downward slack negative
p_slackDirection(upwardSlack) = 1;
p_slackDirection(downwardSlack) = -1;

632
* --- Using default value for reserves update frequency -----------------------
633

634
loop(m,
635
636
637
    p_groupReserves(group, restype, 'update_frequency')${  not p_groupReserves(group, restype, 'update_frequency')  }
        = mSettings(m, 't_jump');
    p_gnReserves(grid, node, restype, 'update_frequency')${  not p_gnReserves(grid, node, restype, 'update_frequency')  }
638
639
        = mSettings(m, 't_jump');
);
640

641
642
643
644
645
646
* --- Include 't_start' as a realized time step -------------------------------

loop(msf(m, s, f)${ mf_realization(m, f) },
    // Initial values included into previously realized time steps
    ft_realizedNoReset(f, t_full(t))${ ord(t) = mSettings(m, 't_start') }
        = yes;
Erkka Rinne's avatar
Erkka Rinne committed
647
648
    sft_realizedNoReset(s, f, t_full(t))${ ord(t) = mSettings(m, 't_start') }
        = yes;
649
650
651
652
    msft_realizedNoReset(m, s, f, t_full(t))${ ord(t) = mSettings(m, 't_start') }
        = yes;
); // END loop(m, s, f)

653
654
655
* =============================================================================
* --- Model Parameter Validity Checks -----------------------------------------
* =============================================================================
656

657
658
loop(m, // Not ideal, but multi-model functionality is not yet implemented

659
* --- Prefect foresight not longer than forecast length
660
    if(mSettings(m, 't_perfectForesight')
661
662
663
       > max(mSettings(m, 't_forecastLengthUnchanging'),
             mSettings(m, 't_forecastLengthDecreasesFrom')),
        put log "!!! Error in model ", m.tl:0 /;
664
        put log "!!! Error: t_perfectForesight > max(t_forecastLengthUnchanging, t_forecastLengthDecreasesFrom)"/;
665
666
667
        abort "Period of perfect foresight cannot be longer than forecast horizon";
    );

668
669
* --- Reserve structure checks ------------------------------------------------

670
    loop(restypeDirectionGroup(restype, up_down, group),
671
        // Check that 'update_frequency' is longer than 't_jump'
672
673
        if(p_groupReserves(group, restype, 'update_frequency') < mSettings(m, 't_jump'),
            put log '!!! Error occurred on p_groupReserves ' group.tl:0 ',' restype.tl:0 /;
674
            put log '!!! Abort: The update_frequency parameter should be longer than or equal to t_jump!' /;
675
            abort "The 'update_frequency' parameter should be longer than or equal to 't_jump'!";
676
677
678
        ); // END if('update_frequency' < 't_jump')

        // Check that 'update_frequency' is divisible by 't_jump'
679
680
        if(mod(p_groupReserves(group, restype, 'update_frequency'), mSettings(m, 't_jump')) <> 0,
            put log '!!! Error occurred on p_groupReserves ' group.tl:0 ',' restype.tl:0 /;
681
            put log '!!! Abort: The update_frequency parameter should be divisible by t_jump!' /;
682
            abort "The 'update_frequency' parameter should be divisible by 't_jump'!";
683
684
685
        ); // END if(mod('update_frequency'))

        // Check if the first interval is long enough for proper commitment of reserves
686
687
        if(mInterval(m, 'lastStepInIntervalBlock', 'c000') < p_groupReserves(group, restype, 'update_frequency') + p_groupReserves(group, restype, 'gate_closure'),
            put log '!!! Error occurred on p_groupReserves ' group.tl:0 ',' restype.tl:0 /;
Niina Helistö's avatar
Niina Helistö committed
688
689
            put log '!!! Abort: The first interval block should not be shorter than update_frequency + gate_closure for proper commitment of reserves!' /;
            abort "The first interval block should not be shorter than 'update_frequency' + 'gate_closure' for proper commitment of reserves!";
690
        ); // END if
691
    ); // END loop(restypeDirectionGroup)
692
693
694
695
696

* --- Check that there aren't more effLevels defined than exist in data -------

    if( smax(effLevel, ord(effLevel)${mSettingsEff(m, effLevel)}) > smax(effLevelGroupUnit(effLevel, effSelector, unit), ord(effLevel)),
        put log '!!! Error occurred on mSettingsEff' /;
697
698
        put log '!!! Abort: There are insufficient effLevels in the effLevelGroupUnit data for all the defined mSettingsEff!' /;
        abort "There are insufficient effLevels in the effLevelGroupUnit data for all the defined mSettingsEff!";
699
700
    ); // END if(smax)

701
702
* --- Check if time intervals are aggregated before 't_trajectoryHorizon' -----

Niina Helistö's avatar
Niina Helistö committed
703
704
    if (mInterval(m, 'lastStepInIntervalBlock', 'c000') < mSettings(m, 't_trajectoryHorizon')
       OR (mInterval(m, 'stepsPerInterval', 'c000') > 1 and mSettings(m, 't_trajectoryHorizon') > 0),
705
706
707
        put log '!!! Warning: Trajectories used on aggregated time steps! This could result in significant distortion of the trajectories.';
    ); // END if()

708
709
710
711
712
713
714
715
716
717
718
719
* --- Check that the first interval block is compatible with t_jump' ----------

    if (mod(mSettings(m, 't_jump'), mInterval(m, 'stepsPerInterval', 'c000')) <> 0,
        put log '!!! Abort: t_jump should be divisible by the first interval!' /;
        abort "'t_jump' should be divisible by the first interval!";
    ); // END if()

    if (mInterval(m, 'lastStepInIntervalBlock', 'c000') < mSettings(m, 't_jump'),
        put log '!!! Abort: The first interval block should not be shorter than t_jump!' /;
        abort "The first interval block should not be shorter than 't_jump'!";
    ); // END if()

720
); // END loop(m)
721