3a_periodicInit.gms 35.3 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
59
60
    // 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;
    );

61
62
63
64
65
66
67
68
69
    // 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;
        );
    );

70
71
72
73
74
    // 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
    );

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

79
80
81
    // Select combinations of models, samples and forecasts to be solved
    msf(m, s, f_solve(f))$(ms(m, s) and mf(m, f)) = yes;

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

            // 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')
                        }
95
            = sum(mf_realization(m, f_), ord(f_) - ord(f));
96
97

            // Increase the t_skip counter
98
99
100
            t_skip_counter = t_skip_counter + 1;
        );
    );
101

102
* --- Intervals and Time Series -----------------------------------------------
Juha Kiviluoma's avatar
Juha Kiviluoma committed
103

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

118
119
120
121
    // 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
122
        put log '!!! Warning: mSettings(m, dataLength) is not defined! Calculating dataLength based on ts_influx and ts_node.' /;
123
        // Calculate the length of the time series data (based on realized forecast)
124
125
126
        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
127
                tt(t)${ts_influx(grid, node, f, t)} = yes;
128
                loop(param_gnBoundaryTypes,
129
                    tt(t)${ts_node(grid, node, param_gnBoundaryTypes, f, t)} = yes;
130
131
132
133
                ); // 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.
134
    ); // END if(mSettings(dataLength))
135

136
); // END loop(m)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
137

138
139
140
* --- Calculate Time Series Length and Circular Time Displacement -------------

// Circular displacement of time index for data loop
141
142
143
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));
144

145
146
147
148
* =============================================================================
* --- Initialize Unit Efficiency Approximations -------------------------------
* =============================================================================

149
loop(m,
150
151

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

153
154
    unitAggregator_unit(unit, unit_)$sum(effLevel$(mSettingsEff(m, effLevel)), unitUnitEffLevel(unit, unit_, effLevel)) = yes;

Topi Rasku's avatar
Topi Rasku committed
155
    // Define unit aggregation sets
156
157
158
159
160
161
162
163
    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
164
165
    // Process data for unit aggregations
    // Aggregate maxGen as the sum of aggregated maxGen
166
167
168
169
    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
170
    // Aggregate maxCons as the sum of aggregated maxCons
171
172
173
174
175
    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
176
177
* --- Calculate 'lastStepNotAggregated' for aggregated units and aggregator units ---

178
179
180
181
182
183
184
185
    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_);
        );
    );
186
187
);

Topi Rasku's avatar
Topi Rasku committed
188
* --- Ensure that efficiency levels extend to the end of the model horizon and do not go beyond ---
189
190
191

loop(m,
    // First check how many efficiency levels there are and cut levels going beyond the t_horizon
192
    tmp = 0;
193
    loop(effLevel$mSettingsEff(m, effLevel),
194
195
196
197
198
        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;
        );
199
        if (mSettingsEff(m, effLevel) > mSettings(m, 't_horizon'),
200
201
202
203
204
205
206
207
208
209
            // 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, ')' /;
            );
210
211
        );
    );
212
213
214
215
    // 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 /;
216
    );
217
    // Remove effLevels with same end time step (keep the last one)
218
    loop(effLevel$mSettingsEff(m, effLevel),
219
220
        loop(effLevel_${mSettingsEff(m, effLevel_) and ord(effLevel) <> ord(effLevel_)},
            if (mSettingsEff(m, effLevel_) = mSettingsEff(m, effLevel),
221
222
                mSettingsEff(m, effLevel) = no;
                put log '!!! Removed mSettingsEff(', m.tl:0, ', ', effLevel.tl:0, ')' /;
223
            );
224
225
226
227
        );
    );
);

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
* --- 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);
        );
    );
);

244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
* --- 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;

259
260
261
262
263
    // effSelector using IncHR
    effGroup(effIncHR(effSelector)) = yes;
    effGroupSelector(effIncHR(effSelector), effSelector) = yes;
    effGroupSelectorUnit(effIncHR(effSelector), unit, effSelector) = yes;

