3b_periodicLoop.gms 29.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
* --- Initialize unnecessary parameters and variables in order to save memory -
20
* =============================================================================
21

22
23
* --- Variables ---------------------------------------------------------------

24
25
26
27
28
// Free Variables
Option clear = v_gen;
Option clear = v_state;
Option clear = v_genRamp;
Option clear = v_transfer;
29
Option clear = v_ICramp;
30
31
32
33
// Integer Variables
Option clear = v_online_MIP;
Option clear = v_invest_MIP;
Option clear = v_investTransfer_MIP;
34
35
// Binary Variables
Option clear = v_help_inc;
36
37
38
// SOS2 Variables
Option clear = v_sos2;
// Positive Variables
39
40
Option clear = v_startup_LP;
Option clear = v_startup_MIP;
41
42
Option clear = v_shutdown_LP;
Option clear = v_shutdown_MIP;
43
44
45
46
47
48
49
50
51
52
Option clear = v_genRampUpDown;
Option clear = v_spill;
Option clear = v_transferRightward;
Option clear = v_transferLeftward;
Option clear = v_resTransferRightward;
Option clear = v_resTransferLeftward;
Option clear = v_reserve;
Option clear = v_investTransfer_LP;
Option clear = v_online_LP;
Option clear = v_invest_LP;
53
Option clear = v_gen_inc;
54
55
56
57
58
// Feasibility control
Option clear = v_stateSlack;
Option clear = vq_gen;
Option clear = vq_resDemand;
Option clear = vq_resMissing;
Niina Helistö's avatar
Niina Helistö committed
59
Option clear = vq_capacity;
60

61
62
* --- Equations ---------------------------------------------------------------

63
64
65
66
// Objective Function, Energy Balance, and Reserve demand
Option clear = q_obj;
Option clear = q_balance;
Option clear = q_resDemand;
67
68
69
70
Option clear = q_resDemandLargestInfeedUnit;
Option clear = q_rateOfChangeOfFrequencyUnit;
Option clear = q_rateOfChangeOfFrequencyTransfer;
Option clear = q_resDemandLargestInfeedTransfer;
71
72
73

// Unit Operation
Option clear = q_maxDownward;
Niina Helistö's avatar
Niina Helistö committed
74
Option clear = q_maxDownwardOfflineReserve;
75
Option clear = q_maxUpward;
Niina Helistö's avatar
Niina Helistö committed
76
Option clear = q_maxUpwardOfflineReserve;
77
Option clear = q_reserveProvision;
Niina Helistö's avatar
Niina Helistö committed
78
Option clear = q_reserveProvisionOnline;
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
Option clear = q_startshut;
Option clear = q_startuptype;
Option clear = q_onlineLimit;
Option clear = q_onlineMinUptime;
Option clear = q_onlineCyclic;
Option clear = q_onlineOnStartUp;
Option clear = q_offlineAfterShutdown;
Option clear = q_genRamp;
Option clear = q_rampUpLimit;
Option clear = q_rampDownLimit;
Option clear = q_rampUpDown;
Option clear = q_rampSlack;
Option clear = q_conversionDirectInputOutput;
Option clear = q_conversionSOS2InputIntermediate;
Option clear = q_conversionSOS2Constraint;
Option clear = q_conversionSOS2IntermediateOutput;
95
Option clear = q_conversionIncHR;
96
Option clear = q_conversionIncHRMaxOutput;
97
98
99
Option clear = q_conversionIncHRBounds;
Option clear = q_conversionIncHR_help1;
Option clear = q_conversionIncHR_help2;
100
Option clear = q_unitEqualityConstraint;
Niina Helistö's avatar
Niina Helistö committed
101
Option clear = q_unitGreaterThanConstraint;
102
*Option clear = q_commodityUseLimit;
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126

// Energy Transfer
Option clear = q_transfer;
Option clear = q_transferRightwardLimit;
Option clear = q_transferLeftwardLimit;
Option clear = q_resTransferLimitRightward;
Option clear = q_resTransferLimitLeftward;
Option clear = q_reserveProvisionRightward;
Option clear = q_reserveProvisionLeftward;

// State Variables
Option clear = q_stateSlack;
Option clear = q_stateUpwardLimit;
Option clear = q_stateDownwardLimit;
Option clear = q_boundStateMaxDiff;
Option clear = q_boundCyclic;

// Policy
Option clear = q_inertiaMin;
Option clear = q_instantaneousShareMax;
Option clear = q_constrainedOnlineMultiUnit;
Option clear = q_capacityMargin;
Option clear = q_constrainedCapMultiUnit;
Option clear = q_emissioncap;
127
Option clear = q_energyMax;
128
129
Option clear = q_energyShareMax;
Option clear = q_energyShareMin;
130
Option clear = q_ReserveShareMax;
131
132
133

* --- Temporary Time Series ---------------------------------------------------

134
135
136
137
138
139
140
141
142
// Initialize temporary time series
Option clear = ts_unit_;
*Option clear = ts_effUnit_;
*Option clear = ts_effGroupUnit_;
Option clear = ts_influx_;
Option clear = ts_cf_;
Option clear = ts_unit_;
Option clear = ts_reserveDemand_;
Option clear = ts_node_;
143
144
Option clear = ts_vomCost_;
Option clear = ts_startupCost_;
145

146

147
* =============================================================================
148
* --- Determine the forecast-intervals included in the current solve ----------
149
* =============================================================================
150

151
// Determine the time steps of the current solve
152
tSolveFirst = ord(tSolve);  // tSolveFirst: the start of the current solve, t0 used only for initial values
153

154
155
156
* --- Build the forecast-time structure using the intervals -------------------

// Initializing forecast-time structure sets
157
Option clear = p_stepLength;
158
Option clear = msft;
159
Option clear = mft;
160
Option clear = ft;
161
Option clear = sft;
162
Option clear = mst_start, clear = mst_end;
163
$ifthen declared scenario
164
if(mSettings(mSolve, 'scenarios'),  // Only clear these if using long-term scenarios
165
166
    Options clear = s_active, clear = s_scenario, clear = ss,
            clear = p_msProbability, clear = ms_central;
167
);
168
$endif
169

170

171
// Initialize the set of active t:s, counters and interval time steps
172
Option clear = t_active;
173
Option clear = dt_active;
174
Option clear = tt_block;
175
Option clear = cc;
176
tCounter = 1;
177
count_sample = 1;
178

179
// Determine the set of active interval counters (or blocks of intervals)
180
cc(counter)${ mInterval(mSolve, 'stepsPerInterval', counter) }
181
    = yes;
182

183
184
185
186
187
188
// Update tForecastNext
tForecastNext(mSolve)
    ${ tSolveFirst >= tForecastNext(mSolve) }
    = tForecastNext(mSolve) + mSettings(mSolve, 't_forecastJump');

// Calculate forecast length
189
currentForecastLength
190
    = max(  mSettings(mSolve, 't_forecastLengthUnchanging'),  // Unchanging forecast length would remain the same
191
            mSettings(mSolve, 't_forecastLengthDecreasesFrom') - [mSettings(mSolve, 't_forecastJump') - {tForecastNext(mSolve) - tSolveFirst}] // While decreasing forecast length has a fixed horizon point and thus gets shorter
192
            );   // Larger forecast horizon is selected
193

194
195
196
// Is there any case where t_forecastLength should be larger than t_horizon? Could happen if one doesn't want to join forecasts at the end of the solve horizon.
// If not, add a check for currentForecastLength <= mSettings(mSolve, 't_horizon')
// and change the line below to 'tSolveLast = ord(tSolve) + mSettings(mSolve, 't_horizon');'
197
tSolveLast = ord(tSolve) + mSettings(mSolve, 't_horizon');  // tSolveLast: the end of the current solve
198
Option clear = t_current;
199
200
201
202
t_current(t_full(t))
    ${  ord(t) >= tSolveFirst
        and ord (t) <= tSolveLast
        }
203
204
    = yes;

