1e_inputs.gms 19 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
20
21
* =============================================================================
* --- Load Input Data ---------------------------------------------------------
* =============================================================================

22
$gdxin  'input/inputData.gdx'
Juha Kiviluoma's avatar
Juha Kiviluoma committed
23
$loaddc grid
24
$loaddc node
Juha Kiviluoma's avatar
Juha Kiviluoma committed
25
$loaddc flow
26
$loaddc unittype
27
$loaddc unit
28
$loaddc unitUnittype
29
$loaddc fuel
30
$loaddc unitUnit_aggregate
Juha Kiviluoma's avatar
Juha Kiviluoma committed
31
32
$loaddc uFuel
$loaddc effLevelGroupUnit
33
34
35
$loaddc p_gn
$loaddc p_gnn
$loaddc p_gnu
36
$loaddc p_gnuBoundaryProperties
Juha Kiviluoma's avatar
Juha Kiviluoma committed
37
$loaddc p_unit
38
$loaddc ts_unit
39
40
41
$loaddc restype
$loaddc restypeDirection
$loaddc restypeReleasedForRealization
42
$loaddc p_nReserves
43
$loaddc p_nuReserves
44
$loaddc ts_reserveDemand
45
$loaddc p_gnBoundaryPropertiesForStates
46
$loaddc p_gnPolicy
Juha Kiviluoma's avatar
Juha Kiviluoma committed
47
$loaddc p_uFuel
48
$loaddc flowUnit
Juha Kiviluoma's avatar
Juha Kiviluoma committed
49
50
$loaddc gngnu_fixedOutputRatio
$loaddc gngnu_constrainedOutputRatio
51
$loaddc emission
52
$loaddc p_fuelEmission
53
$loaddc ts_cf
54
$loaddc ts_fuelPriceChange
55
$loaddc ts_influx
56
$loaddc ts_nodeState
57
$loaddc t_invest
58
$loaddc group
59
60
61
62
$loaddc uGroup
$loaddc gnuGroup
$loaddc gn2nGroup
$loaddc gnGroup
63
$loaddc p_groupPolicy
64
$loaddc p_groupPolicy3D
65
66
67
68
69
70
71
72
73
74
75
$gdxin

$ontext
 * Load stochastic scenarios
 $batinclude 'inc/gdxload_fluctuation.inc' wind
 $batinclude 'inc/gdxload_fluctuation.inc' solar
 $ifthen exist 'input/scenarios_hydro.gdx'
    $$gdxin 'input/scenarios_hydro.gdx'
 $endif
 $gdxin
$offtext
76
77
78
$ifthen exist 'input/changes.inc'
   $$include 'input/changes.inc'
$endif
79

80

81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
* =============================================================================
* --- Initialize Unit Related Sets & Parameters Based on Input Data -----------
* =============================================================================

* --- Unit Aggregation --------------------------------------------------------

// Define unit aggregation sets
unit_aggregate(unit)${ sum(unit_, unitUnit_aggregate(unit, unit_)) }
    = yes; // Set of aggregate units
unit_noAggregate(unit)${ unit(unit) - unit_aggregate(unit) - sum(unit_, unitUnit_aggregate(unit_, unit))}
    = yes; // Set of units that are not aggregated into any aggregate, or are not aggregates themselves

// Process data for unit aggregations
// Aggregate maxGen as the sum of aggregated maxGen
p_gnu(grid, node, unit_aggregate(unit), 'maxGen')
    = sum(unit_${unitUnit_aggregate(unit, unit_)},
        + p_gnu(grid, node, unit_, 'maxGen')
        );
// Aggregate maxCons as the sum of aggregated maxCons
p_gnu(grid, node, unit_aggregate(unit), 'maxCons')
    = sum(unit_${unitUnit_aggregate(unit, unit_)},
        + p_gnu(grid, node, unit_, 'maxCons')
        );

* --- Generate Unit Related Sets ----------------------------------------------

// Set of all existing gnu
gnu(grid, node, unit)${ p_gnu(grid, node, unit, 'maxGen')
                        or p_gnu(grid, node, unit, 'maxCons')
                        or p_gnu(grid, node, unit, 'unitSizeGen')
                        or p_gnu(grid, node, unit, 'unitSizeCons')
                        or p_gnu(grid, node, unit, 'maxGenCap')
                        or p_gnu(grid, node, unit, 'maxConsCap')
                        }
    = yes;
// Reduce the grid dimension
nu(node, unit) = sum(grid, gnu(grid, node, unit));

// Separation of gnu into inputs and outputs
gnu_output(gnu(grid, node, unit))${ p_gnu(grid, node, unit, 'maxGen')
                                    or p_gnu(grid, node, unit, 'maxGenCap')
                                    or p_gnu(grid, node, unit, 'unitSizeGen')
                                    }
    = yes;
gnu_input(gnu(grid, node, unit))${  p_gnu(grid, node, unit, 'maxCons')
                                    or p_gnu(grid, node, unit, 'maxConsCap')
                                    or p_gnu(grid, node, unit, 'unitSizeCons')
                                    }
    = yes;

// Units connecting gn-gn pairs
gn2gnu(grid, node_input, grid_, node_output, unit)${    gnu_input(grid, node_input, unit)
                                                        and gnu_output(grid_, node_output, unit)
                                                        }
    = yes;

// Units with reserve provision capabilities
nuRescapable(restypeDirection(restype, up_down), nu(node, unit))${ p_nuReserves(node, unit, restype, up_down) }
    = yes;

// Units with minimum load requirements
unit_minload(unit)${    p_unit(unit, 'op00') > 0 // If the first defined operating point is between 0 and 1, then the unit is considered to have a min load limit
                        and p_unit(unit, 'op00') < 1
                        }
    = yes;
146

147
148
149
150
151
152
// Units with flows/fuels
unit_flow(unit)${ sum(flow, flowUnit(flow, unit)) }
    = yes;
unit_fuel(unit)${ sum(fuel, uFuel(unit, 'main', fuel)) }
    = yes;

153
154
// Units with special startup properties
// All units can cold start (default start category)
155
156
157
158
159
unitStarttype(unit, starttype('cold'))${ p_unit(unit, 'startCostCold')
                                         or p_unit(unit, 'startFuelConsCold')
                                         or p_unit(unit, 'rampSpeedToMinLoad')
                                       }
    = yes;
160
// Units with parameters regarding hot/warm starts
161
unitStarttype(unit, starttypeConstrained)${ p_unit(unit, 'startWarmAfterXhours')
162
163
164
165
                                            or p_unit(unit, 'startCostHot')
                                            or p_unit(unit, 'startFuelConsHot')
                                            or p_unit(unit, 'startCostWarm')
                                            or p_unit(unit, 'startFuelConsWarm')
166
                                            or p_unit(unit, 'startColdAfterXhours')
167
168
169
                                            }
    = yes;

170
// Units with investment variables
171
unit_investLP(unit)${  not p_unit(unit, 'investMIP')
172
                       and p_unit(unit, 'maxUnitCount')
173
174
175
176
177
178
179
180
181
182
183
184
185
186
                        }
    = yes;
unit_investMIP(unit)${  p_unit(unit, 'investMIP')
                        and p_unit(unit, 'maxUnitCount')
                        }
    = yes;

* --- Unit Related Parameters -------------------------------------------------

// Assume values for critical unit related parameters, if not provided by input data
// If the unit does not have efficiency set, it is 1
p_unit(unit, 'eff00')${ not p_unit(unit, 'eff00') }
    = 1;

187
// In case number of units has not been defined it is 1 except for units with investments allowed.
188
p_unit(unit, 'unitCount')${ not p_unit(unit, 'unitCount')
189
190
                            and not unit_investMIP(unit)
                            and not unit_investLP(unit)
191
192
193
194
195
196
                            }
    = 1;

// By default add outputs in order to get the total capacity of the unit
p_unit(unit, 'outputCapacityTotal')${ not p_unit(unit, 'outputCapacityTotal') }
    = sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'maxGen'));
