periodicInit.gms 13.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
$ontext
This file is part of Backbone.

Backbone is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

Backbone is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with Backbone.  If not, see <http://www.gnu.org/licenses/>.
$offtext

18
19
*  Generate model rules from basic patterns defined in the model definition files 

20
21
22
23
24
* Set the time for the next available forecast unless it has been preset.
loop(m,
    tForecastNext(m) = mSettings(m, 't_forecastStart');
);

25
26
* Check the modelSolves for preset patterns for model solve timings
* If not found, then use mSettings to set the model solve timings
27
28
loop(m,
    if(sum(t$modelSolves(m, t), 1) = 0,
29
        t_skip_counter = 0;
30
31
        loop(t$( ord(t) = mSettings(m, 't_start') + mSettings(m, 't_jump') * t_skip_counter and ord(t) <= mSettings(m, 't_end') ),
            modelSolves(m, t)=yes;
32
33
            t_skip_counter = t_skip_counter + 1;
        );
34
        p_stepLengthNoReset(m, f, t) = no;
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
    );
);

* Select samples for the model
loop(m,
    // Set samples in use for the models
    if (not sum(s, ms(m, s)),  // unless they have been provided as input
        ms(m, s)$(ord(s) <= mSettings(m, 'samples')) = yes;
        if (mSettings(m, 'samples') = 0,     // Use all samples if mSettings/samples is 0
            ms(m, s) = yes;
        );
    );
    // Set forecasts in use for the models
    if (not sum(f, mf(m, f)),  // unless they have been provided as input
        mf(m, f)$(ord(f) <= 1 + mSettings(m, 'forecasts')) = yes;  // realization needs one f, therefore 1 + number of forecasts
    );
51
    msf(m, s, f)$(ms(m, s) and mf(m, f)) = yes;
52
53
54
);


55
* Calculate the length of the time series
56
continueLoop = 1;
57
tmp = 0;
58
loop(mSolve${continueLoop},
59
    loop(gn(grid, node),
60
61
        tmp = max(sum(t${ts_influx(grid, node, 'f00', t)}, 1), tmp); // Find the maximum length of the given influx time series
        tmp = max(sum(t${ts_nodeState(grid, node, 'reference', 'f00', t)}, 1), tmp); // Find the maximum length of the given node state time series
62
    );
63
    ts_length = tmp;
64
65
    ct(t)$(
            ord(t) > ts_length
66
67
        and ord(t) <= ts_length + max(mSettings(mSolve, 't_forecastLength'), mSettings(mSolve, 't_horizon'))
        and ord(t) <= mSettings(mSolve, 't_end') + max(mSettings(mSolve, 't_forecastLength'), mSettings(mSolve, 't_horizon'))
68
69
    ) = -ts_length;
    continueLoop = 0;
Juha Kiviluoma's avatar
Juha Kiviluoma committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
);

loop(m,
    tmp = 0;
    continueLoop = 1;
    loop(counter$continueLoop,
        if(not mInterval(m, 'intervalEnd', counter),
            continueLoop = 0;
        else
            abort$(mod(mInterval(m, 'intervalEnd', counter) - mInterval(m, 'intervalEnd', counter-1), mInterval(m, 'intervalLength', counter))) "IntervalLength is not evenly divisible within the interval", m, continueLoop;
            continueLoop = continueLoop + 1;
        );
    );
);


86
* Parse through effLevelGroupUnit and convert selected effSelectors into sets representing those selections
Juha Kiviluoma's avatar
Juha Kiviluoma committed
87
88
loop(unit,
    loop(effLevel$sum(m, mSettingsEff(m, effLevel)),
89
90
91
92
93
94
        // effSelector using DirectOff
        if (sum(effDirectOff, effLevelGroupUnit(effLevel, effDirectOff, unit)),
            loop(effDirectOff${effLevelGroupUnit(effLevel, effDirectOff, unit)},
                    effGroup(effDirectOff) = yes;
                    effGroupSelector(effDirectOff, effDirectOff) = yes;
                    effGroupSelectorUnit(effDirectOff, unit, effDirectOff) = yes;
Juha Kiviluoma's avatar
Juha Kiviluoma committed
95
96
            );
        );
97
98
99
100
101
102
103

        // effSelector using DirectOn
        if (sum(effDirectOn, effLevelGroupUnit(effLevel, effDirectOn, unit)),
            loop(effDirectOn${effLevelGroupUnit(effLevel, effDirectOn, unit)},
                    effGroup(effDirectOn) = yes;
                    effGroupSelector(effDirectOn, effDirectOn) = yes;
                    effGroupSelectorUnit(effDirectOn, unit, effDirectOn) = yes;
Juha Kiviluoma's avatar
Juha Kiviluoma committed
104
105
106
107
            );
        );

        // effSelectors using Lambda
108
109
110
111
112
113
        if (sum(effLambda, effLevelGroupUnit(effLevel, effLambda, unit)),
            loop(effLambda${effLevelGroupUnit(effLevel, effLambda, unit)},
                    effGroup(effLambda) = yes;
                    loop(effLambda_${ord(effLambda_) <= ord(effLambda)},
                            effGroupSelector(effLambda, effLambda_) = yes;
                            effGroupSelectorUnit(effLambda, unit, effLambda_) = yes;
Juha Kiviluoma's avatar
Juha Kiviluoma committed
114
115
116
117
118
                    );
            );
        );

        // Calculate parameters for units using direct input output conversion without online variable
119
        loop(effGroupSelectorUnit(effDirectOff, unit, effDirectOff_)$effLevelGroupUnit(effLevel, effDirectOff, unit),
120
            p_effUnit(effDirectOff, unit, effDirectOff_, 'lb') = 0; // No min load for the DirectOff approximation
121
            p_effUnit(effDirectOff, unit, effDirectOff_, 'op') = 1; // No max load for the DirectOff approximation
122
            p_effUnit(effDirectOff, unit, effDirectOff_, 'slope') = 1 / smax(eff${p_unit(unit, eff)}, p_unit(unit, eff)); // Uses maximum found (nonzero) efficiency.
123
124
125
126
            p_effUnit(effDirectOff, unit, effDirectOff_, 'section') = 0; // No section for the DirectOff approximation
        );

        // Calculate parameters for units using direct input output conversion with online variable
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
        loop(effGroupSelectorUnit(effDirectOn, unit, effDirectOn_)$effLevelGroupUnit(effLevel, effDirectOn, unit),
            p_effUnit(effDirectOn, unit, effDirectOn_, 'lb') = p_unit(unit, 'op00'); // op00 contains the possible min load of the unit
            p_effUnit(effDirectOn, unit, effDirectOn_, 'op') = smax(op, p_unit(unit, op)); // Maximum load determined by the largest 'op' parameter found in data
            tmp_count_op = 0;
            loop(op$p_unit(unit, op),  // Determine the last operating point in use
                tmp_count_op = ord(op);
            );
            loop(op__$(ord(op__) = tmp_count_op), // Find the maximum defined 'op'.
                loop(eff__${ord(eff__) = ord(op__)},                     // ...  and the corresponding 'eff'.
                     if(p_unit(unit, 'op00') = 0,  // If the minimum operating point is at zero, then the section and slope are calculated with the assumption that the efficiency curve crosses at opFirstCross
                         heat_rate =  // heat rate at the cross between real efficiency curve and approximated efficiency curve
                            1 / {
                                  + p_unit(unit, 'eff00') * [ p_unit(unit, op__) - p_unit(unit, 'opFirstCross') ] / [ p_unit(unit, op__) - p_unit(unit, 'op00') ]
                                  + p_unit(unit, eff__) * [ p_unit(unit, 'opFirstCross') - p_unit(unit, 'op00') ] / [ p_unit(unit, op__) - p_unit(unit, 'op00') ]
                                };
                         p_effGroupUnit(effDirectOn, unit, 'section') = p_unit(unit, 'section');
                         p_effGroupUnit(effDirectOn, unit, 'section')$(not p_effGroupUnit(effDirectOn, unit, 'section')) =  // Unless section has been defined, it is calculated based on the opFirstCross
                              p_unit(unit, 'opFirstCross')
                                * ( heat_rate - 1 / p_unit(unit, eff__) )
                                / ( p_unit(unit, op__) - p_unit(unit, 'opFirstCross') );
                         p_effUnit(effDirectOn, unit, effDirectOn_, 'slope') =
                             1 / p_unit(unit, eff__) - p_effGroupUnit(effDirectOn, unit, 'section') / p_unit(unit, op__);
                     else  // If the minimum operating point is above zero, then the approximate efficiency curve crosses the real efficiency curve at minimum and maximum.
                          // Calculating the slope based on the first nonzero and the last defined data points.
                         p_effUnit(effDirectOn, unit, effDirectOn_, 'slope') =
                           + (p_unit(unit, op__) / p_unit(unit, eff__) - p_unit(unit, 'op00') / p_unit(unit, 'eff00'))
                                / (p_unit(unit, op__) - p_unit(unit, 'op00'));
                         // Calculating the section based on the slope and the last defined point.
155
                         p_effGroupUnit(effDirectOn, unit, 'section') =
156
157
158
                             ( 1 / p_unit(unit, eff__) - p_effUnit(effDirectOn, unit, effDirectOn_, 'slope') )
                                * p_unit(unit, op__);
                     );
159
160
                );
            );
Juha Kiviluoma's avatar
Juha Kiviluoma committed
161
162
163
        );

        // Calculate lambdas
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
        loop(effGroupSelectorUnit(effLambda, unit, effLambda_)$effLevelGroupUnit(effLevel, effLambda, unit),
            p_effUnit(effLambda, unit, effLambda_, 'lb') = p_unit(unit, 'op00'); // op00 contains the min load of the unit
            tmp_count_op = 0;
            loop(op$p_unit(unit, op),  // Determine the last operating point in use
                tmp_count_op = ord(op);
            );
            tmp_op = (ord(effLambda_)-1) / (ord(effLambda) - 1); //  Calculate the relative location of the operating point in the lambdas
            p_effUnit(effLambda, unit, effLambda_, 'op') = tmp_op;  // Copy the operating point to the p_effUnit
            tmp_op$(ord(effLambda_) = 1 and not p_unit(unit, 'op00') and not p_unit(unit, 'section')) = p_unit(unit, 'opFirstCross');  // Copy the cross between the p_unit efficiency curve and opFirstCross to the temporary op parameter for further calculations
            // tmp_op falls between two p_unit defined operating points or then it is equal to one of them
            loop((op_, op__)$({[tmp_op > p_unit(unit, op_) and tmp_op < p_unit(unit, op__) and ord(op_) = ord(op__) - 1] or [p_unit(unit, op_) = tmp_op and ord(op_) = ord(op__)]} and ord(op__) <= tmp_count_op),
                loop((eff_, eff__)$(ord(op_) = ord(eff_) and ord(op__) = ord(eff__)),
                    tmp_dist = p_unit(unit, op__) - p_unit(unit, op_); // Calculate the distance between the operating points (zero if the points are the same)
                    if (tmp_dist, // If the operating points are not the same
                        heat_rate =  // Heat rate is a weighted average of the heat rates at the p_unit operating points
                            1 / {
                                  + p_unit(unit, eff_) * [ p_unit(unit, op__) - tmp_op ] / tmp_dist
                                  + p_unit(unit, eff__) * [ tmp_op - p_unit(unit, op_) ] / tmp_dist
                                };
                    else  // If the operating point is the same, the the heat rate can be used directly
                        heat_rate = 1 / p_unit(unit, eff_);
185
                    );
186
187
188
189
190
                    if (ord(effLambda_) = 1,
                        if(p_unit(unit, 'op00') or p_unit(unit, 'section'),  // If the min. load of the unit is not zero or the section has been pre-defined, then section is copied directly from the unit properties
                            p_effGroupUnit(effLambda, unit, 'section') = p_unit(unit, 'section');
                        else
                            p_effGroupUnit(effLambda, unit, 'section') =  // Calculate section based on the opFirstCross, which has been calculated into p_effUnit(effLambda, unit, effLambda_, 'op')
191
                                tmp_op
192
193
                                  * ( heat_rate - 1 / p_unit(unit, eff__) )
                                  / ( p_unit(unit, 'op01') - tmp_op );
194
195
                        );
                    );
196
197
                    p_effUnit(effLambda, unit, effLambda_, 'slope')$(ord(effLambda_) > 1 or p_unit(unit, 'op00')) =
                        heat_rate - p_effGroupUnit(effLambda, unit, 'section') / tmp_op;
198
199
                );
            );
Juha Kiviluoma's avatar
Juha Kiviluoma committed
200
        );
201
202
    ); // END LOOP OVER effLevel
); // END LOOP OVER unit
203
204
205
206
207
208


