Commit fc4e9380 authored by Erkka Rinne's avatar Erkka Rinne
Browse files

Reduce scenarios based on total available energy

Added new parameter `ts_energy_(s)` to which total energy is
calculated. This parameter has the sample dimension last, as required
by the SCENRED2 tool. This allowed changing all other time series
to have indices `s, f, t` instead of previous `f, t, s` and dropping
the set `fts` completely.

Fixes #105
parent 7efcad0b
......@@ -12,6 +12,7 @@ All notable changes to this project will be documented in this file.
### Changed
- Updated tool defintions for Sceleton Titan and Spine Toolbox
- The program will now stop looping in case of execution errors.
- Scenario reduction is done based on total available energy
### Fixed
- Removed hard-coded `elec grids` from *setVariableLimits* and *rampSched files*
......
......@@ -110,7 +110,6 @@ Sets
msf(mType, s, f) "Combination of samples and forecasts in the models"
msft(mType, s, f, t) "Combination of models, samples, forecasts and t's"
sft(s, f, t) "Combination of samples, forecasts and t's in the current model solve"
fts(f, t, s) "Like `sft` but different order"
sft_realized(s, f, t) "Realized sft"
sft_realizedNoReset(s, f, t) "Full set of realized sft, facilitates calculation of results"
msft_realizedNoReset(mType, s, f, t) "Combination of realized samples, forecasts and t:s in the current model solve and previously realized t:s"
......
......@@ -159,10 +159,10 @@ Parameters
// Aliases used in the equations after interval aggregation
// NOTE: Sample dimension has to be last because of the scenario reduction algorithm
ts_influx_(grid, node, f, t, s) "Mean external power inflow/outflow during a time step (MWh/h)"
ts_cf_(flow, node, f, t, s) "Mean available capacity factor time series (p.u.)"
ts_influx_(grid, node, s, f, t) "Mean external power inflow/outflow during a time step (MWh/h)"
ts_cf_(flow, node, s, f, t) "Mean available capacity factor time series (p.u.)"
ts_reserveDemand_(restype, up_down, node, f, t) "Mean reserve demand in region in the time step (MW)"
ts_node_(grid, node, param_gnBoundaryTypes, f, t, s) "Mean value of ts_node"
ts_node_(grid, node, param_gnBoundaryTypes, s, f, t) "Mean value of ts_node"
ts_fuelPrice_(fuel, t) "Mean fuel price time during time step (EUR/MWh)"
// Aliases used for updating data in inputsLoop.gms
......@@ -185,6 +185,9 @@ Parameters
// Bounds for scenario smoothening
p_tsMinValue(*, node, timeseries) "Minimum allowed value of timeseries for grid/flow and node"
p_tsMaxValue(*, node, timeseries) "Maximum allowed value of timeseries in grid/flow and node"
// Help parameters for scenario reduction
ts_energy_(s) "Total energy available from inflow and other flows (MWh)"
;
* --- Other time dependent parameters -----------------------------------------
......
......@@ -78,7 +78,7 @@ q_balance(gn(grid, node), msft(m, s, f, t)) // Energy/power balance dynamics sol
- v_spill(grid, node, s, f, t)${node_spill(node)}
// Power inflow and outflow timeseries to/from the node
+ ts_influx_(grid, node, f, t, s) // Incoming (positive) and outgoing (negative) absolute value time series
+ ts_influx_(grid, node, s, f, t) // Incoming (positive) and outgoing (negative) absolute value time series
// Dummy generation variables, for feasibility purposes
+ vq_gen('increase', grid, node, s, f, t) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
......@@ -293,7 +293,7 @@ q_maxDownward(gnu(grid, node, unit), msft(m, s, f, t))
* [
// Capacity factors for flow units
+ sum(flowUnit(flow, unit),
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
) // END sum(flow)
+ 1${not unit_flow(unit)}
] // END * p_unit(availability)
......@@ -359,7 +359,7 @@ q_maxDownwardOfflineReserve(gnu(grid, node, unit), msft(m, s, f, t))
* [
// Capacity factors for flow units
+ sum(flowUnit(flow, unit),
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
) // END sum(flow)
+ 1${not unit_flow(unit)}
] // END * p_unit(availability)
......@@ -438,7 +438,7 @@ q_maxUpward(gnu(grid, node, unit), msft(m, s, f, t))
* [
// Capacity factor for flow units
+ sum(flowUnit(flow, unit),
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
) // END sum(flow)
+ 1${not unit_flow(unit)}
] // END * p_unit(availability)
......@@ -524,7 +524,7 @@ q_maxUpwardOfflineReserve(gnu(grid, node, unit), msft(m, s, f, t))
* [
// Capacity factor for flow units
+ sum(flowUnit(flow, unit),
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
) // END sum(flow)
+ 1${not unit_flow(unit)}
] // END * p_unit(availability)
......@@ -572,7 +572,7 @@ q_reserveProvision(nuRescapable(restypeDirectionNode(restype, up_down, node), un
* [
// ... and capacity factor for flow units
+ sum(flowUnit(flow, unit),
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
) // END sum(flow)
+ 1${not unit_flow(unit)}
] // How to consider reserveReliability in the case of investments when we typically only have "realized" time steps?
......@@ -604,7 +604,7 @@ q_reserveProvisionOnline(nuRescapable(restypeDirectionNode(restype, up_down, nod
* [
// ... and capacity factor for flow units
+ sum(flowUnit(flow, unit),
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
) // END sum(flow)
+ 1${not unit_flow(unit)}
] // How to consider reserveReliability in the case of investments when we typically only have "realized" time steps?
......@@ -2033,7 +2033,7 @@ q_stateSlack(gn_stateSlack(grid, node), slack, sft(s, f, t))
* [
+ v_state(grid, node, s, f, t)
- p_gnBoundaryPropertiesForStates(grid, node, slack, 'constant')$p_gnBoundaryPropertiesForStates(grid, node, slack, 'useConstant')
- ts_node_(grid, node, slack, f, t, s)${ p_gnBoundaryPropertiesForStates(grid, node, slack, 'useTimeSeries') }
- ts_node_(grid, node, slack, s, f, t)${ p_gnBoundaryPropertiesForStates(grid, node, slack, 'useTimeSeries') }
] // END * p_slackDirection
;
......@@ -2049,7 +2049,7 @@ q_stateUpwardLimit(gn_state(grid, node), msft(m, s, f, t))
+ [
// Upper boundary of the variable
+ p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'constant')${p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'useConstant')}
+ ts_node_(grid, node, 'upwardLimit', f, t, s)${ p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'useTimeseries') }
+ ts_node_(grid, node, 'upwardLimit', s, f, t)${ p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'useTimeseries') }
// Investments
+ sum(gnu(grid, node, unit),
......@@ -2126,7 +2126,7 @@ q_stateDownwardLimit(gn_state(grid, node), msft(m, s, f, t))
// Lower boundary of the variable
- p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'constant')${p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'useConstant')}
- ts_node_(grid, node, 'downwardLimit', f, t, s)${ p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'useTimeseries') }
- ts_node_(grid, node, 'downwardLimit', s, f, t)${ p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'useTimeseries') }
] // END Headroom
* [
// Conversion to energy
......@@ -2372,7 +2372,7 @@ q_instantaneousShareMax(group, sft(s, f, t))
* [
// External power inflow/outflow
- sum(gnGroup(grid, node, group),
+ ts_influx_(grid, node, f, t, s)
+ ts_influx_(grid, node, s, f, t)
) // END sum(gnGroup)
// Consumption of units
......@@ -2469,7 +2469,7 @@ q_capacityMargin(gn(grid, node), sft(s, f, t))
},
// Capacity factors for flow units
+ sum(flowUnit(flow, unit)${ unit_flow(unit) },
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
) // END sum(flow)
* p_unit(unit, 'availability')
* [
......@@ -2518,7 +2518,7 @@ q_capacityMargin(gn(grid, node), sft(s, f, t))
) // END sum(gnu_input)
// Energy influx
+ ts_influx_(grid, node, f, t, s)
+ ts_influx_(grid, node, s, f, t)
// Capacity margin feasibility dummy variables
+ vq_capacity(grid, node, s, f, t)
......@@ -2630,7 +2630,7 @@ q_energyShareMax(group)
- p_groupPolicy(group, 'energyShareMax')
* [
- sum(gnGroup(grid, node, group),
+ ts_influx_(grid, node, f, t, s)
+ ts_influx_(grid, node, s, f, t)
) // END sum(gnGroup)
- sum(gnu_input(grid, node, unit)${ p_gnu(grid, node, unit, 'unitSizeCons')
and gnGroup(grid, node, group)
......@@ -2670,7 +2670,7 @@ q_energyShareMin(group)
- p_groupPolicy(group, 'energyShareMin')
* [
- sum(gnGroup(grid, node, group),
+ ts_influx_(grid, node, f, t, s)
+ ts_influx_(grid, node, s, f, t)
) // END sum(gnGroup)
- sum(gnu_input(grid, node, unit)${ p_gnu(grid, node, unit, 'unitSizeCons')
and gnGroup(grid, node, group)
......
......@@ -151,7 +151,7 @@ Option clear = p_stepLength;
Option clear = msft;
Option clear = mft;
Option clear = ft;
Option clear = sft, clear = fts;
Option clear = sft;
Option clear = mst_start, clear = mst_end;
$ifthen declared scenario
if(mSettings(mSolve, 'scenarios'), // Only clear these if using long-term scenarios
......@@ -333,7 +333,6 @@ 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;
fts(ft(f, t), s)$sft(s, f, t) = yes;
* Build stochastic tree by definfing previous samples
$ifthen defined scenario
......
......@@ -319,7 +319,7 @@ $ontext
)
/ mInterval(mSolve, 'stepsPerInterval', counter);
$offtext
ts_influx_(gn(grid, node), fts(f, tt_interval(t), s))
ts_influx_(gn(grid, node), sft(s, f, tt_interval(t)))
= sum(tt_aggregate(t, t_),
ts_influx(grid, node,
f + ( df_realization(f, t)$(not gn_forecasts(grid, node, 'ts_influx'))
......@@ -328,7 +328,7 @@ $offtext
+ dt_circular(t_)$(not gn_scenarios(grid, node, 'ts_influx'))))
)
/ mInterval(mSolve, 'stepsPerInterval', counter);
ts_cf_(flowNode(flow, node), fts(f, tt_interval(t), s))
ts_cf_(flowNode(flow, node), sft(s, f, tt_interval(t)))
= sum(tt_aggregate(t, t_),
ts_cf(flow, node,
f + ( df_realization(f, t)$(not gn_forecasts(flow, node, 'ts_cf'))
......@@ -347,7 +347,7 @@ $offtext
t_+ dt_circular(t_))
)
/ mInterval(mSolve, 'stepsPerInterval', counter);
ts_node_(gn_state(grid, node), param_gnBoundaryTypes, fts(f, tt_interval(t), s))
ts_node_(gn_state(grid, node), param_gnBoundaryTypes, sft(s, f, tt_interval(t)))
$p_gnBoundaryPropertiesForStates(grid, node, param_gnBoundaryTypes, 'useTimeseries')
// Take average if not a limit type
= (sum(tt_aggregate(t, t_),
......
......@@ -57,14 +57,14 @@ v_state.fx(gn_state(grid, node), sft(s, f,t))${ mft_lastSteps(mSolve, f, t)
v_state.up(gn_state(grid, node), sft(s, f, t))${ p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'useTimeSeries')
and not df_central(f,t)
}
= ts_node_(grid, node, 'upwardLimit', f, t, s)
= ts_node_(grid, node, 'upwardLimit', s, f, t)
* p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'multiplier')
;
// Lower bound
v_state.lo(gn_state(grid, node), sft(s, f, t))${ p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'useTimeSeries')
and not df_central(f,t)
}
= ts_node_(grid, node, 'downwardLimit', f, t, s)
= ts_node_(grid, node, 'downwardLimit', s, f, t)
* p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'multiplier')
;
// Fixed value
......@@ -72,7 +72,7 @@ v_state.fx(gn_state(grid, node), sft(s, f, t))${ p_gn(grid, node, 'boundAll')
and p_gnBoundaryPropertiesForStates(grid, node, 'reference', 'useTimeSeries')
and not df_central(f,t)
}
= ts_node_(grid, node, 'reference', f, t, s)
= ts_node_(grid, node, 'reference', s, f, t)
* p_gnBoundaryPropertiesForStates(grid, node, 'reference', 'multiplier')
;
// BoundEnd to a timeseries value
......@@ -80,7 +80,7 @@ v_state.fx(gn_state(grid, node), sft(s, f,t))${ mft_lastSteps(mSolve, f, t)
and p_gn(grid, node, 'boundEnd')
and p_gnBoundaryPropertiesForStates(grid, node, 'reference', 'useTimeSeries')
}
= ts_node_(grid, node, 'reference', f, t, s)
= ts_node_(grid, node, 'reference', s, f, t)
* p_gnBoundaryPropertiesForStates(grid, node, 'reference', 'multiplier');
// BoundStartToEnd: bound the last interval in the horizon to the value just before the horizon
......@@ -115,7 +115,7 @@ v_spill.lo(gn(grid, node_spill), sft(s, f, t))${ p_gnBoundaryPropertiesForSta
* p_gnBoundaryPropertiesForStates(grid, node_spill, 'minSpill', 'multiplier')
;
v_spill.lo(gn(grid, node_spill), sft(s, f, t))${ p_gnBoundaryPropertiesForStates(grid, node_spill, 'minSpill', 'useTimeSeries') }
= ts_node_(grid, node_spill, 'minSpill', f, t, s)
= ts_node_(grid, node_spill, 'minSpill', s, f, t)
* p_gnBoundaryPropertiesForStates(grid, node_spill, 'minSpill', 'multiplier')
;
v_spill.up(gn(grid, node_spill), sft(s, f, t))${ p_gnBoundaryPropertiesForStates(grid, node_spill, 'maxSpill', 'constant') }
......@@ -123,7 +123,7 @@ v_spill.up(gn(grid, node_spill), sft(s, f, t))${ p_gnBoundaryPropertiesForSta
* p_gnBoundaryPropertiesForStates(grid, node_spill, 'maxSpill', 'multiplier')
;
v_spill.up(gn(grid, node_spill), sft(s, f, t))${ p_gnBoundaryPropertiesForStates(grid, node_spill, 'maxSpill', 'useTimeSeries') }
= ts_node_(grid, node_spill, 'maxSpill', f, t, s)
= ts_node_(grid, node_spill, 'maxSpill', s, f, t)
* p_gnBoundaryPropertiesForStates(grid, node_spill, 'maxSpill', 'multiplier')
;
......@@ -143,7 +143,7 @@ v_gen.up(gnu_output(grid, node, unit_flow), sft(s, f, t))${gnuft(grid, node, uni
= sum(flow${ flowUnit(flow, unit_flow)
and nu(node, unit_flow)
},
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
* p_gnu(grid, node, unit_flow, 'maxGen')
* p_unit(unit_flow, 'availability')
) // END sum(flow)
......@@ -170,7 +170,7 @@ v_gen.lo(gnu_input(grid, node, unit_flow), sft(s, f, t))${gnuft(grid, node, unit
= - sum(flow${ flowUnit(flow, unit_flow)
and nu(node, unit_flow)
},
+ ts_cf_(flow, node, f, t, s)
+ ts_cf_(flow, node, s, f, t)
* p_gnu(grid, node, unit_flow, 'maxCons')
* p_unit(unit_flow, 'availability')
) // END sum(flow)
......
......@@ -186,8 +186,8 @@ d_capacityFactor(flowNode(flow, node), sft(s, f_solve(f), t_current(t)))
and t_active(t)
and sum(flowUnit(flow, unit), nu(node, unit))
}
= ts_cf_(flow, node, f, t, s)
+ ts_cf(flow, node, f, t + dt_scenarioOffset(flow, node, 'ts_cf', s))${ not ts_cf_(flow, node, f, t, s) }
= ts_cf_(flow, node, s, f, t)
+ ts_cf(flow, node, f, t + dt_scenarioOffset(flow, node, 'ts_cf', s))${ not ts_cf_(flow, node, s, f, t) }
+ Eps
;
// Temperature forecast for examining the error
......@@ -196,8 +196,8 @@ d_nodeState(gn_state(grid, node), param_gnBoundaryTypes, sft(s, f_solve(f), t_cu
and t_active(t)
and msf(mSolve, s, f)
}
= 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)}
= ts_node_(grid, node, param_gnBoundaryTypes, s, f, t)
+ ts_node(grid, node, param_gnBoundaryTypes, f, t)${ not ts_node_(grid, node, param_gnBoundaryTypes, s, f, t)}
+ Eps
;
// Influx forecast for examining the errors
......@@ -205,8 +205,8 @@ d_influx(gn(grid, node), sft(s, f_solve(f), t_current(t)))
${ 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)}
= ts_influx_(grid, node, s, f, t)
+ ts_influx(grid, node, f, t)${ not ts_influx_(grid, node, s, f, t)}
+ Eps
;
// Scenario values for time series
......@@ -216,8 +216,8 @@ loop(s_scenario(s, scenario),
d_state(gn_state(grid, node), scenario, f, t) = v_state.l(grid, node, s, f, t);
);
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;
d_ts_scenarios('ts_influx', gn, scenario, ft)$sft(s, ft) = ts_influx_(gn, s, ft) + eps;
d_ts_scenarios('ts_cf', flowNode, scenario, ft)$sft(s, ft) = ts_cf_(flowNode, s, ft) + eps;
);
$endif.diag
......
......@@ -19,10 +19,24 @@ ScenRedParms('visual_init') = 1;
ScenRedParms('visual_red') = 1;
$endif
* Calculate data for scenario reduction
Option clear = ts_energy_;
ts_energy_(s)
= sum(sft(s, f, t)$mf_central(mSolve, f),
sum(gn(grid, node),
ts_influx_(grid, node, s, f, t)
+ sum((flowNode(flow, node), flowUnit(flow, unit))
$gnu_output(grid, node, unit),
ts_cf_(flow, node, s, f, t)
* p_gnu(grid, node, unit, 'maxGen')
)
) * p_stepLength(mSolve, f, t)
);
* Export data
execute_unload 'srin.gdx', ScenRedParms,
s_active, ss, p_sProbability,
ts_influx_, ts_cf_;
ts_energy_;
* Choose right null device
$ifthen %system.filesys% == 'MSNT' $set nuldev NUL
$else $set nuldev /dev/null
......@@ -34,7 +48,7 @@ execute 'scenred2 inc/scenred.cmd > %nuldev%';
if(errorLevel,
put log "!!! Scenario reduction (SCENRED2) failed. ";
put log "See file 'sr.log' for details."/;
put_utility log, 'Click' / 'sr.log';
put_utility log, 'Click' / 'sr.log';
putclose;
execError = execError + 1;
else
......@@ -78,11 +92,10 @@ else
msf(mSolve, s, f)$msf(mSolve, s, f) = ms(mSolve, s);
msft(mSolve, s, f, t)$msft(mSolve, s, f, t) = msf(mSolve, s, f);
sft(s, f, t)$sft(s, f, t) = yes$p_msProbability(mSolve, s);
fts(f, t, s)$fts(f, t, s) = sft(s, f, t);
// Clear data from removed samples
ts_influx_(gn, ft, s_active(s))$(not p_sProbability(s)) = 0;
ts_cf_(flowNode, ft, s_active(s))$(not p_sProbability(s)) = 0;
ts_influx_(gn, s_active(s), ft)$(not p_sProbability(s)) = 0;
ts_cf_(flowNode, s_active(s), ft)$(not p_sProbability(s)) = 0;
s_active(s)$s_active(s) = yes$p_sProbability(s);
);
......
......@@ -28,12 +28,12 @@ $else
$setlocal linking_set gn
$endif
%timeseries%_(%linking_set%, fts(f, t, s))$[p_autocorrelation(%linking_set%,
%timeseries%_(%linking_set%, sft(s, f, t))$[p_autocorrelation(%linking_set%,
'%timeseries%')
and ms_central(mSolve, s)]
= min(p_tsMaxValue(%linking_set%, '%timeseries%'),
max(p_tsMinValue(%linking_set%, '%timeseries%'),
%timeseries%_(%linking_set%, f, t, s)
%timeseries%_(%linking_set%, s, f, t)
+ (%timeseries%(%linking_set%,
f_ + (df_realization(f_, t_)
$(not gn_forecasts(%linking_set%, '%timeseries%'))),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment