4a_outputVariant.gms 12.4 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
* --- Recording realized parameter values -------------------------------------
20
* =============================================================================
21

22
* --- Result arrays required by model dynamics --------------------------------
Erkka Rinne's avatar
Erkka Rinne committed
23
24
25
26
27
28
29
30
31
32
if(tSolveFirst >= mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod') - mSettings(mSolve, 't_jump')
   and firstResultsOutputSolve,
    loop(msf(mSolve, s, f_solve),
        firstResultsOutputSolve = 0;
        r_state(gn_state(grid, node), f_solve, t) $[ord(t) = mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
          = v_state.l(grid, node, s, f_solve, t);
        r_online(unit, f_solve, t)$[unit_online(unit) and ord(t) = mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
          = v_online_LP.l(unit, s, f_solve, t)$unit_online_LP(unit)
              + v_online_MIP.l(unit, s, f_solve, t)$unit_online_MIP(unit);
    );
33
);
34
// Realized state history
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
loop(sft_realized(s, f, t),
    r_state(gn_state(grid, node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = v_state.l(grid, node, s, f, t);

    // Realized state history - initial state values in samples
    r_state(gn_state(grid, node), f_solve(f), t+dt(t))${   ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')
                                                           and sum(ms(mSolve, s), mst_start(mSolve, s, t))
                                                       }
        = v_state.l(grid, node, s, f, t+dt(t))
    ;
    // Realized unit online history
    r_online(uft_online(unit, f, t))$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = v_online_LP.l(unit, s, f, t)${ uft_onlineLP(unit, f, t)    }
            + v_online_MIP.l(unit, s, f, t)${  uft_onlineMIP(unit, f, t)   }
    ;
    // Unit startup and shutdown history
    r_startup(unit, starttype, f, t)${ uft_online(unit, f, t) and [ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')] }
52
53
        = v_startup_LP.l(unit, starttype, s, f, t)${ uft_onlineLP(unit, f, t) }
            + v_startup_MIP.l(unit, starttype, s, f, t)${ uft_onlineMIP(unit, f, t) }
54
55
    ;
    r_shutdown(uft_online(unit, f, t))$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
56
57
        = v_shutdown_LP.l(unit, s, f, t)${ uft_onlineLP(unit, f, t) }
            + v_shutdown_MIP.l(unit, s, f, t)${ uft_onlineMIP(unit, f, t) }
58
    ;
59
);
60

61
62
63
* --- Reserve results ---------------------------------------------------------

// Loop over reserve horizon, as the reserve variables use a different ft-structure due to commitment
64
loop((restypeDirectionNode(restype, up_down, node), sft(s, f, t))
65
    ${  ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')
66
67
68
69
70
        and ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')
        },

    // Reserve provisions of units
    r_reserve(nuRescapable(restype, up_down, node, unit), f+df_reserves(node, restype, f, t), t)
71
        ${  not [   restypeReleasedForRealization(restype)
72
                    and sft_realized(s, f+df_reserves(node, restype, f, t), t)
73
74
                    ]
            }
75
        = + v_reserve.l(restype, up_down, node, unit, s, f+df_reserves(node, restype, f, t), t)
76
          + sum(restype_$p_nuRes2Res(node, unit, restype_, up_down, restype),
77
              + v_reserve.l(restype_, up_down, node, unit, s, f+df_reserves(node, restype_, f, t), t)
78
                  * p_nuRes2Res(node, unit, restype_, up_down, restype)
79
            );
80

81
    // Reserve requirement due to N-1 reserve constraint
82
    r_resDemandLargestInfeedUnit(grid, f+df_reserves(node, restype, f, t), t, restype, up_down, node)
83
        ${ sum(unit_fail, nuRescapable(restype, up_down, node, unit_fail)) } // Calculate only for nodes with units that can fail.
84
85
        = smax(unit_fail(unit_), v_gen.l(grid, node, unit_, s, f, t) * p_nuReserves(node, unit_, restype, 'portion_of_infeed_to_reserve'));

86
    // Reserve transfer capacity for links defined out from this node
87
    r_resTransferRightward(restype, up_down, node, to_node, f+df_reserves(node, restype, f, t), t)
88
89
        ${  sum(grid, gn2n_directional(grid, node, to_node))
            and restypeDirectionNodeNode(restype, up_down, node, to_node)
90
            }
91
        = v_resTransferRightward.l(restype, up_down, node, to_node, s, f+df_reserves(node, restype, f, t), t);
92

93
94
95
    r_resTransferLeftward(restype, up_down, node, to_node, f+df_reserves(to_node, restype, f, t), t)
        ${  sum(grid, gn2n_directional(grid, node, to_node))
            and restypeDirectionNodeNode(restype, up_down, to_node, node)
96
            }
97
        = v_resTransferLeftward.l(restype, up_down, node, to_node, s, f+df_reserves(to_node, restype, f, t), t);
98
99

    // Dummy reserve demand changes
100
101
    r_qResDemand(restype, up_down, group, f+df_reservesGroup(group, restype, f, t), t)
        = vq_resDemand.l(restype, up_down, group, s, f+df_reservesGroup(group, restype, f, t), t);
102

103
104
    r_qResMissing(restype, up_down, group, f+df_reservesGroup(group, restype, f, t), t)
        = vq_resMissing.l(restype, up_down, group, s, f+df_reservesGroup(group, restype, f, t), t);
105

106
); // END loop(restypeDirectionNode, sft)
107

108
109
* --- Interesting results -----------------------------------------------------

110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
loop(sft_realized(s, f, t),
    // Unit generation and consumption
    r_gen(gnuft(grid, node, unit, f, t))$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = v_gen.l(grid, node, unit, s, f, t)
    ;
    // Fuel use of units
    r_fuelUse(fuel, uft(unit_fuel, f, t))$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = v_fuelUse.l(fuel, unit_fuel, s, f, t)
    ;
    // Transfer of energy between nodes
    r_transfer(gn2n(grid, from_node, to_node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = v_transfer.l(grid, from_node, to_node, s, f, t)
    ;
    // Energy spilled from nodes
    r_spill(gn(grid, node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = v_spill.l(grid, node, s, f, t)
    ;
127
128
);

129
// Total Objective function
130
131
r_totalObj(tSolve)
    = r_totalObj(tSolve - mSettings(mSolve, 't_jump')) + v_obj.l
132
;
133

134
// q_balance marginal values
135
loop(sft_realized(s, f, t),
136
137
    r_balanceMarginal(gn(grid, node), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = q_balance.m(grid, node, mSolve, s, f, t)
138
139
    ;
    // q_resDemand marginal values
140
141
    r_resDemandMarginal(restypeDirectionGroup(restype, up_down, group), f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = q_resDemand.m(restype, up_down, group, s, f, t)
142
143
144
145
146
    ;
    // v_stateSlack values for calculation of realized costs later on
    r_stateSlack(gn_stateSlack(grid, node), slack, f, t)$[ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')]
        = v_stateSlack.l(grid, node, slack, s, f, t)
    ;
147
);
148
// Unit investments
Niina Helistö's avatar
Niina Helistö committed
149
r_invest(unit)${unit_investLP(unit) or unit_investMIP(unit)}
150
151
    = sum(t_invest, v_invest_LP.l(unit, t_invest) + v_invest_MIP.l(unit, t_invest))
;
Niina Helistö's avatar
Niina Helistö committed
152
153
154
155
156
157
158
159
// Link investments
r_investTransfer(grid, node, node_, t_invest(t))${ p_gnn(grid, node, node_, 'transferCapInvLimit')
*                                                   and t_current(t)
                                                   and ord(t) <= tSolveFirst + mSettings(mSolve, 't_jump')
                                                   }
    = v_investTransfer_LP.l(grid, node, node_, t)
        + v_investTransfer_MIP.l(grid, node, node_, t) * p_gnn(grid, node, node_, 'unitSize')
;
160
161

* --- Feasibility results -----------------------------------------------------
162
loop(sft_realized(s, f, t),
163
// Dummy generation & consumption
164
r_qGen(inc_dec, gn(grid, node), f, t)
165
166
    ${  ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')
        }
167
    = vq_gen.l(inc_dec, grid, node, s, f, t)
168
;
169
170
171
172
173
174
// Dummy capacity
r_qCapacity(gn(grid, node), f, t)
    ${  ord(t) > mSettings(mSolve, 't_start') + mSettings(mSolve, 't_initializationPeriod')
        }
    = vq_capacity.l(grid, node, s, f, t)
;
175
);
176

177
* =============================================================================
178
* --- Diagnostics Results -----------------------------------------------------
179
* =============================================================================
180

181
// Only include these if '--diag=yes' given as a command line argument
182
$iftheni.diag %diag% == 'yes'
183
// Capacity factors for examining forecast errors
184
d_capacityFactor(flowNode(flow, node), sft(s, f_solve(f), t_current(t)))
185
186
187
188
    ${  msf(mSolve, s, f)
        and t_active(t)
        and sum(flowUnit(flow, unit), nu(node, unit))
        }
189
    = ts_cf_(flow, node, f, t, s)
190
        + ts_cf(flow, node, f, t + dt_scenarioOffset(flow, node, 'ts_cf', s))${ not ts_cf_(flow, node, f, t, s) }
191
        + Eps
Topi Rasku's avatar
Topi Rasku committed
192
193
;
// Temperature forecast for examining the error
194
d_nodeState(gn_state(grid, node), param_gnBoundaryTypes, sft(s, f_solve(f), t_current(t)))
195
196
197
198
    ${  p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')
        and t_active(t)
        and msf(mSolve, s, f)
        }
199
200
    = ts_node_(grid, node, param_gnBoundaryTypes, f, t, s)
        + ts_node(grid, node, param_gnBoundaryTypes, f, t)${ not ts_node_(grid, node, param_gnBoundaryTypes, f, t, s)}
201
202
203
        + Eps
;
// Influx forecast for examining the errors
204
d_influx(gn(grid, node), sft(s, f_solve(f), t_current(t)))
205
206
207
208
209
210
    ${  msf(mSolve, s, f)
        and t_active(t)
        }
    = ts_influx_(grid, node, f, t, s)
        + ts_influx(grid, node, f, t)${ not ts_influx_(grid, node, f, t, s)}
        + Eps
Topi Rasku's avatar
Topi Rasku committed
211
;
212
213
// Scenario values for time series
Options clear = d_state, clear = d_ts_scenarios; // Only keep latest results
214
215
216
217
loop(s_scenario(s, scenario),
    loop(mft_start(mSolve, f, t)$ms_initial(mSolve, s),
        d_state(gn_state(grid, node), scenario, f, t) = v_state.l(grid, node, s, f, t);
    );
218
219
220
    d_state(gn_state, scenario, ft)$sft(s, ft) = v_state.l(gn_state, s, ft) + eps;
    d_ts_scenarios('ts_influx', gn, scenario, ft)$sft(s, ft) = ts_influx_(gn, ft, s) + eps;
    d_ts_scenarios('ts_cf', flowNode, scenario, ft)$sft(s, ft) = ts_cf_(flowNode, ft, s) + eps;
221
);
222
$endif.diag
223
224

* --- Model Solve & Status ----------------------------------------------------
225

226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
// Model/solve status
if (mSolve('schedule'),
    r_solveStatus(tSolve,'modelStat')=schedule.modelStat;
    r_solveStatus(tSolve,'solveStat')=schedule.solveStat;
    r_solveStatus(tSolve,'totalTime')=schedule.etSolve;
    r_solveStatus(tSolve,'iterations')=schedule.iterUsd;
    r_solveStatus(tSolve,'nodes')=schedule.nodUsd;
    r_solveStatus(tSolve,'numEqu')=schedule.numEqu;
    r_solveStatus(tSolve,'numDVar')=schedule.numDVar;
    r_solveStatus(tSolve,'numVar')=schedule.numVar;
    r_solveStatus(tSolve,'numNZ')=schedule.numNZ;
    r_solveStatus(tSolve,'sumInfes')=schedule.sumInfes;
    r_solveStatus(tSolve,'objEst')=schedule.objEst;
    r_solveStatus(tSolve,'objVal')=schedule.objVal;
);
241
242
243
244
245
246
247
248
249
250
251
252
253
254
if (mSolve('invest'),
    r_solveStatus(tSolve,'modelStat')=invest.modelStat;
    r_solveStatus(tSolve,'solveStat')=invest.solveStat;
    r_solveStatus(tSolve,'totalTime')=invest.etSolve;
    r_solveStatus(tSolve,'iterations')=invest.iterUsd;
    r_solveStatus(tSolve,'nodes')=invest.nodUsd;
    r_solveStatus(tSolve,'numEqu')=invest.numEqu;
    r_solveStatus(tSolve,'numDVar')=invest.numDVar;
    r_solveStatus(tSolve,'numVar')=invest.numVar;
    r_solveStatus(tSolve,'numNZ')=invest.numNZ;
    r_solveStatus(tSolve,'sumInfes')=invest.sumInfes;
    r_solveStatus(tSolve,'objEst')=invest.objEst;
    r_solveStatus(tSolve,'objVal')=invest.objVal;
);
255

256
257