197
198
p_unit(unit, 'unitOutputCapacityTotal')
    = sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'unitSizeGen'));
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215

// Assume unit sizes based on given maximum capacity parameters and unit counts if able
p_gnu(grid, node, unit, 'unitSizeGen')${    p_gnu(grid, node, unit, 'maxGen')
                                            and p_unit(unit, 'unitCount')
                                            }
    = p_gnu(grid, node, unit, 'maxGen') / p_unit(unit, 'unitCount');  // If maxGen and unitCount are given, calculate unitSizeGen based on them.
p_gnu(grid, node, unit, 'unitSizeCons')${   p_gnu(grid, node, unit, 'maxCons')
                                            and p_unit(unit, 'unitCount')
                                            }
    = p_gnu(grid, node, unit, 'maxCons') / p_unit(unit, 'unitCount');  // If maxCons and unitCount are given, calculate unitSizeCons based on them.
p_gnu(grid, node, unit, 'unitSizeTot')
    = p_gnu(grid, node, unit, 'unitSizeGen') + p_gnu(grid, node, unit, 'unitSizeCons');
p_gnu(grid, node, unit, 'unitSizeGenNet')
    = p_gnu(grid, node, unit, 'unitSizeGen') - p_gnu(grid, node, unit, 'unitSizeCons');

// Determine unit startup parameters based on data
// Hot startup parameters
216
p_uNonoperational(unitStarttype(unit, 'hot'), 'min')
217
    = p_unit(unit, 'minShutdownHours');
218
p_uNonoperational(unitStarttype(unit, 'hot'), 'max')
219
    = p_unit(unit, 'startWarmAfterXhours');
220
p_uStartup(unitStarttype(unit, 'hot'), 'cost', 'unit')
221
222
    = p_unit(unit, 'startCostHot')
        * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'unitSizeGen'));
223
p_uStartup(unitStarttype(unit, 'hot'), 'consumption', 'unit')
224
225
226
227
    = p_unit(unit, 'startFuelConsHot')
        * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'unitSizeGen'));

// Warm startup parameters
228
p_uNonoperational(unitStarttype(unit, 'warm'), 'min')
229
    = p_unit(unit, 'startWarmAfterXhours');
230
p_uNonoperational(unitStarttype(unit, 'warm'), 'max')
231
    = p_unit(unit, 'startColdAfterXhours');
232
p_uStartup(unitStarttype(unit, 'warm'), 'cost', 'unit')
233
234
    = p_unit(unit, 'startCostWarm')
        * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'unitSizeGen'));
235
p_uStartup(unitStarttype(unit, 'warm'), 'consumption', 'unit')
236
237
238
239
    = p_unit(unit, 'startFuelConsWarm')
        * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'unitSizeGen'));

// Cold startup parameters
240
241
p_uNonoperational(unitStarttype(unit, 'cold'), 'min')
    = p_unit(unit, 'startColdAfterXhours');
242
p_uStartup(unit, 'cold', 'cost', 'unit')
243
    = p_unit(unit, 'startCostCold')
244
245
        * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'unitSizeGen'));
p_uStartup(unit, 'cold', 'consumption', 'unit')
246
    = p_unit(unit, 'startFuelConsCold')
247
248
249
        * sum(gnu_output(grid, node, unit), p_gnu(grid, node, unit, 'unitSizeGen'));

// Determine unit emission costs
250
p_unitFuelEmissionCost(unit_fuel, fuel, emission)${ sum(param_fuel, uFuel(unit_fuel, param_fuel, fuel)) }
251
252
    = p_fuelEmission(fuel, emission)
        / 1e3 // NOTE!!! Conversion to t/MWh from t/GWh in data