205
// Loop over the defined blocks of intervals
206
loop(cc(counter),
207
208
209
210
211
212
213

    // Abort if stepsPerInterval is less than one
    if(mInterval(mSolve, 'stepsPerInterval', counter) < 1,
        abort "stepsPerInterval < 1 is not defined!";
    );  // END IF stepsPerInterval

    // Time steps within the current block
214
215
    option clear = tt;
    tt(t_current(t))
216
        ${ord(t) >= tSolveFirst + tCounter
217
218
          and ord(t) <= min(tSolveFirst
                            + mInterval(mSolve, 'lastStepInIntervalBlock', counter),
219
220
                            tSolveLast)
         } = yes;
221

222
223
224
    // Store the interval time steps for each interval block (counter)
    tt_block(counter, tt) = yes;

Erkka Rinne's avatar
Erkka Rinne committed
225
226
227
    // Initialize tInterval
    Option clear = tt_interval;

Erkka Rinne's avatar
Erkka Rinne committed
228
229
    // If stepsPerInterval equals one, simply use all the steps within the block
    if(mInterval(mSolve, 'stepsPerInterval', counter) = 1,
230
231
        // Include all time steps within the block
        tt_interval(tt(t)) = yes;
Erkka Rinne's avatar
Erkka Rinne committed
232
233
234

    // If stepsPerInterval exceeds 1 (stepsPerInterval < 1 not defined)
    elseif mInterval(mSolve, 'stepsPerInterval', counter) > 1,
235
236
237
238
239
240

        // Calculate the displacement required to reach the corresponding active time step from any time step
        dt_active(tt(t)) = - (mod(ord(t) - tSolveFirst - tCounter, mInterval(mSolve, 'stepsPerInterval', counter)));

        // Select the active time steps within the block
        tt_interval(tt(t))${ not dt_active(t) } = yes;
241

242
    ); // END ELSEIF intervalLenght
243

244
245
246
247
248
249
250
    // Calculate the interval length in hours
    p_stepLength(mf(mSolve, f_solve), tt_interval(t))
      = mInterval(mSolve, 'stepsPerInterval', counter) * mSettings(mSolve, 'stepLengthInHours');
    p_stepLengthNoReset(mf(mSolve, f_solve), tt_interval(t)) = p_stepLength(mSolve, f_solve, t);

    // Determine the combinations of forecasts and intervals
    // Include the t_jump for the realization
251
    ft(f_solve, tt_interval(t))
252
253
254
       ${ord(t) <= tSolveFirst + max(mSettings(mSolve, 't_jump'),
                                     min(mSettings(mSolve, 't_perfectForesight'),
                                         currentForecastLength))
255
256
257
258
         and mf_realization(mSolve, f_solve)
        } = yes;

    // Include the full horizon for the central forecast
259
    ft(f_solve, tt_interval(t))
260
261
262
      ${ord(t) > tSolveFirst + max(mSettings(mSolve, 't_jump'),
                                   min(mSettings(mSolve, 't_perfectForesight'),
                                       currentForecastLength))
263
        and (mf_central(mSolve, f_solve)
264
             or mSettings(mSolve, 'forecasts') = 0)
265
266
267
       } = yes;

    // Include up to forecastLength for remaining forecasts
268
    ft(f_solve, tt_interval(t))
269
270
      ${not mf_central(mSolve, f_solve)
        and not mf_realization(mSolve, f_solve)
271
272
273
        and ord(t) > tSolveFirst + max(mSettings(mSolve, 't_jump'),
                                       min(mSettings(mSolve, 't_perfectForesight'),
                                           currentForecastLength))
274
275
276
277
278
279
        and ord(t) <= tSolveFirst + currentForecastLength
       } = yes;

    // Update tActive
    t_active(tt_interval) = yes;

280
281
282
283
284
285
286
287
288
289
290
    // Update tCounter for the next block of intervals
    tCounter = mInterval(mSolve, 'lastStepInIntervalBlock', counter) + 1;

); // END loop(counter)

// Reset initial sample start and end times if using scenarios
if(mSettings(mSolve, 'scenarios'),
    Option clear = msStart, clear = msEnd;
    msStart(ms_initial) = 1;
    msEnd(ms_initial) = currentForecastLength + 1;
);
291