264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
    // 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
285
        p_effUnit(effSelector, unit, effSelector, 'op') = smax(op, p_unit(unit, op)); // Maximum operating point
286
287
288
289
290
291
292
293
294
295
296
297
298
        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,
299

300
                    // Heat rate at the cross between real efficiency curve and approximated efficiency curve
301
302
                    // !!! NOTE !!! It is advised not to define opFirstCross as any of the op points to avoid accidental division by zero!
                    heat_rate = 1 / [
303
304
305
306
307
308
309
310
311
312
313
314
315
                                    + 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__) )
316
                            / ( p_unit(unit, op__) - p_unit(unit, 'op00') );
317
318
319
320
321
322
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
                    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')
382
383
                            = p_unit(unit, 'opFirstCross')
                                * ( heat_rate - 1 / p_unit(unit, 'eff01') )
384
385
386
387
388
                                / ( p_unit(unit, 'op01') - tmp_op );
                    ); // END if(p_unit)
                ); // END if(ord(effSelector))

                // Calculate the slope
389
390
                p_effUnit(effSelector, unit, effSelector_, 'slope')
                    = heat_rate - p_effGroupUnit(effSelector, unit, 'section') / [tmp_op + 1${not tmp_op}];
391
392
393
            ); // END loop(eff_,eff__)
        ); // END loop(op_,op__)
    ); // END if(effLambda)
394
395
396
397
398

// 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
399
400
        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
401
402
403
404
405
406
407
408
409
410
411
412
413

        // 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)
414
    ); // END if(effIncHR)
415
); // END loop(effGroupSelectorUnit)
416
417

// Calculate unit wide parameters for each efficiency group
418
419
420
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'));
421
    p_effGroupUnit(effGroup, unit, 'slope') = smin(effGroupSelectorUnit(effGroup, unit, effSelector), p_effUnit(effGroup, unit, effSelector, 'slope'));
422
); // END loop(effLevelGroupUnit)
Juha Kiviluoma's avatar
Juha Kiviluoma committed
423
424


425
426
427
428
429
* =============================================================================
* --- Initialize Unit Startup and Shutdown Counters ---------------------------
* =============================================================================

* --- Unit Start-up Generation Levels -----------------------------------------
430

431
432
loop(m,
    loop(unit$(p_unit(unit, 'rampSpeedToMinLoad') and p_unit(unit,'op00')),
433

434
        // Calculate time intervals needed for the run-up phase
435
        tmp = [ p_unit(unit,'op00') / (p_unit(unit, 'rampSpeedToMinLoad') * 60) ] / mSettings(m, 'stepLengthInHours');
436
        p_u_runUpTimeIntervals(unit) = tmp;
437
438
439
440
441
442
443
444
        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

445
446
        // Calculate minimum output during the run-up phase; partial intervals calculated using weighted averaging with min load
        p_uCounter_runUpMin(runUpCounter(unit, counter))
447
448
449
450
451
452
453
            = + 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]
454
455
456
457
458
459
460
461
462
463
              + 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) );
464

465
    ); // END loop(unit)
466
467
); // END loop(m)

468
469
* --- Unit Shutdown Generation Levels -----------------------------------------

470
471
472
loop(m,
    loop(unit$(p_unit(unit, 'rampSpeedFromMinLoad') and p_unit(unit,'op00')),
        // Calculate time intervals needed for the shutdown phase
473
        tmp = [ p_unit(unit,'op00') / (p_unit(unit, 'rampSpeedFromMinLoad') * 60) ] / mSettings(m, 'stepLengthInHours');
474
        p_u_shutdownTimeIntervals(unit) = tmp;
475
476
        p_u_shutdownTimeIntervalsCeil(unit) = ceil(p_u_shutdownTimeIntervals(unit));
        shutdownCounter(unit, counter) // Store the required number of shutdown intervals for each unit
477
            ${ ord(counter) <= p_u_shutDownTimeIntervalsCeil(unit)}
478
479
480
481
482
            = yes;
        dt_trajectory(counter)
            ${ shutdownCounter(unit, counter) }
            = - ord(counter) + 1; // Shutdown starts immediately at v_shutdown

483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
        // 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) );
502

503
    ); // END loop(unit)
504
505
); // END loop(m)

