Commit 0ed10d22 authored by Juha Kiviluoma's avatar Juha Kiviluoma
Browse files

A working version of forecasts and dispatch loop within the scheduling loop

parent a5178141
......@@ -17,3 +17,9 @@ diffile.gdx
defModels/schedule.gms
defModels/storage.gms
EfficiencyPiecewise.xlsx
*.pptx
*.docx
*.pdf
*.gdx
*.xlsx
stosschver
......@@ -55,7 +55,7 @@ files log /''/, gdx, f_info /'output\info.txt'/;
options
optca = 0
optcr = 0.00001
profile = 8
* profile = 8
solvelink = %Solvelink.Loadlibrary%
bratio = 1
solveopt = merge
......@@ -82,17 +82,12 @@ $include 'input\3a_modelsInit.gms' // Sets that are defined over the whole mode
* === Simulation ==============================================================
loop(modelSolves(mSolve, tSolve),
$$include 'input\3b_modelsLoop.gms' // Set sets that define model scope
$$include 'inc\3c_setVariableLimits.gms' // Set new variable limits (.lo and .up)
$$ifi '%debug%' == 'yes' execute_unload 'output\debug.gdx'; // Output debugging information
$$include 'inc\3d_solve.gms' // Solve model(s)
put log;
put schedule.resGen;
put tSolveFirst;
putclose log;
$$include 'inc\3b_inputsLoop.gms' // Read input data that is updated within the loop
$$include 'input\3c_modelsLoop.gms' // Set sets that define model scope
$$include 'inc\3d_setVariableLimits.gms' // Set new variable limits (.lo and .up)
$$include 'inc\3e_solve.gms' // Solve model(s)
$$include 'inc\4a_outputVariant.gms' // Store results from the loop
$$ifi '%debug%' == 'yes' execute_unload 'output\debug.gdx'; // Output debugging information
* $$ifi.debug '%debug%' == 'yes'
* putclose gdx;
* put_utility 'gdxout' / 'output\'mSolve.tl:0, '-', tSolve.tl:0, '.gdx';
......
* Generate model rules from basic patterns defined in the model definition files
* Set the time for the next available forecast unless it has been preset.
loop(m,
tForecastNext(m) = mSettings(m, 't_forecastStart');
);
* Check the modelSolves for preset patterns for model solve timings
* If not found, then use mSettings to set the model solve timings
loop(m,
......@@ -116,7 +121,6 @@ loop(unit,
+ 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') ]
};
put log unit.tl:20, heat_rate /;
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')
......
......@@ -6,6 +6,8 @@
p_stepLength(mSolve, f, t) = no;
ft_new(f,t) = no;
tForecastNext(mSolve)$(ord(tSolve) >= tForecastNext(mSolve)) = tForecastNext(mSolve) + mSettings(mSolve, 't_ForecastJump');
$offOrder
// Calculate parameters affected by model intervals
if(sum(counter, mInterval(mSolve, 'intervalLength', counter)), // But only if interval data has been set in the model definition file
......@@ -16,7 +18,7 @@ $offOrder
tInterval(t)$( ord(t) >= tSolveFirst + tCounter
and ord(t) < min(tSolveFirst + mInterval(mSolve, 'intervalEnd', counter), tSolveLast)
) = yes; // Set of t's where the interval is 1 (right border defined by intervalEnd) - but do not go beyond tSolveLast
ts_influx_(grid, node, fSolve, t)$tInterval(t) = ts_influx(grid, node, fSolve, t+ct(t));
ts_influx_(grid, node, fSolve, t)$tInterval(t) = ts_influx(grid, node, 'f00', t+ct(t));
ts_cf_(flow, node, fSolve, t)$tInterval(t) = ts_cf(flow, node, fSolve, t+ct(t));
ts_nodeState_(gn_state(grid, node), param_gnBoundaryTypes, fSolve, t)$tInterval(t) = ts_nodeState(grid, node, param_gnBoundaryTypes, fSolve, t+ct(t));
ts_unit_(unit, param_unit, fSolve, t)$tInterval(t) = ts_unit(unit, param_unit, fSolve, t+ct(t));
......@@ -24,7 +26,6 @@ $offOrder
p_stepLength(mf(mSolve, fSolve), t)$tInterval(t) = mSettings(mSolve, 'intervalInHours'); // p_stepLength will hold the length of the interval in hours in model equations
p_stepLengthNoReset(mSolve, fSolve, t)$tInterval(t) = mSettings(mSolve, 'intervalInHours');
pt(t + 1)$tInterval(t) = -1;
mftLastSteps(mf(mSolve,fSolve),t)$[ord(t)-1 = mInterval(mSolve,'intervalEnd',counter)] = yes;
elseif mInterval(mSolve, 'intervalLength', counter) > 1, // intervalLength > 1 (not defined if intervalLength < 1)
loop(t$[ord(t) >= tSolveFirst + tCounter and ord(t) < min(tSolveFirst + mInterval(mSolve, 'intervalEnd', counter), tSolveLast)], // Loop t relevant for each interval (interval borders defined by intervalEnd) - but do not go beyond tSolveLast
if (not mod(tCounter - mInterval(mSolve, 'intervalEnd', counter), mInterval(mSolve, 'intervalLength', counter)), // Skip those t's that are not at the start of any interval
......@@ -45,14 +46,6 @@ $offOrder
// Set the previous time step displacement
pt(t+intervalLength) = -intervalLength;
);
if (mInterval(mSolve, 'intervalEnd', counter) <= mSettings(mSolve, 't_forecastLength'),
mftLastForecast(mf(mSolve,fSolve),t_) = no;
mftLastForecast(mf(mSolve,fSolve),t)$[ord(t) = tSolveFirst + tCounter] = yes;
);
if (mInterval(mSolve, 'intervalEnd', counter) <= tSolveLast,
mftLastSteps(mf(mSolve,fSolve),t_) = no;
mftLastSteps(mf(mSolve,fSolve),t+intervalLength)$[ord(t) = tSolveFirst + tCounter] = yes;
);
); // end if that skips t's not at the start of any interval
tCounter = tCounter + 1;
); // end loop t
......@@ -69,92 +62,27 @@ $offOrder
$$ifi '%rampSched%' == 'yes' ts_unit_(unit, param_unit, fSolve, t)${ord(t) = tSolveLast} = ts_unit(unit, param_unit, fSolve, t+ct(t));
$onOrder
// Set mft for the modelling period and model forecasts
mft(mSolve,f,t) = no;
mft(mSolve, f, t)$( p_stepLength(mSolve, f, t) and ord(t) < tSolveFirst + mSettings(mSolve, 't_horizon' ) ) = yes;
* mft(mSolve, f, t)${ [ord(t) >= ord(tSolve)]
* $$ifi '%rampSched%' == 'yes' and [ord(t) <=
* $$ifi not '%rampSched%' == 'yes' and [ord(t) <
* ord(tSolve) + mSettings(mSolve, 't_forecastLength')]
* and mf(mSolve, f)
* } = yes;
mftStart(mSolve,f,t) = no;
mftStart(mSolve,fRealization,t)$[ord(t) = ord(tSolve)] = yes;
mftBind(mSolve,f,t) = no;
mft_bind(mSolve,f,t) = no;
mt_bind(mSolve,t) = no;
* mftBind(mft(mSolve,f,t))$[ord(t) = ord(tSolve) + mSettings(mSolve, 't_forecastLength')] = yes;
* mft_bind(mft(mSolve,f,t))$[ord(t) = ord(tSolve) + mSettings(mSolve, 't_forecastLength')] = 1 - ord(f);
* mt_bind(mSolve,t)$[ord(t) = ord(tSolve) + mSettings(mSolve, 't_forecastLength')] = -1;
msft(mSolve, s, f, t) = no;
msft(mSolve, 's000', f, t) = mft(mSolve,f,t);
msft(mSolve, 's000', fRealization(f), t)${ [ord(t) >= ord(tSolve) + mSettings(mSolve, 't_forecastLength')]
$$ifi '%rampSched%' == 'yes' and [ord(t) <=
$$ifi not '%rampSched%' == 'yes' and [ord(t) <
ord(tSolve) + mSettings(mSolve, 't_horizon')]
and mf(mSolve, f)
} = yes;
ft(f,t) = no;
ft(f,t) = mft(mSolve, f, t);
ft_dynamic(f,t) = ft(f,t);
ft_dynamic(f,tSolve) = no;
ft_full(f,t) = no;
ft_full(f,t) = ft(f,t);
loop(counter$mInterval(mSolve, 'intervalLength', counter),
lastCounter = ord(counter);
);
loop(counter$(ord(counter) = lastCounter),
ft_dynamic(f,t)$(mf(mSolve, f) and ord(t) = min(tSolveFirst + mInterval(mSolve, 'intervalEnd', counter), tSolveLast)) = yes;
ft_full(f,t)$ (mf(mSolve, f) and ord(t) = min(tSolveFirst + mInterval(mSolve, 'intervalEnd', counter), tSolveLast)) = yes;
);
ft_realized(f,t) = no;
ft_realized(f,t)$[fRealization(f) and ord(t) >= ord(tSolve) and ord(t) <= ord(tSolve) + mSettings(mSolve, 't_jump')] = yes;
ft_realizedLast(f,t) = no;
ft_realizedLast(f,t)$[fRealization(f) and ord(t) = ord(tSolve) + mSettings(mSolve, 't_jump')] = yes;
// Defining unit aggregations and ramps
uft(unit, f, t) = no;
uft(unit, f, t)$[ ft(f, t)
and ord(t) <= tSolveFirst + mSettings(mSolve, 't_aggregate') - 1
and not unit_aggregate(unit)
] = yes;
uft(unit, f, t)$[ ft(f, t)
and ord(t) > tSolveFirst + mSettings(mSolve, 't_aggregate') - 1
and (unit_aggregate(unit) or unit_noAggregate(unit))
] = yes;
nuft(node, unit, f, t) = no;
nuft(node, unit, f, t)$(nu(node, unit) and uft(unit, f, t)) = yes;
gnuft(grid, node, unit, f, t) = no;
gnuft_ramp(grid, node, unit, f, t) = no;
gnuft(grid, node, unit, f, t)$(gn(grid, node) and nuft(node, unit, f, t)) = yes;
gnuft_ramp(gnuft(grid, node, unit, f, t))${ p_gnu(grid, node, unit, 'maxRampUp')
OR p_gnu(grid, node, unit, 'maxRampDown')
OR p_gnu(grid, node, unit, 'rampUpCost')
OR p_gnu(grid, node, unit, 'rampDownCost') } = yes;
// Defining unit efficiency groups etc.
suft(effGroup, unit, f, t) = no;
loop(effLevel$mSettingsEff(mSolve, effLevel),
tInterval(t) = no;
tInterval(t)$(ord(effLevel) = 1 and ord(t) = tSolveFirst) = yes;
tInterval(t)$( ord(t) >= tSolveFirst + mSettingsEff(mSolve, effLevel)
and ord(t) < tSolveFirst + mSettingsEff(mSolve, effLevel+1)
) = yes;
loop(effLevelGroupUnit(effLevel, effGroup, unit)$(not sum(flow$flowUnit(flow, unit), 1)),
suft(effGroup, unit, f, tInterval(t))$(effLevelGroupUnit(effLevel, effGroup, unit) and uft(unit, f, tInterval)) = yes;
);
);
sufts(effGroup, unit, f, t, effSelector) = no;
sufts(effGroup, unit, f, t, effSelector)$(effGroupSelector(effGroup, effSelector) and suft(effGroup, unit, f, t)) = yes;
uft_online(unit, f, t) = no; // Initialize the 'uft_online' set so that no units use online variables by default.
loop(suft(effOnline, uft(unit, f, t)), // Determine the time steps when units need to have online variables.
uft_online(unit, f, t) = yes;
);
$include 'defModels\periodicLoopDispatch.gms';
ft_limits(f,t) = no;
ft_limits(f,t) = ft_full(f,t) + ft_realized(f,t);
uft_limits(unit, f, t) = no;
uft_limits(unit, f, t)$[ ft_limits(f, t)
and ord(t) <= tSolveFirst + mSettings(mSolve, 't_aggregate') - 1
and not unit_aggregate(unit)
] = yes;
uft_limits(unit, f, t)$[ ft_limits(f, t)
and ord(t) > tSolveFirst + mSettings(mSolve, 't_aggregate') - 1
and (unit_aggregate(unit) or unit_noAggregate(unit))
] = yes;
uft_limits_online(unit, f, t) = no;
loop(suft(effOnline, uft_limits(unit, f, t)), // Determine the time steps when units need to have online variables.
uft_limits_online(unit, f, t) = yes;
);
pf(ft(f,t))$(ord(t) eq ord(tSolve) + 1) = 1 - ord(f);
// Calculate time series for unit parameters when necessary and/or possible
loop(unit${p_unit(unit, 'useTimeseries')},
......
// Set mft for the modelling period and model forecasts
mft(mSolve,f,t) = no;
mft(mSolve, f, t)$( p_stepLength(mSolve, f, t) and ord(t) < tSolveFirst + mSettings(mSolve, 't_horizon') and
{
[ ord(f) = 1 and ord(t) = tSolveFirst + tDispatchCurrent - 1]
or
[ ord(f) > 1 and ord(t) <= tSolveFirst + mSettings(mSolve, 't_forecastLength') and ord(t) >= tSolveFirst + tDispatchCurrent ]
or
[ fCentral(f) and ord(t) > mSettings(mSolve, 't_forecastLength') ]
} ) = yes;
* mft(mSolve, f, t)${ [ord(t) >= ord(tSolve)]
* $$ifi '%rampSched%' == 'yes' and [ord(t) <=
* $$ifi not '%rampSched%' == 'yes' and [ord(t) <
* ord(tSolve) + mSettings(mSolve, 't_forecastLength')]
* and mf(mSolve, f)
* } = yes;
mftStart(mSolve,f,t) = no;
mftStart(mSolve,f,t)${ tDispatchCurrent = 0 and fSolve(f) and ord(f) > 1 and ord(t) = tSolveFirst + tDispatchCurrent} = yes;
mftStart(mSolve,f,t)${ tDispatchCurrent >= 1 and ord(f) = 1 and ord(t) = tSolveFirst + tDispatchCurrent - 1} = yes;
//mftStart(mSolve,fRealization,t)$[ord(t) = ord(tSolve)] = yes;
mftBind(mSolve,f,t) = no;
mft_bind(mSolve,f,t) = no;
mt_bind(mSolve,t) = no;
// Connect the realization of the dispatch to the branches of the forecast tree
mftBind(mft(mSolve,f,t))$[ord(t) = tSolveFirst + tDispatchCurrent] = yes;
mft_bind(mft(mSolve,f,t))$[ord(t) = tSolveFirst + tDispatchCurrent] = 1 - ord(f);
mt_bind(mSolve,t)$[ord(t) = tSolveFirst + tDispatchCurrent] = 1;
// Connect branches of the forecast tree to the central forecast at the end of the forecast period
mftBind(mft(mSolve,f,t))$[ord(t) = ord(tSolve) + mSettings(mSolve, 't_forecastLength')] = yes;
mft_bind(mft(mSolve,f,t))$[ord(t) = ord(tSolve) + mSettings(mSolve, 't_forecastLength')] = sum(f_$(fCentral(f_)), ord(f_)) - ord(f);
mt_bind(mSolve,t)$[ord(t) = ord(tSolve) + mSettings(mSolve, 't_forecastLength')] = -1;
msft(mSolve, s, f, t) = no;
msft(mSolve, 's000', f, t) = mft(mSolve,f,t);
* msft(mSolve, 's000', fRealization(f), t)${ [ord(t) >= ord(tSolve) + mSettings(mSolve, 't_forecastLength')]
* $$ifi '%rampSched%' == 'yes' and [ord(t) <=
* $$ifi not '%rampSched%' == 'yes' and [ord(t) <
* ord(tSolve) + mSettings(mSolve, 't_horizon')]
* and mf(mSolve, f)
* } = yes;
ft(f,t) = no;
ft(f,t) = mft(mSolve, f, t);
ft_dynamic(f,t) = no;
ft_dynamic(f,t) = ft(f,t);
//fft_dynamic(f,f_,t) = no;
//fft_dynamic(f,f_,t)$(ord(f) = ord(f_) and ft(f,t) and [ord(t) < tSolveFirst + mSettings(mSolve, 't_forecastLength') + 1 or ord(t) > tSolveFirst + mSettings(mSolve, 't_forecastLength') + 1]) = yes;
//fft_dynamic(f,f_,t)$(ord(f_) > 1 and mf(mSolve, f_) and ft(f,t) and ord(t) = tSolveFirst + mSettings(mSolve, 't_forecastLength') + 1) = yes;
//fft_dynamic(f,f_,tSolve+tDispatchCurrent) = no;
continueLoop = 1;
cf(f,t) = no;
loop(t$(ord(t) > tSolveFirst + mSettings(mSolve,'t_forecastLength') and continueLoop), // Loop through all the different intervals set in the model definition file
if(sum(f,p_stepLength(mSolve,f,t)),
ft_dynamic(f,t)$(ord(f) > 1) = yes;
cf(ft_dynamic(f,t)) = sum(f_$(fCentral(f_)), ord(f_)) - ord(f);
continueLoop = 0;
);
);
ft_dynamic(f,tSolve+max(0,(tDispatchCurrent -1))) = no;
ft_dynamic('f00',tSolve+tDispatchCurrent)$(tDispatchCurrent > 0) = yes;
ft_dynamic(f,t)$(ord(f) = sum(f_$(fCentral(f_)), ord(f_)) and ord(t) = tSolveFirst + mSettings(mSolve, 't_horizon')) = yes;
ft_full(f,t) = no;
ft_full(f,t) = ft(f,t) + ft_dynamic(f,t);
ft_realized(f,t) = no;
ft_realized(f,t)$[fRealization(f) and ord(t) >= ord(tSolve) and ord(t) <= ord(tSolve) + mSettings(mSolve, 't_jump')] = yes;
ft_realizedLast(f,t) = no;
ft_realizedLast(f,t)$[fRealization(f) and ord(t) = ord(tSolve) + mSettings(mSolve, 't_jump')] = yes;
// Defining unit aggregations and ramps
uft(unit, f, t) = no;
uft(unit, f, t)$[ ft(f, t)
and ord(t) <= tSolveFirst + mSettings(mSolve, 't_aggregate') - 1
and not unit_aggregate(unit)
] = yes;
uft(unit, f, t)$[ ft(f, t)
and ord(t) > tSolveFirst + mSettings(mSolve, 't_aggregate') - 1
and (unit_aggregate(unit) or unit_noAggregate(unit))
] = yes;
uft_limits(unit, f, t) = no;
nuft(node, unit, f, t) = no;
nuft(node, unit, f, t)$(nu(node, unit) and uft(unit, f, t)) = yes;
gnuft(grid, node, unit, f, t) = no;
gnuft_ramp(grid, node, unit, f, t) = no;
gnuft(grid, node, unit, f, t)$(gn(grid, node) and nuft(node, unit, f, t)) = yes;
gnuft_ramp(gnuft(grid, node, unit, f, t))${ p_gnu(grid, node, unit, 'maxRampUp')
OR p_gnu(grid, node, unit, 'maxRampDown')
OR p_gnu(grid, node, unit, 'rampUpCost')
OR p_gnu(grid, node, unit, 'rampDownCost') } = yes;
// Defining unit efficiency groups etc.
suft(effGroup, unit, f, t) = no;
loop(effLevel$mSettingsEff(mSolve, effLevel),
tInterval(t) = no;
tInterval(t)$(ord(effLevel) = 1 and ord(t) = tSolveFirst) = yes;
tInterval(t)$( ord(t) >= tSolveFirst + mSettingsEff(mSolve, effLevel)
and ord(t) < tSolveFirst + mSettingsEff(mSolve, effLevel+1)
) = yes;
loop(effLevelGroupUnit(effLevel, effGroup, unit)$(not sum(flow$flowUnit(flow, unit), 1)),
suft(effGroup, unit, f, tInterval(t))$(effLevelGroupUnit(effLevel, effGroup, unit) and uft(unit, f, tInterval)) = yes;
);
);
sufts(effGroup, unit, f, t, effSelector) = no;
sufts(effGroup, unit, f, t, effSelector)$(effGroupSelector(effGroup, effSelector) and suft(effGroup, unit, f, t)) = yes;
uft_online(unit, f, t) = no; // Initialize the 'uft_online' set so that no units use online variables by default.
loop(suft(effOnline, uft(unit, f, t)), // Determine the time steps when units need to have online variables.
uft_online(unit, f, t) = yes;
);
pf(f,t) = no;
pf(ft(f,t))$(ord(t) eq ord(tSolve) + tDispatchCurrent) = 1 - ord(f);
$offOrder
loop(counter$mInterval(mSolve, 'intervalLength', counter), // Loop through all the different intervals set in the model definition file
if (mSettings(mSolve, 't_horizon') > mInterval(mSolve, 'intervalEnd', counter-1) and mSettings(mSolve, 't_horizon') <= mInterval(mSolve, 'intervalEnd', counter) ,
mftLastSteps(mf(mSolve,fSolve),t) = no;
if (mInterval(mSolve, 'intervalEnd', counter) <= mSettings(mSolve, 't_forecastLength'),
mftLastSteps(mf(mSolve,fSolve),t)$[ord(fSolve) > 1 and ord(t) = tSolveFirst + mSettings(mSolve, 't_horizon')] = yes;
mftLastForecast(mf(mSolve,fSolve),t) = no;
mftLastForecast(mf(mSolve,fSolve),t)$[ord(fSolve) > 1 and ord(t) = tSolveFirst + mSettings(mSolve, 't_forecastLength')] = yes;
else
mftLastSteps(mf(mSolve,fSolve),t)$[fCentral(fSolve) and ord(t) = tSolveFirst + mSettings(mSolve, 't_horizon')] = yes;
);
);
);
$onOrder;
......@@ -4,7 +4,7 @@ if (mType('schedule'),
// Define the temporal structure of the model in time indeces
mSettings('schedule', 'intervalInHours') = 1; // Define the duration of a single time-step in hours
mInterval('schedule', 'intervalLength', 'c000') = 1;
mInterval('schedule', 'intervalEnd', 'c000') = 168;
mInterval('schedule', 'intervalEnd', 'c000') = 48;
mInterval('schedule', 'intervalLength', 'c001') = 3;
mInterval('schedule', 'intervalEnd', 'c001') = 336;
mInterval('schedule', 'intervalLength', 'c002') = 6;
......@@ -16,29 +16,36 @@ if (mType('schedule'),
// Define the model execution parameters in time indeces
mSettings('schedule', 't_start') = 1; // Ord of first solve (i.e. >0)
mSettings('schedule', 't_horizon') = 8760;
mSettings('schedule', 't_jump') = 48;
mSettings('schedule', 't_forecastLength') = 8760;
mSettings('schedule', 't_end') = 80;
mSettings('schedule', 't_horizon') = 168;
mSettings('schedule', 't_jump') = 12;
mSettings('schedule', 't_forecastStart') = 1; // Ord of first forecast available
mSettings('schedule', 't_forecastLength') = 72;
mSettings('schedule', 't_forecastJump') = 24;
mSettings('schedule', 't_end') = 30;
// Define unit aggregation and efficiency levels starting indeces
mSettings('schedule', 't_aggregate') = 72;
mSettingsEff('schedule', 'level1') = 1;
mSettingsEff('schedule', 'level2') = 24;
mSettingsEff('schedule', 'level3') = 48;
mSettingsEff('schedule', 'level4') = 168;
mSettingsEff('schedule', 'level2') = 6;
mSettingsEff('schedule', 'level3') = 12;
mSettingsEff('schedule', 'level4') = 18;
//mSettingsEff('schedule', 'level1') = 1;
//mSettingsEff('schedule', 'level2') = 24;
//mSettingsEff('schedule', 'level3') = 48;
//mSettingsEff('schedule', 'level4') = 168;
// Define active model features
active('storageValue') = yes;
// Define model stochastic parameters
mSettings('schedule', 'samples') = 1;
mSettings('schedule', 'forecasts') = 0;
mSettings('schedule', 'forecasts') = 3;
mSettings('schedule', 'readForecastsInTheLoop') = 1;
mf('schedule', f)$[ord(f)-1 <= mSettings('schedule', 'forecasts')] = yes;
fRealization(f) = no;
fRealization('f00') = yes;
fCentral(f) = no;
fCentral('f01') = yes;
fCentral('f02') = yes;
sInitial(s) = no;
sInitial('s000') = yes;
sCentral(s) = no;
......@@ -49,6 +56,9 @@ if (mType('schedule'),
p_sProbability('s000') = 1;
p_fProbability(f) = 0;
p_fProbability(fRealization) = 1;
p_fProbability('f01') = 0.2;
p_fProbability('f02') = 0.6;
p_fProbability('f03') = 0.2;
);
Model schedule /
......@@ -61,7 +71,6 @@ Model schedule /
q_startup
q_genRamp
q_genRampChange
q_bindOnline
q_conversionDirectInputOutput
q_conversionSOS2InputIntermediate
q_conversionSOS2Constraint
......@@ -72,7 +81,32 @@ Model schedule /
q_stateUpwardLimit
q_stateDownwardLimit
q_boundState
q_boundStateMaxDiff
q_boundCyclic
q_bidirectionalTransfer
/;
Model schedule_dispatch /
q_obj
q_balance
q_resDemand
q_resTransfer
q_maxDownward
q_maxUpward
q_startup
q_genRamp
q_genRampChange
q_conversionDirectInputOutput
q_conversionSOS2InputIntermediate
q_conversionSOS2Constraint
q_conversionSOS2IntermediateOutput
q_outputRatioFixed
q_outputRatioConstrained
q_stateSlack
q_stateUpwardLimit
q_stateDownwardLimit
q_boundState
q_boundCyclic
q_bidirectionalTransfer
/;
......@@ -3,7 +3,7 @@ Sets // Model related selections
mType "model types in the Backbone"
/invest, storage, schedule, realtime, building/
mSetting "setting categories for models"
/t_start, t_jump, t_horizon, t_forecastLength, t_end, samples, forecasts, intervalEnd, intervalLength, IntervalInHours, t_aggregate/
/t_start, t_jump, t_horizon, t_forecastStart, t_forecastLength, t_forecastJump, t_end, samples, forecasts, readForecastsInTheLoop, intervalEnd, intervalLength, IntervalInHours, t_aggregate/
counter "general counter set"
/c000*c999/
......@@ -88,6 +88,7 @@ param_gn "Possible parameters for grid, node" /
boundAll "A flag to bound the state to the reference in all time steps"
boundStartToEnd "Force the last states to equal the first state"
boundCyclic "A flag to impose cyclic bounds for the first and the last states"
forecastLength "Length of forecasts in use for the node (hours). After this, the node will use the central forecast."
/
param_gnBoundaryTypes "Types of boundaries that can be set for a node with a state variable" /
......
......@@ -73,9 +73,11 @@ Sets
mstStart(mType, s, t) "Start point of samples"
ft(f, t) "Combination of forecasts and time periods in the current model"
ft_dynamic(f, t) "ft without first t and with tLast+1 (moved right)"
//fft_dynamic(f, f, t) "ft without first t and with tLast+1 (moved right) and additional forecasts when decreasing the number of scenarios in the forecast tree"
ft_full(f, t) "ft with all t's in the solve including tSolve and tLast+1"
ft_realized(f, t) "Realized ft"
ft_realizedLast(f, t) "Last realized ft"
ft_limits(f, t) "All ft for which variable limits are set within the tSolve loop"
ft_new(f, t) "Newly introduced f,t to be used in calculating parameter/variable values"
mft(mType, f, t) "Combination of forecasts and time periods in the models"
mft_(mType, f, t) "Combination of forecasts and time periods in the models"
......@@ -90,10 +92,13 @@ Sets
mftLastSteps(mType, f, t) "Last time periods of the model (can be end of forecasts or end of samples)"
modelSolves(mType, t) "when different models are to be solved"
fSolve(f) "forecasts in the model to be solved next"
tSolveDispatch(t)
* --- Sets used for the changing unit aggregation and efficiency approximations
uft(unit, f, t) "Enables aggregation of units for later time periods"
uft_online(unit, f, t) "Units with online and startup variables on time periods"
uft_limits(unit, f, t) "Used to set limits in the tSolve loop"
uft_limits_online(unit, f, t) "Used to set limits in the tSolve loop"
nuft(node, unit, f, t) "Enables aggregation of nodes and units for later time periods"
gnuft(grid, node, unit, f, t) "Enables aggregation of nodes and units for later time periods"
gnuft_ramp(grid, node, unit, f, t) "Units with ramp requirements or costs"
......
......@@ -3,6 +3,7 @@ Scalars
errorcount /0/
tSolveFirst "counter (ord) for the first t in the solve"
tSolveLast "counter for the last t in the solve"
tDispatchCurrent "counter for the current t in the dispatch loop" /0/
tCounter "counter for t" /0/
lastCounter "last member in use of the general counter"
ts_length "Length of time series (t)"
......@@ -53,12 +54,14 @@ Scalar p_sWeightSum "Sum of sample weights";
Parameters
pt(t) "Displacement needed to reach the previous time period (in time periods)"
pf(f, t) "Displacement needed to reach the previous forecast (in forecasts)"
cf(f, t) "Displacement needed to reach the current forecast (in forecasts) - this is needed when the forecast tree gets reduced in dynamic equations"
ct(t) "Circular t displacement if the time series data is not long enough to cover the model horizon"
t_bind(t) "Displacement to reach the binding time period in the parent sample (in time periods). Can skip with aggregated steps as well as when connecting samples."
ft_bind(f, t) "Displacement to reach the binding forecast (in forecasts) in the current model"
mt_bind(mType, t) "Displacement to reach the binding time period in the parent sample (in time periods) in the models"
mft_bind(mType, f, t) "Displacement to reach the binding forecast (in forecasts) in the models"
p_slackDirection(slack) "+1 for upward slacks and -1 for downward slacks"
tForecastNext(mType) "When the next forecast will be availalbe (ord time)"
;
* --- Stochastic data parameters ----------------------------------------------
......@@ -75,6 +78,7 @@ Parameters
ts_influx_(grid, node, f, t)
ts_cf_(flow, node, f, t)
ts_nodeState_(grid, node, param_gnBoundaryTypes, f, t)
ts_forecast(flow, node, t, f, t)
;
* --- Other time dependent parameters -----------------------------------------
......
......@@ -5,7 +5,6 @@ equations
q_resTransfer(grid, node, node, f, t) "Transfer of energy and capacity reservations are less than the transfer capacity"
q_maxDownward(grid, node, unit, f, t) "Downward commitments will not undercut power plant minimum load constraints or maximum elec. consumption"
q_maxUpward(grid, node, unit, f, t) "Upward commitments will not exceed maximum available capacity or consumed power"
q_bindOnline(unit, mType, f, t) "Couple online variable when joining forecasts or when joining sample time periods"
q_startup(unit, f, t) "Capacity started up is greater than the difference of online cap. now and in the previous time step"
q_genRamp(grid, node, mType, unit, f, t) "Record the ramps of units with ramp restricitions or costs"
q_genRampChange(grid, node, mType, unit, f, t) "Record the ramp rates of units with ramping costs"
......@@ -139,7 +138,7 @@ q_obj ..
q_balance(gn(grid, node), m, ft_dynamic(f, t))${p_stepLength(m, f+pf(f,t), t+pt(t)) and not p_gn(grid, node, 'boundAll')} .. // Energy/power balance dynamics solved using implicit Euler discretization
// The left side of the equation is the change in the state (will be zero if the node doesn't have a state)
+ p_gn(grid, node, 'energyStoredPerUnitOfState')$gn_state(grid, node) // Unit conversion between v_state of a particular node and energy variables (defaults to 1, but can have node based values if e.g. v_state is in Kelvins and each node has a different heat storage capacity)
* ( + v_state(grid, node, f, t) // The difference between current
* ( + v_state(grid, node, f+cf(f,t), t) // The difference between current
- v_state(grid, node, f+pf(f,t), t+pt(t)) // ... and previous state of the node
)
=E=
......@@ -148,20 +147,20 @@ q_balance(gn(grid, node), m, ft_dynamic(f, t))${p_stepLength(m, f+pf(f,t), t+pt(
+ (
// Self discharge out of the model boundaries
- p_gn(grid, node, 'selfDischargeLoss')$gn_state(grid, node) * (
+ v_state(grid, node, f, t) // The current state of the node
+ v_state(grid, node, f+cf(f,t), t) // The current state of the node
$$ifi '%rampSched%' == 'yes' + v_state(grid, node, f+pf(f,t), t+pt(t)) // and possibly averaging with the previous state of the node
)
// Energy diffusion from this node to neighbouring nodes
- sum(to_node$(gnn_state(grid, node, to_node)),
+ p_gnn(grid, node, to_node, 'diffCoeff') * (
+ v_state(grid, node, f, t)
+ v_state(grid, node, f+cf(f,t), t)
$$ifi '%rampSched%' == 'yes' + v_state(grid, node, f+pf(f,t), t+pt(t))
)
)
// Energy diffusion from neighbouring nodes to this node
+ sum(from_node$(gnn_state(grid, from_node, node)),
+ p_gnn(grid, from_node, node, 'diffCoeff') * (
+ v_state(grid, from_node, f, t) // Incoming diffusion based on the state of the neighbouring node
+ v_state(grid, from_node, f+cf(f,t), t) // Incoming diffusion based on the state of the neighbouring node
$$ifi '%rampSched%' == 'yes' + v_state(grid, from_node, f+pf(f,t), t+pt(t)) // Ramp schedule averaging, NOTE! State and other terms use different indeces for non-ramp-schedule!
)
)
......@@ -169,31 +168,31 @@ q_balance(gn(grid, node), m, ft_dynamic(f, t))${p_stepLength(m, f+pf(f,t), t+pt(
+ sum(from_node$(gn2n(grid, from_node, node)),
+ (1 - p_gnn(grid, from_node, node, 'transferLoss')) * ( // Include transfer losses
+ v_transfer(grid, from_node, node, f+pf(f,t), t+pt(t))
$$ifi '%rampSched%' == 'yes' + v_transfer(grid, from_node, node, f, t) // Ramp schedule averaging, NOTE! State and other terms use different indeces for non-ramp-schedule!
$$ifi '%rampSched%' == 'yes' + v_transfer(grid, from_node, node, f+cf(f,t), t) // Ramp schedule averaging, NOTE! State and other terms use different indeces for non-ramp-schedule!
)
)
// Controlled energy transfer to other nodes from this one
- sum(to_node$(gn2n(grid, node, to_node)),
+ v_transfer(grid, node, to_node, f+pf(f,t), t+pt(t)) // Transfer losses accounted for in the previous term
$$ifi '%rampSched%' == 'yes' + v_transfer(grid, node, to_node, f, t) // Ramp schedule averaging
$$ifi '%rampSched%' == 'yes' + v_transfer(grid, node, to_node, f+cf(f,t), t) // Ramp schedule averaging
)
// Interactions between the node and its units
+ sum(gnuft(grid, node, unit, f+pf(f,t), t+pt(t)),
+ v_gen(grid, node, unit, f+pf(f,t), t+pt(t)) // Unit energy generation and consumption
$$ifi '%rampSched%' == 'yes' + v_gen(grid, node, unit, f, t)
$$ifi '%rampSched%' == 'yes' + v_gen(grid, node, unit, f+cf(f,t), t)
)
// Spilling energy out of the endogenous grids in the model
- v_spill(grid, node, f+pf(f,t), t+pt(t))$node_spill(node)
$$ifi '%rampSched%' == 'yes' - v_spill(grid, node, f, t)$node_spill(node)
$$ifi '%rampSched%' == 'yes' - v_spill(grid, node, f+cf(f,t), t)$node_spill(node)
// Power inflow and outflow timeseries to/from the node
+ ts_influx_(grid, node, f+pf(f,t), t+pt(t)) // Incoming (positive) and outgoing (negative) absolute value time series
$$ifi '%rampSched%' == 'yes' + ts_influx_(grid, node, f, t)
$$ifi '%rampSched%' == 'yes' + ts_influx_(grid, node, f+cf(f,t), t)
// Dummy generation variables, for feasibility purposes
+ vq_gen('increase', grid, node, f+pf(f,t), t+pt(t)) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
$$ifi '%rampSched%' == 'yes' + vq_gen('increase', grid, node, f, t)
$$ifi '%rampSched%' == 'yes' + vq_gen('increase', grid, node, f+cf(f,t), t)
- vq_gen('decrease', grid, node, f+pf(f,t), t+pt(t)) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
$$ifi '%rampSched%' == 'yes' - vq_gen('decrease', grid, node, f, t)
$$ifi '%rampSched%' == 'yes' - vq_gen('decrease', grid, node, f+cf(f,t), t)
) * p_stepLength(m, f+pf(f,t), t+pt(t)) // Multiply by time step to get energy terms
)
$$ifi '%rampSched%' == 'yes' / 2 // Averaging all the terms on the right side of the equation over the timestep here.
......@@ -272,7 +271,7 @@ q_maxUpward(gnuft(grid, node, unit, f, t))${ [uft_online(unit, f, t) and p_
q_startup(uft_online(unit, ft_dynamic(f, t))) ..
+ v_startup(unit, f+pf(f,t), t+pt(t))
=G=
+ v_online(unit, f, t)
+ v_online(unit, f+cf(f,t), t)
- v_online(unit, f+pf(f,t), t+pt(t)) // This reaches to tFirstSolve when pt = -1
;
* -----------------------------------------------------------------------------
......@@ -282,32 +281,26 @@ q_genRamp(gn(grid, node), m, unit, ft_dynamic(f, t))${gnuft_ramp(grid, node, uni
* p_stepLength(m, f+pf(f,t), t+pt(t))
=E=
// Change in generation over the time step
+ v_gen(grid, node, unit, f, t)
+ v_gen(grid, node, unit, f+cf(f,t), t)
- v_gen(grid, node, unit, f+pf(f,t), t+pt(t))
// Correction term to account for online variables and min loads
- (
+ v_online(unit, f, t)${uft_online(unit, f, t)}
+ v_online(unit, f+cf(f,t), t)${uft_online(unit, f, t)}
- v_online(unit, f+pf(f,t), t+pt(t))${uft_online(unit, f+pf(f,t), t+pt(t))}
)
/ p_unit(unit, 'unitCount')
* ( p_gnu(grid, node, unit, 'maxGen') - p_gnu(grid, node, unit, 'maxCons') )
* sum(suft(effGroup, unit, f, t), p_effGroupUnit(effGroup, unit, 'lb')) // Newly started units are assumed to start to their minload.
* sum(suft(effGroup, unit, f+cf(f,t), t), p_effGroupUnit(effGroup, unit, 'lb')) // Newly started units are assumed to start to their minload.
;
* -----------------------------------------------------------------------------
q_genRampChange(gn(grid, node), m, unit, ft_dynamic(f, t))${ gnuft_ramp(grid, node, unit, f, t) AND [ p_gnu(grid, node, unit, 'rampUpCost') OR p_gnu(grid, node, unit, 'rampDownCost') ]} ..
+ v_genRampChange(grid, node, unit, 'up', f, t)
- v_genRampChange(grid, node, unit, 'down', f, t)
+ v_genRampChange(grid, node, unit, 'up', f+cf(f,t), t)
- v_genRampChange(grid, node, unit, 'down', f+cf(f,t), t)
=E=
+ v_genRamp(grid, node, unit, f, t)
+ v_genRamp(grid, node, unit, f+cf(f,t), t)
- v_genRamp(grid, node, unit, f+pf(f,t), t+pt(t))
;
* -----------------------------------------------------------------------------
q_bindOnline(unit, mftBind(m, f, t))${uft_online(unit, f, t) and uft_online(unit, f+mft_bind(m,f,t), t+mt_bind(m,t))} ..
+ v_online(unit, f, t)
=E=