292
$ifthen defined scenario
293
294
295
296
// Create stochastic programming scenarios
// Select root sample and central forecast
loop((ms_initial(mSolve, s_), mf_central(mSolve, f)),
    s_active(s_) = yes;
297
    p_msProbability(mSolve, s_)$mSettings(mSolve, 'scenarios') = 1;
298
    loop(scenario $p_scenProbability(scenario),
299
300
301
302
303
304
305
306
307
308
309
        s_scenario(s_, scenario) = yes;
        if(mSettings(mSolve, 'scenarios') > 1,
            loop(ft(f, t)$(ord(t) >= msEnd(mSolve, s_) + tSolveFirst),
                loop(s$(ord(s) = mSettings(mSolve, 'samples') + count_sample),
                    s_active(s) = yes;
                    ms_central(mSolve, s) = yes;
                    s_scenario(s, scenario) = yes;
                    p_msProbability(mSolve, s) = p_scenProbability(scenario);
                    msStart(mSolve, s) = ord(t) - tSolveFirst;
                    msEnd(mSolve, s) = ord(t) - tSolveFirst
                                              + p_stepLength(mSolve, f, t);
310
                );
311
312
313
314
315
316
317
318
319
320
321
                count_sample = count_sample + 1;
            );
        elseif mSettings(mSolve, 'scenarios') = 1,
            loop(ms(mSolve, s)$(not sameas(s, s_)),
                s_active(s) = yes;
                ms_central(mSolve, s) = yes;
                p_msProbability(mSolve, s) = 1;
                s_scenario(s, scenario) = yes;
                msStart(mSolve, s) = msEnd(mSolve, s_);
                msEnd(mSolve, s) = msStart(mSolve, s_)
                                   + mSettings(mSolve, 't_horizon');
322
323
324
            );
        );
    );
325
326
327
    ms(ms_central(mSolve, s)) = yes;
    msf(ms_central(mSolve, s), f) = yes;
);
328
329
$endif

330
331
// Loop over defined samples
loop(msf(mSolve, s, f)$msStart(mSolve, s),
332
                      // Move the samples along with the dispatch if scenarios are used
333
334
    sft(s, ft(f, t))${ord(t) > msStart(mSolve, s) + tSolveFirst - 1
                      and ord(t) < msEnd(mSolve, s) + tSolveFirst
335
336
337
338
339
340
                      and mSettings(mSolve, 'scenarios')
                     } = yes;
                      // Otherwise do not move the samples along with the rolling horizon
    sft(s, ft(f, t))${ord(t) > msStart(mSolve, s)
                      and ord(t) <= msEnd(mSolve, s)
                      and not mSettings(mSolve, 'scenarios')
341
342
343
                     } = yes;
);

344
345
346
347
348
// Update the model specific sets and the reversed dimension set
mft(mSolve, ft(f, t)) = yes;
ms(mSolve, s)$ms(mSolve, s) = s_active(s);
msf(mSolve, s, f)$msf(mSolve, s, f) = s_active(s);
msft(mSolve, sft(s, f, t)) = yes;
349

350
* Build stochastic tree by definfing previous samples
351
$ifthen defined scenario
352
Option clear = s_prev;
353
loop(scenario $p_scenProbability(scenario),
354
355
356
357
358
    loop(s_scenario(s, scenario),
        if(not ms_initial(mSolve, s), ss(s, s_prev) = yes);
        Option clear = s_prev; s_prev(s) = yes;
    );
);
359
360
$endif

361
362
363
364
365

* --- Define sample offsets for creating stochastic scenarios -----------------

Option clear = dt_scenarioOffset;

366
$ifthen defined scenario
367
368
369
370
371
372
373
374
375
376
377
loop(s_scenario(s, scenario)$(ord(s) > 1 and ord(scenario) > 1),
    loop(gn_scenarios(grid, node, timeseries),
         dt_scenarioOffset(grid, node, timeseries, s)
             = (ord(scenario) - 1) * mSettings(mSolve, 'scenarioLength');
    );

    loop(gn_scenarios(flow, node, timeseries),
        dt_scenarioOffset(flow, node, timeseries, s)
            = (ord(scenario) - 1) * mSettings(mSolve, 'scenarioLength');
    );
);
378
$endif
379
380


381
382
383
384
385
386
387
* --- Determine various other forecast-time sets required for the model -------

// Set of realized intervals in the current solve
Option clear = ft_realized;
ft_realized(ft(f_solve, t))
    ${  mf_realization(mSolve, f_solve)
        and ord(t) <= tSolveFirst + mSettings(mSolve, 't_jump')
Topi Rasku's avatar
Topi Rasku committed
388
389
        }
    = yes;
390

391
392
393
Option clear = sft_realized;
sft_realized(sft(s, ft_realized(f_solve, t))) = yes;

394
395
// Update the set of realized intervals in the whole simulation so far
ft_realizedNoReset(ft_realized(f, t)) = yes;
396
sft_realizedNoReset(sft_realized(s, f, t)) = yes;
397
398
399
400
401
402
403
404
405
406
407
// Update the set of realized intervals in the whole simulation so far, including model and sample dimensions
msft_realizedNoReset(msft(mSolve, s, ft_realized(f, t))) = yes;