506
* --- Unit Starttype, Uptime and Downtime Counters ----------------------------
Topi Rasku's avatar
Topi Rasku committed
507

508
509
510
loop(m,
    // Loop over units with online approximations in the model
    loop(effLevelGroupUnit(effLevel, effOnline(effGroup), unit)${mSettingsEff(m, effLevel)},
511
        // Loop over the constrained start types
512
        loop(starttypeConstrained(starttype),
513
            // Find the time step displacements needed to define the start-up time frame
514
            Option clear = cc;
515
516
            cc(counter)${   ord(counter) <= p_uNonoperational(unit, starttype, 'max') / mSettings(m, 'stepLengthInHours')
                            and ord(counter) > p_uNonoperational(unit, starttype, 'min') / mSettings(m, 'stepLengthInHours')
517
518
                            }
                = yes;
Topi Rasku's avatar
Topi Rasku committed
519
            unitCounter(unit, cc(counter)) = yes;
520
521
522
            dt_starttypeUnitCounter(starttype, unit, cc(counter)) = - ord(counter);
        ); // END loop(starttypeConstrained)

523
        // Find the time step displacements needed to define the downtime requirements (include run-up phase and shutdown phase)
524
        Option clear = cc;
525
        cc(counter)${   ord(counter) <= ceil(p_unit(unit, 'minShutdownHours') / mSettings(m, 'stepLengthInHours'))
526
527
                                        + ceil(p_u_runUpTimeIntervals(unit)) // NOTE! Check this
                                        + ceil(p_u_shutdownTimeIntervals(unit)) // NOTE! Check this
528
                        }
529
            = yes;
Topi Rasku's avatar
Topi Rasku committed
530
        unitCounter(unit, cc(counter)) = yes;
531
532
533
534
        dt_downtimeUnitCounter(unit, cc(counter)) = - ord(counter);

        // Find the time step displacements needed to define the uptime requirements
        Option clear = cc;
535
        cc(counter)${ ord(counter) <= ceil(p_unit(unit, 'minOperationHours') / mSettings(m, 'stepLengthInHours'))}
536
            = yes;
Topi Rasku's avatar
Topi Rasku committed
537
        unitCounter(unit, cc(counter)) = yes;
538
539
540
541
        dt_uptimeUnitCounter(unit, cc(counter)) = - ord(counter);
    ); // END loop(effLevelGroupUnit)
); // END loop(m)

542
543
544
545
// Initialize tmp_dt based on the first model interval
loop(m,
    tmp_dt = -mInterval(m, 'stepsPerInterval', 'c000');
); // END loop(m)
546
// Estimate the maximum amount of history required for the model (very rough estimate atm, just sums all possible delays together)
547
loop(unit,
548
    tmp_dt = min(   tmp_dt, // dt operators have negative values, thus use min instead of max
549
550
551
                    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))
552
553
                    - 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
554
555
556
                    );
); // END loop(starttype, unitCounter)

557
* =============================================================================
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
* --- Disable reserves according to model definition --------------------------
* =============================================================================

loop(m,
    // Disable node reserve requirements
    restypeDirectionNode(restype, up_down, node)
        ${  not mSettingsReservesInUse(m, restype, up_down)
            }
        = no;

    // Disable node-node reserve connections
    restypeDirectionNodeNode(restype, up_down, node, node_)
        ${  not mSettingsReservesInUse(m, restype, up_down)
            }
      = no;

    // Disable reserve provision capability from units
    nuRescapable(restype, up_down, node, unit)
        ${  not mSettingsReservesInUse(m, restype, up_down)
            }
      = no;
); // END loop(m)

* =============================================================================
582
583
584
* --- Various Initial Values and Calculations ---------------------------------
* =============================================================================

