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
127
128

// 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;
Option clear = q_energyShareMax;
Option clear = q_energyShareMin;
129
Option clear = q_ReserveShareMax;
130
131
132

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

133
134
135
136
137
138
139
140
141
// 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_;
142
143
Option clear = ts_vomCost_;
Option clear = ts_startupCost_;
144

145

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

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

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

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

169

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

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

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

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

193
194
195
// 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');'
196
tSolveLast = ord(tSolve) + mSettings(mSolve, 't_horizon');  // tSolveLast: the end of the current solve
197
Option clear = t_current;
198
199
200
201
t_current(t_full(t))
    ${  ord(t) >= tSolveFirst
        and ord (t) <= tSolveLast
        }
202
203
    = yes;

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

    // 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
213
214
    option clear = tt;
    tt(t_current(t))
215
        ${ord(t) >= tSolveFirst + tCounter
216
217
          and ord(t) <= min(tSolveFirst
                            + mInterval(mSolve, 'lastStepInIntervalBlock', counter),
218
219
                            tSolveLast)
         } = yes;
220

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

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

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

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

        // 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;
240

241
    ); // END ELSEIF intervalLenght
242

243
244
245
246
247
248
249
    // 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
250
    ft(f_solve, tt_interval(t))
251
252
253
       ${ord(t) <= tSolveFirst + max(mSettings(mSolve, 't_jump'),
                                     min(mSettings(mSolve, 't_perfectForesight'),
                                         currentForecastLength))
254
255
256
257
         and mf_realization(mSolve, f_solve)
        } = yes;

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

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

    // Update tActive
    t_active(tt_interval) = yes;

279
280
281
282
283
284
285
286
287
288
289
    // 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;
);
290

291
$ifthen defined scenario
292
293
294
295
// Create stochastic programming scenarios
// Select root sample and central forecast
loop((ms_initial(mSolve, s_), mf_central(mSolve, f)),
    s_active(s_) = yes;
296
    p_msProbability(mSolve, s_)$mSettings(mSolve, 'scenarios') = 1;
297
    loop(scenario $p_scenProbability(scenario),
298
299
300
301
302
303
304
305
306
307
308
        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);
309
                );
310
311
312
313
314
315
316
317
318
319
320
                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');
321
322
323
            );
        );
    );
324
325
326
    ms(ms_central(mSolve, s)) = yes;
    msf(ms_central(mSolve, s), f) = yes;
);
327
328
$endif

329
330
// Loop over defined samples
loop(msf(mSolve, s, f)$msStart(mSolve, s),
331
                      // Move the samples along with the dispatch if scenarios are used
332
333
    sft(s, ft(f, t))${ord(t) > msStart(mSolve, s) + tSolveFirst - 1
                      and ord(t) < msEnd(mSolve, s) + tSolveFirst
334
335
336
337
338
339
                      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')
340
341
342
                     } = yes;
);

343
344
345
346
347
// 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;
348

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

360
361
362
363
364

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

Option clear = dt_scenarioOffset;

365
$ifthen defined scenario
366
367
368
369
370
371
372
373
374
375
376
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');
    );
);
377
$endif
378
379


380
381
382
383
384
385
386
* --- 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
387
388
        }
    = yes;
389

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

393
394
// Update the set of realized intervals in the whole simulation so far
ft_realizedNoReset(ft_realized(f, t)) = yes;
395
sft_realizedNoReset(sft_realized(s, f, t)) = yes;
396
397
398
399
400
401
402
403
404
405
406
// 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
407

408
// Time step displacement to reach previous time step
409
option clear = dt;
Topi Rasku's avatar
Topi Rasku committed
410
option clear = dt_next;
411
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
412
413
loop(t_active(t),
    dt(t) = tmp - ord(t);
Topi Rasku's avatar
Topi Rasku committed
414
    dt_next(t+dt(t)) = -dt(t);
415
416
417
    tmp = ord(t);
); // END loop(t_active)

418
419
420
421
422
423
424
425
426
427
428
429
// 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
;

430
431
432
433
434
435
436
437
// 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
438
        );
439
440
441
442
443
444
445
446
447
448
449
450
        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
451
452
453
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
454

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

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

Niina Helistö's avatar
Niina Helistö committed
491
// Forecast index displacement between realized and forecasted intervals, required for locking reserves ahead of (dispatch) time.
492
Option clear = df_reserves;
493
494
495
496
497
498
499
500
501
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'))
502
503
        }
    = 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.
504
505
506
507
508
509
510
511
512
513
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'))
514
515
516
        }
    = 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.

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

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

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

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

568
// Units are not active before or after their lifetime
569
570
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') ]
571
572
573
                        }
    = no;

574
575
576
577
578
579
580
581
582
583
584
585
// 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;
);

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

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

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

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

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

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

631
// Determine the intervals when units need to have online variables.
632
633
loop(effOnline(effSelector),
    uft_online(uft(unit, f, t))${ suft(effOnline, unit, f, t) }
634
        = yes;
635
636
637
638
); // 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);
639

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

644
// Determine the intervals when units need to follow start-up and shutdown trajectories.
645
646
647
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') }
648
        = yes;
649
650
651
652
); // 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') }
653
        = yes;
654
); // END loop(shutdownCounter)
655

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

660
661
* --- Start-up decisions ------------------------------------------------------

Niina Helistö's avatar
Niina Helistö committed
662
663
// 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
664
Option clear = dt_toStartup;
665
loop(runUpCounter(unit, 'c000'), // Loop over units with meaningful run-ups
666
    dt_toStartup(unit, t_active(t))$(ord(t) <= tSolveFirst + mSettings(mSolve, 't_trajectoryHorizon'))
667
        = - p_u_runUpTimeIntervalsCeil(unit) + dt_active(t - p_u_runUpTimeIntervalsCeil(unit));
668
); // END loop(runUpCounter)
669

670
* --- Shutdown decisions ------------------------------------------------------
671
672

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

681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
* --- 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)

700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
// 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)