// Include the necessary amount of historical timesteps to the active time step set of the current solve
loop(ft_realizedNoReset(f, t),
    t_active(t)
        ${  ord(t) <= tSolveFirst
            and ord(t) > tSolveFirst + tmp_dt // Strict inequality accounts for tSolvefirst being one step before the first ft step.
            }
        = yes;
); // END loop(ft_realizedNoReset
408

409
// Time step displacement to reach previous time step
410
option clear = dt;
Topi Rasku's avatar
Topi Rasku committed
411
option clear = dt_next;
412
tmp = max(tSolveFirst + tmp_dt, 1); // The ord(t) of the first time step in t_active, cannot decrease below 1 to avoid referencing time steps before t000000
413
414
loop(t_active(t),
    dt(t) = tmp - ord(t);
Topi Rasku's avatar
Topi Rasku committed
415
    dt_next(t+dt(t)) = -dt(t);
416
417
418
    tmp = ord(t);
); // END loop(t_active)

419
420
421
422
423
424
425
426
427
428
429
430
// First model ft
Option clear = mft_start;
mft_start(mf_realization(mSolve, f), tSolve)
    = yes
;
// Last model fts
Option clear = mft_lastSteps;
mft_lastSteps(mSolve, ft(f,t))
    ${ not dt_next(t) }
    = yes
;

431
432
433
434
435
436
437
438
// Sample start and end intervals
loop(ms(mSolve, s),
    tmp = 1;
    tmp_ = 1;
    loop(t_active(t),
        if(tmp and ord(t) - tSolveFirst + 1 > msStart(mSolve, s),
            mst_start(mSolve, s, t) = yes;
            tmp = 0;
Niina Helistö's avatar
Niina Helistö committed
439
        );
440
441
442
443
444
445
446
447
448
449
450
451
        if(tmp_ and ord(t) - tSolveFirst + 1 > msEnd(mSolve, s),
            mst_end(mSolve, s, t+dt(t)) = yes;
            tmp_ = 0;
        );
    ); // END loop(t_active)
    // If the last interval of a sample is in mft_lastSteps, the method above does not work
    if(tmp_,
        mst_end(mSolve, s, t)${sum(f_solve, mft_lastSteps(mSolve, f_solve, t))} = yes;
    );
); // END loop(ms)
// Displacement from the first interval of a sample to the previous interval is always -1,
// except for stochastic samples
452
453
454
dt(t_active(t))
    ${ sum(ms(mSolve, s)$(not ms_central(mSolve, s)), mst_start(mSolve, s, t)) }
    = -1;
Niina Helistö's avatar
Niina Helistö committed
455

Niina Helistö's avatar
Niina Helistö committed
456
// Forecast index displacement between realized and forecasted intervals
457
// NOTE! This set cannot be reset without references to previously solved time steps in the stochastic tree becoming ill-defined!
458
459
460
df(f_solve(f), t_active(t))${ ord(t) <= tSolveFirst + max(mSettings(mSolve, 't_jump'),
                                                          min(mSettings(mSolve, 't_perfectForesight'),
                                                              currentForecastLength))}
461
    = sum(mf_realization(mSolve, f_), ord(f_) - ord(f));
462
// Displacement to reach the realized forecast
463
Option clear = df_realization;
464
loop(mf_realization(mSolve, f_),
465
    df_realization(ft(f, t))$[ord(t) <= tSolveFirst + currentForecastLength]
466
467
      = ord(f_) - ord(f);
);
468
// Central forecast for the long-term scenarios comes from a special forecast label
469
Option clear = df_scenario;
470
if(mSettings(mSolve, 'scenarios') >= 1,
471
472
    loop((msft(ms_central(mSolve, s), f, t), mf_scenario(mSolve, f_)),
        df_scenario(ft(f, t)) = ord(f_) - ord(f);
473
474
    );
);
475
476
// Check that df_forecast and df_scenario do not overlap
loop(ft(f, t),
477
  if(df_realization(f, t) <> 0 and df_scenario(f, t) <> 0,
478
479
      put log "!!! Overlapping period of using realization and scenarios"/;
      put log "!!! Check forecast lengths, `gn_scenarios` and `gn_forecasts`"/;
480
      execError = execError + 1;
481
482
  );
);
483