253
254
255
256
257
258
259
260
261
262
263
        * sum(gnu_output(grid, node, unit_fuel),
            + p_gnPolicy(grid, node, 'emissionTax', emission)  // Weighted average of emission costs from different output energy types
                * [ + p_gnu(grid, node, unit_fuel, 'maxGen')
                    + p_gnu(grid, node, unit_fuel, 'unitSizeGen')${not p_gnu(grid, node, unit_fuel, 'maxGen')}
                    ] // END * p_gnPolicy
        ) // END sum(gnu_output)
        / sum(gnu_output(grid, node, unit_fuel), // Weighted average of emission costs from different output energy types
            + p_gnu(grid, node, unit_fuel, 'maxGen')
            + p_gnu(grid, node, unit_fuel, 'unitSizeGen')$(not p_gnu(grid, node, unit_fuel, 'maxGen'))
        ) // END sum(gnu_output)
;
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324

* =============================================================================
* --- Generate Node Related Sets Based on Input Data --------------------------
* =============================================================================

* --- Node Connectivity -------------------------------------------------------

// Node pairs connected via transfer links
gn2n(grid, from_node, to_node)${    p_gnn(grid, from_node, to_node, 'transferCap')
                                    or p_gnn(grid, from_node, to_node, 'transferLoss')
                                    or p_gnn(grid, from_node, to_node, 'transferCapBidirectional')
                                    or p_gnn(grid, to_node, from_node, 'transferCapBidirectional')
                                    or p_gnn(grid, from_node, to_node, 'transferCapInvLimit')
                                    }
    = yes;

// Node pairs with relatively bound states
gnn_boundState(grid, node, node_)${ p_gnn(grid, node, node_, 'boundStateMaxDiff') }
    = yes;

// Node pairs connected via energy diffusion
gnn_state(grid, node, node_)${  p_gnn(grid, node, node_, 'diffCoeff')
                                or gnn_boundState(grid, node, node_)
                                }
    = yes;

// Generate the set for transfer links where the order of the first node must be smaller than the order of the second node
Option clear = gn2n_directional;
gn2n_directional(gn2n(grid, node, node_))${ ord(node) < ord(node_) }
    = yes;
gn2n_directional(gn2n(grid, node, node_))${ ord(node) > ord(node_)
                                            and not gn2n(grid, node_, node)
                                            }
    = yes;

* --- Node States -------------------------------------------------------------

// States with slack variables
gn_stateSlack(grid, node)${ sum((slack, useConstantOrTimeSeries), p_gnBoundaryPropertiesForStates(grid, node, slack, useConstantOrTimeSeries)) }
    = yes;

// Nodes with states
gn_state(grid, node)${  gn_stateSlack(grid, node)
                        or p_gn(grid, node, 'energyStoredPerUnitOfState')
                        or sum((stateLimits, useConstantOrTimeSeries), p_gnBoundaryPropertiesForStates(grid, node, stateLimits, useConstantOrTimeSeries))
                        or sum(useConstantOrTimeSeries, p_gnBoundaryPropertiesForStates(grid, node, 'reference', useConstantOrTimeSeries))
                        }
    = yes;

// Existing grid-node pairs
gn(grid, node)${    sum(unit, gnu(grid, node, unit)
                    or gn_state(grid, node))
                    }
    = yes;

// Nodes with spill permitted
node_spill(node)${ sum((grid, spillLimits, useConstantOrTimeSeries), p_gnBoundaryPropertiesForStates(grid, node, spillLimits, useConstantOrTimeSeries)) }
    = yes;

// Assume values for critical node related parameters, if not provided by input data
// Boundary multiplier
325
326
327
p_gnBoundaryPropertiesForStates(gn(grid, node), param_gnBoundaryTypes, 'multiplier')${  not p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'multiplier')
                                                                                        and sum(param_gnBoundaryProperties, p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, param_gnBoundaryProperties))
    } = 1; // If multiplier has not been set, set it to 1 by default
328

329
330
331
* --- Other Node Properties ---------------------------------------------------

// Nodes with flows
332
333
334
flowNode(flow, node)${  sum((f, t), ts_cf(flow, node, f, t))
                        and sum(grid, gn(grid, node))
                        }
335
336
    = yes;

337
338
339
* =============================================================================
* --- Reserves Sets & Parameters ----------------------------------------------
* =============================================================================
340
341