// Calculate unit wide parameters for each efficiency group
loop(unit,
    loop(effLevel$sum(m, mSettingsEff(m, effLevel)),
        loop(effLevelGroupUnit(effLevel, effGroup, unit),
209
            p_effGroupUnit(effGroup, unit, 'op') = smax(effSelector$effGroupSelectorUnit(effGroup, unit, effSelector), p_effUnit(effGroup, unit, effSelector, 'op'));
210
            p_effGroupUnit(effGroup, unit, 'lb') = smin(effSelector${effGroupSelectorUnit(effGroup, unit, effSelector)}, p_effUnit(effGroup, unit, effSelector, 'lb'));
211
            p_effGroupUnit(effGroup, unit, 'slope') = smin(effSelector${effGroupSelectorUnit(effGroup, unit, effSelector)}, p_effUnit(effGroup, unit, effSelector, 'slope')); // NOTE! Uses maximum efficiency for the group.
212
        );
Juha Kiviluoma's avatar
Juha Kiviluoma committed
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
    );
);


* Ensure that efficiency levels extend to the end of the model horizon
loop(m,
    continueLoop = 0;
    loop(effLevel$mSettingsEff(m, effLevel),
        continueLoop = continueLoop + 1;
    );
    loop(effLevel$(ord(effLevel) = continueLoop),
        if (mSettingsEff(m, effLevel) <> mSettings(m, 't_horizon'),
            mSettingsEff(m, effLevel + 1) = mSettings(m, 't_horizon');
        );
    );
);

230
231
232
233

* Set slack direction
p_slackDirection(upwardSlack) = 1;
p_slackDirection(downwardSlack) = -1;