Niina Helistö's avatar
Niina Helistö committed
484
// Forecast displacement between central and forecasted intervals at the end of forecast horizon
485
Option clear = df_central; // This can be reset.
486
487
df_central(ft(f,t))${   ord(t) > tSolveFirst + currentForecastLength - p_stepLength(mSolve, f, t) / mSettings(mSolve, 'stepLengthInHours')
                        and ord(t) <= tSolveFirst + currentForecastLength
488
489
                        and not mf_realization(mSolve, f)
                        }
490
    = sum(mf_central(mSolve, f_), ord(f_) - ord(f));
491

Niina Helistö's avatar
Niina Helistö committed
492
// Forecast index displacement between realized and forecasted intervals, required for locking reserves ahead of (dispatch) time.
493
Option clear = df_reserves;
494
495
496
497
498
499
500
501
502
df_reserves(grid, node, restype, ft(f, t))
    ${  p_gnReserves(grid, node, restype, 'update_frequency')
        and p_gnReserves(grid, node, restype, 'gate_closure')
        and ord(t) <= tSolveFirst + p_gnReserves(grid, node, restype, 'gate_closure')
                                  + p_gnReserves(grid, node, restype, 'update_frequency')
                                  - mod(tSolveFirst - 1 + p_gnReserves(grid, node, restype, 'gate_closure')
                                                    + p_gnReserves(grid, node, restype, 'update_frequency')
                                                    - p_gnReserves(grid, node, restype, 'update_offset'),
                                    p_gnReserves(grid, node, restype, 'update_frequency'))
503
504
        }
    = sum(f_${ mf_realization(mSolve, f_) }, ord(f_) - ord(f)) + Eps; // The Eps ensures that checks to see if df_reserves exists return positive even if the displacement is zero.
505
506
507
508
509
510
511
512
513
514
Option clear = df_reservesGroup;
df_reservesGroup(group, restype, ft(f, t))
    ${  p_groupReserves(group, restype, 'update_frequency')
        and p_groupReserves(group, restype, 'gate_closure')
        and ord(t) <= tSolveFirst + p_groupReserves(group, restype, 'gate_closure')
                                  + p_groupReserves(group, restype, 'update_frequency')
                                  - mod(tSolveFirst - 1 + p_groupReserves(group, restype, 'gate_closure')
                                                    + p_groupReserves(group, restype, 'update_frequency')
                                                    - p_groupReserves(group, restype, 'update_offset'),
                                    p_groupReserves(group, restype, 'update_frequency'))
515
516
517
        }
    = sum(f_${ mf_realization(mSolve, f_) }, ord(f_) - ord(f)) + Eps; // The Eps ensures that checks to see if df_reserves exists return positive even if the displacement is zero.

518
// Set of ft-steps where the reserves are locked due to previous commitment
519
Option clear = ft_reservesFixed;
520
ft_reservesFixed(group, restype, f_solve(f), t_active(t))
521
522
    ${  mf_realization(mSolve, f)
        and not tSolveFirst = mSettings(mSolve, 't_start') // No reserves are locked on the first solve!
523
524
525
526
        and p_groupReserves(group, restype, 'update_frequency')
        and p_groupReserves(group, restype, 'gate_closure')
        and ord(t) <= tSolveFirst + p_groupReserves(group, restype, 'gate_closure')
                                  + p_groupReserves(group, restype, 'update_frequency')
Erkka Rinne's avatar
Erkka Rinne committed
527
                                  - mod(tSolveFirst - 1
528
                                          + p_groupReserves(group, restype, 'gate_closure')
Erkka Rinne's avatar
Erkka Rinne committed
529
                                          - mSettings(mSolve, 't_jump')
530
531
532
                                          + p_groupReserves(group, restype, 'update_frequency')
                                          - p_groupReserves(group, restype, 'update_offset'),
                                        p_groupReserves(group, restype, 'update_frequency'))
Erkka Rinne's avatar
Erkka Rinne committed
533
                                  - mSettings(mSolve, 't_jump')
534
        and not [   restypeReleasedForRealization(restype) // Free desired reserves for the to-be-realized time steps
535
536
                    and ft_realized(f, t)
                    ]
537
538
        }
    = yes;
539

540
541
542
// Form a temporary clone of t_current
option clear = tt;
tt(t_current) = yes;
543
544
// Group each full time step under each active time step for time series aggregation.
option clear = tt_aggregate;
545
tt_aggregate(t_current(t+dt_active(t)), tt(t))
546
547
    = yes;