// Nodes with reserve requirements
342
343
344
345
restypeDirectionNode(restypeDirection(restype, up_down), node)${    p_nReserves(node, restype, up_down)
                                                                    or p_nReserves(node, restype, 'use_time_series')
                                                                    }
    = yes;
346

347
348
// Assume values for critical reserve related parameters, if not provided by input data
// Reserve contribution "reliability" assumed to be perfect if not provided in data
349
350
p_nuReserves(nu(node, unit), restype, 'reserveContribution')${  not p_nuReserves(node, unit, restype, 'reserveContribution')
                                                                and sum(up_down, nuRescapable(restype, up_down, node, unit))
351
352
353
354
355
356
                                                                }
    = 1;

* =============================================================================
* --- Data Integrity Checks ---------------------------------------------------
* =============================================================================
357

358
* Check the integrity of node connection related data
359
Option clear = count;
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
loop(gn2n(grid, node, node_),
    count = count + 1; // Count the gn2n indeces to make finding the errors easier.
    // Check if the bidirectional transfer parameter exists for this link.
    if(p_gnn(grid, node, node_, 'transferCapBidirectional'),
        // Check for conflicting bidirectional transfer capacities.
        if(p_gnn(grid, node, node_, 'transferCapBidirectional') <> p_gnn(grid, node_, node, 'transferCapBidirectional'),
            put log '!!! Error occurred on gn2n link #' count;
            abort "Conflicting 'transferCapBidirectional' parameters!"
        );
        // Check for conflicting one-directional and bidirectional transfer capacities.
        if(p_gnn(grid, node, node_, 'transferCapBidirectional') < p_gnn(grid, node, node_, 'transferCap') OR (p_gnn(grid, node, node_, 'transferCapBidirectional') < p_gnn(grid, node_, node, 'transferCap')),
            put log '!!! Error occurred on gn2n link #' count;
            abort "Parameter 'transferCapBidirectional' must be greater than or equal to defined one-directional transfer capacities!"
        );
    );
);

377
* Check the integrity of efficiency approximation related data
378
Option clear = tmp; // Log the unit index for finding the error easier.
379
loop( unit,
380
    tmp = ord(unit); // Increase the unit counter
381
    // Check that 'op' is defined correctly
382
    Option clear = count; // Initialize the previous op to zero
383
    loop( op,
384
        abort${p_unit(unit, op) + 1${not p_unit(unit, op)} < count} "param_unit 'op's must be defined as zero or positive and increasing!", tmp, count;
385
        count = p_unit(unit, op);
386
    );
387
388
    // Check that efficiency approximations have sufficient data
    loop( effLevelGroupUnit(effLevel, effSelector, unit),
389
390
391
392
393
394
395
        loop( op__${p_unit(unit, op__) = smax(op, p_unit(unit, op))}, // Loop over the 'op's to find the last defined data point.
            // Lambda  - Has been commented out, since it is ok to improve the efficiency curve by using extra lambda points.
            //loop( lambda${sameas(lambda, effSelector)}, // Loop over the lambdas to find the 'magnitude' of the approximation
                //display count_lambda, count_lambda2, unit.tl;
            //    if(ord(lambda) > ord(op__), put log '!!! Error occurred on unit ' unit.tl:25 ' with effLevel ' effLevel.tl:10 ' with effSelector ' lambda.tl:8); // Display unit that causes error
            //    abort${ord(lambda) > ord(op__)} "Order of the lambda approximation cannot exceed the number of efficiency data points!"
            // );
396
            // DirectOn
397
398
399
400
401
            loop( op_${p_unit(unit, op_) = smin(op${p_unit(unit, op)}, p_unit(unit, op))}, // Loop over the 'op's to find the first nonzero 'op' data point.
                if(effDirectOn(effSelector) AND ord(op__) = ord(op_) AND not p_unit(unit, 'section') AND not p_unit(unit, 'opFirstCross'),
                    put log '!!! Error occurred on unit #' tmp; // Display unit that causes error, NEEDS WORK
                    abort "directOn requires two efficiency data points with nonzero 'op' or 'section' or 'opFirstCross'!";
                );
402
403
404
            );
        );
    );
405
406
);