585
* --- Calculating fuel price time series --------------------------------------
586

587
loop(fuel${ p_fuelPrice(fuel, 'useTimeSeries') },
588
589
    // Determine the time steps where the prices change
    Option clear = tt;
590
    tt(t)${ ts_fuelPriceChange(fuel ,t) }
591
        = yes;
592
    ts_fuelPrice(fuel, t_full(t)) = sum(tt(t_)${ ord(t_) <= ord(t) }, ts_fuelPriceChange(fuel, t_));
593
); // END loop(fuel)
594

595
* --- Slack Direction ---------------------------------------------------------
596
597
598
599
600

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

601
* --- Using default value for reserves update frequency -----------------------
602

603
604
605
606
loop(m,
    p_nReserves(node, restype, 'update_frequency')${  not p_nReserves(node, restype, 'update_frequency')  }
        = mSettings(m, 't_jump');
);
607

608
609
610
611
612
613
* --- 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
614
615
    sft_realizedNoReset(s, f, t_full(t))${ ord(t) = mSettings(m, 't_start') }
        = yes;
616
617
618
619
    msft_realizedNoReset(m, s, f, t_full(t))${ ord(t) = mSettings(m, 't_start') }
        = yes;
); // END loop(m, s, f)

620
621
622
* =============================================================================
* --- Model Parameter Validity Checks -----------------------------------------
* =============================================================================
623

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

626
627
628
629
630
631
632
633
634
* --- Prefect foresight not longer than forecast length
    if(mSettings(m, 't_jump') + mSettings(m, 't_perfectForesight')
       > max(mSettings(m, 't_forecastLengthUnchanging'),
             mSettings(m, 't_forecastLengthDecreasesFrom')),
        put log "!!! Error in model ", m.tl:0 /;
        put log "!!! Error: t_jump + t_perfectForesight > max(t_forecastLengthUnchanging, t_forecastLengthDecreasesFrom)"/;
        abort "Period of perfect foresight cannot be longer than forecast horizon";
    );

635
636
637
638
639
* --- Reserve structure checks ------------------------------------------------

    loop(restypeDirectionNode(restype, up_down, node),
        // Check that 'update_frequency' is longer than 't_jump'
        if(p_nReserves(node, restype, 'update_frequency') < mSettings(m, 't_jump'),
640
641
            put log '!!! Error occurred on p_nReserves ' node.tl:0 ',' restype.tl:0 /;
            put log '!!! Abort: The update_frequency parameter should be longer than or equal to t_jump!' /;
642
            abort "The 'update_frequency' parameter should be longer than or equal to 't_jump'!";
643
644
645
646
        ); // END if('update_frequency' < 't_jump')

        // Check that 'update_frequency' is divisible by 't_jump'
        if(mod(p_nReserves(node, restype, 'update_frequency'), mSettings(m, 't_jump')) <> 0,
647
648
            put log '!!! Error occurred on p_nReserves ' node.tl:0 ',' restype.tl:0 /;
            put log '!!! Abort: The update_frequency parameter should be divisible by t_jump!' /;
649
            abort "The 'update_frequency' parameter should be divisible by 't_jump'!";
650
651
652
        ); // END if(mod('update_frequency'))

        // Check if the first interval is long enough for proper commitment of reserves
653
        if(mInterval(m, 'lastStepInIntervalBlock', 'c000') < p_nReserves(node, restype, 'update_frequency') + p_nReserves(node, restype, 'gate_closure'),
654
            put log '!!! Error occurred on p_nReserves ' node.tl:0 ',' restype.tl:0 /;
Niina Helistö's avatar
Niina Helistö committed
655
656
            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!";
657
658
        ); // END if
    ); // END loop(restypeDirectionNode)
659
660
661
662
663

* --- 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' /;
664
665
        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!";
666
667
    ); // END if(smax)

668
669
* --- Check if time intervals are aggregated before 't_trajectoryHorizon' -----

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

675
676
677
678
679
680
681
682
683
684
685
686
* --- 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()

687
); // END loop(m)
688