548
* =============================================================================
549
* --- Defining unit aggregations and ramps ------------------------------------
550
* =============================================================================
551

552
// Units active on each ft
553
Option clear = uft;
554
555
556
uft(unit, ft(f, t))${   (   [
                                ord(t) <= tSolveFirst + p_unit(unit, 'lastStepNotAggregated')
                                and (unit_aggregated(unit) or unit_noAggregate(unit)) // Aggregated and non-aggregate units
557
                            ]
558
559
560
561
                            or
                            [
                                ord(t) > tSolveFirst + p_unit(unit, 'lastStepNotAggregated')
                                and (unit_aggregator(unit) or unit_noAggregate(unit)) // Aggregator and non-aggregate units
562
                            ]
563
564
565
                        )
                        and not sameas(unit, 'empty')
                     }
566
// only units with capacities or investment option
567
    = yes;
568

569
// Units are not active before or after their lifetime
570
571
uft(unit, ft(f, t))${   [ ord(t) < p_unit(unit, 'becomeAvailable') and p_unit(unit, 'becomeAvailable') ]
                        or [ ord(t) >= p_unit(unit, 'becomeUnavailable') and p_unit(unit, 'becomeUnavailable') ]
572
573
574
                        }
    = no;

575
576
577
578
579
580
581
582
583
584
585
586
// First ft:s for each aggregator unit
Option clear = uft_aggregator_first;
loop(unit${unit_aggregator(unit)},
    tmp = card(t);
    loop(uft(unit, f, t),
        if(ord(t) < tmp,
            tmp = ord(t)
        );
    );
    uft_aggregator_first(uft(unit, f, t))${ord(t) = tmp} = yes;
);

587
// Active (grid, node, unit) on each ft
588
Option clear = gnuft;
589
gnuft(gn(grid, node), uft(unit, f, t))${    gnu(grid, node, unit)  }
590
591
    = yes
;
592
// Active (grid, node, unit, slack, up_down) on each ft step with ramp restrictions
593
594
595
596
597
Option clear = gnuft_rampCost;
gnuft_rampCost(gnu(grid, node, unit), slack, ft(f, t))${ gnuft(grid, node, unit, f, t)
                                                         and p_gnuBoundaryProperties(grid, node, unit, slack, 'rampCost')
                                                         }
    = yes;
598
// Active (grid, node, unit) on each ft step with ramp restrictions
599
Option clear = gnuft_ramp;
600
601
gnuft_ramp(gnuft(grid, node, unit, f, t))${ p_gnu(grid, node, unit, 'maxRampUp')
                                            OR p_gnu(grid, node, unit, 'maxRampDown')
602
                                            OR sum(slack, gnuft_rampCost(grid, node, unit, slack, f, t))
603
604
605
606
                                            }
    = yes;

* --- Defining unit efficiency groups etc. ------------------------------------
607

608
// Initializing
609
610
Option clear = suft;
Option clear = sufts;
611

612
613
614
// Loop over the defined efficiency groups for units
loop(effLevelGroupUnit(effLevel, effGroup, unit)${ mSettingsEff(mSolve, effLevel) },
    // Determine the used effGroup for each uft
615
616
    suft(effGroup, uft(unit, f, t))${   ord(t) >= tSolveFirst + mSettingsEff(mSolve, effLevel - 1) + 1
                                        and ord(t) <= tSolveFirst + mSettingsEff(mSolve, effLevel) }
617
        = yes;
618
); // END loop(effLevelGroupUnit)
619
620
621

// Determine the efficiency selectors for suft
sufts(suft(effGroup, unit, f, t), effSelector)${    effGroupSelector(effGroup, effSelector) }
622
623
    = yes
;
624

625
// Units with online variables on each ft
626
Option clear = uft_online;
627
628
Option clear = uft_onlineLP;
Option clear = uft_onlineMIP;
629
630
Option clear = uft_onlineLP_withPrevious;
Option clear = uft_onlineMIP_withPrevious;
631

632
// Determine the intervals when units need to have online variables.
633
634
loop(effOnline(effSelector),
    uft_online(uft(unit, f, t))${ suft(effOnline, unit, f, t) }
635
        = yes;
636
637
638
639
); // END loop(effOnline)
uft_onlineLP(uft(unit, f, t))${ suft('directOnLP', unit, f, t) }
    = yes;
uft_onlineMIP(uft_online(unit, f, t)) = uft_online(unit, f, t) - uft_onlineLP(unit, f, t);
640

641
642
643
644
// Units with start-up and shutdown trajectories
Option clear = uft_startupTrajectory;
Option clear = uft_shutdownTrajectory;

645
// Determine the intervals when units need to follow start-up and shutdown trajectories.
646
647
648
loop(runUpCounter(unit, 'c000'), // Loop over units with meaningful run-ups
    uft_startupTrajectory(uft_online(unit, f, t))
        ${ ord(t) <= tSolveFirst + mSettings(mSolve, 't_trajectoryHorizon') }
649
        = yes;
650
651
652
653
); // END loop(runUpCounter)
loop(shutdownCounter(unit, 'c000'), // Loop over units with meaningful shutdowns
    uft_shutdownTrajectory(uft_online(unit, f, t))
        ${ ord(t) <= tSolveFirst + mSettings(mSolve, 't_trajectoryHorizon') }
654
        = yes;
655
); // END loop(shutdownCounter)
656

657
* -----------------------------------------------------------------------------
658
* --- Displacements for start-up and shutdown decisions -----------------------
659
660
* -----------------------------------------------------------------------------

661
662
* --- Start-up decisions ------------------------------------------------------

Niina Helistö's avatar
Niina Helistö committed
663
664
// Calculate dt_toStartup: in case the unit becomes online in the current time interval,
// displacement needed to reach the time interval where the unit was started up
665
Option clear = dt_toStartup;
666
loop(runUpCounter(unit, 'c000'), // Loop over units with meaningful run-ups
667
    dt_toStartup(unit, t_active(t))$(ord(t) <= tSolveFirst + mSettings(mSolve, 't_trajectoryHorizon'))
668
        = - p_u_runUpTimeIntervalsCeil(unit) + dt_active(t - p_u_runUpTimeIntervalsCeil(unit));
669
); // END loop(runUpCounter)
670

671
* --- Shutdown decisions ------------------------------------------------------
672
673

// Calculate dt_toShutdown: in case the generation of the unit becomes zero in
Niina Helistö's avatar
Niina Helistö committed
674
// the current time interval, displacement needed to reach the time interval where
675
676
// the shutdown decisions was made
Option clear = dt_toShutdown;
677
loop(shutdownCounter(unit, 'c000'), // Loop over units with meaningful shutdowns
678
    dt_toShutdown(unit, t_active(t))$(ord(t) <= tSolveFirst + mSettings(mSolve, 't_trajectoryHorizon'))
679
        = - p_u_shutdownTimeIntervalsCeil(unit) + dt_active(t - p_u_shutdownTimeIntervalsCeil(unit))
680
); // END loop(runUpCounter)
681

682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
* --- Historical Unit LP and MIP information ----------------------------------

uft_onlineLP_withPrevious(uft_onlineLP(unit, f, t)) = yes;
uft_onlineMIP_withPrevious(uft_onlineMIP(unit, f, t)) = yes;

// Units with online variables on each active ft starting at t0
loop(mft_start(mSolve, f, t_), // Check the uft_online used on the first time step of the current solve
    uft_onlineLP_withPrevious(unit, f, t_active(t)) // Include all historical t_active
        ${  uft_onlineLP(unit, f, t_+1) // Displace by one to reach the first current time step
            and ord(t) <= tSolveFirst // Include all historical t_active
            }
         = yes;
    uft_onlineMIP_withPrevious(unit, f, t_active(t)) // Include all historical t_active
        ${  uft_onlineMIP(unit, f, t_+1) // Displace by one to reach the first current time step
            and ord(t) <= tSolveFirst // Include all historical t_active
            }
        = yes;
); // END loop(mft_start)

701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
// Historical Unit LP and MIP information for models with multiple samples
// If this is the very first solve
if(tSolveFirst = mSettings(mSolve, 't_start'),
    // Sample start intervals
    loop(mst_start(mSolve, s, t),
        uft_onlineLP_withPrevious(unit, f, t+dt(t)) // Displace by one to reach the time step just before the sample
            ${  uft_onlineLP(unit, f, t)
                }
             = yes;
        uft_onlineMIP_withPrevious(unit, f, t+dt(t)) // Displace by one to reach the time step just before the sample
            ${  uft_onlineMIP(unit, f, t)
                }
            = yes;
    ); // END loop(mst_start)
); // END if(tSolveFirst)