Commit 636c1d80 authored by Juha Kiviluoma's avatar Juha Kiviluoma
Browse files

Bug fixes for v_state and ft_dynamic related issues

parent 48f455d6
fSolve(f) = no;
fSolve(f)$mf(mSolve,f) = yes;
tDispatchCurrent = 0; // Reset dispatch loop
tSolveFirst = ord(tSolve); // tSolveFirst: the start of the current solve
tSolveLast = ord(tSolve) + max(mSettings(mSolve, 't_forecastLength'), mSettings(mSolve, 't_horizon')); // tSolveLast: the end of the current solve
p_stepLength(mSolve, f, t) = no;
......@@ -39,7 +40,7 @@ $offOrder
ft_new(f,t_)$(mf(mSolve, f) and tInterval(t_)) = yes;
p_stepLengthNoReset(mf(mSolve, fSolve), t) = intervalLength * mSettings(mSolve, 'IntervalInHours');
// Aggregates the interval time series data by averaging the power data
ts_influx_(grid, node, fSolve, t) = sum{t_$tInterval(t_), ts_influx(grid, node, fSolve, t_+ct(t_))} / p_stepLength(mSolve, fSolve, t); // Averages the absolute power terms over the interval
ts_influx_(grid, node, fSolve, t) = sum{t_$tInterval(t_), ts_influx(grid, node, 'f00', t_+ct(t_))} / p_stepLength(mSolve, fSolve, t); // Averages the absolute power terms over the interval
ts_cf_(flow, node, fSolve, t) = sum{t_$tInterval(t_), ts_cf(flow, node, fSolve, t_+ct(t_))} / p_stepLength(mSolve, fSolve, t); // Averages the capacity factor over the inverval
ts_nodeState_(gn_state(grid, node), param_gnBoundaryTypes, fSolve, t) = sum(t_${tInterval(t_)}, ts_nodeState(grid, node, param_gnBoundaryTypes, fSolve, t_+ct(t_))) / p_stepLength(mSolve, fSolve, t); // Averages the time-dependent node state boundary conditions over the interval
ts_unit_(unit, param_unit, fSolve, t) = sum(t_${tInterval(t_)}, ts_unit(unit, param_unit, fSolve, t_+ct(t_))) / p_steplength(mSolve, fSolve, t); // Averages the time-dependent unit parameters over the interval
......
......@@ -58,15 +58,20 @@
);
);
ft_dynamic(f,tSolve+max(0,(tDispatchCurrent -1))) = no;
ft_dynamic('f00',tSolve+tDispatchCurrent)$(tDispatchCurrent > 0) = yes;
//ft_dynamic('f00',tSolve+(tDispatchCurrent-1))$(tDispatchCurrent > 1 or (tSolveFirst > 1 and tDispatchCurrent > 0)) = yes;
ft_dynamic('f00',tSolve+(tDispatchCurrent))$(tDispatchCurrent > 0) = yes;
//ft_dynamic(f,tSolve)$(tSolveFirst > 0 and tDispatchCurrent = 0 and ord(f) > 1) = 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_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;
ft_realizedLast(f,t)$[fRealization(f) and ord(t) = ord(tSolve) + mSettings(mSolve, 't_jump') - 1] = yes;
ft_fix(f,t) = no;
ft_fix(f,tSolve + tDispatchCurrent) = yes;
// Defining unit aggregations and ramps
uft(unit, f, t) = no;
......@@ -110,11 +115,15 @@
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;
p_uft_online_last(unit, f, t) = ord(t);
);
pf(f,t) = no;
pf(ft(f,t))$(ord(t) eq ord(tSolve) + tDispatchCurrent) = 1 - ord(f);
p_sft_probability(s, f, t)$(sInitial(s) and ft(f,t)) = p_fProbability(f) / sum(f_$ft(f_,t), p_fProbability(f_));
p_sft_probability(s, f, t)$(sInitial(s) and fCentral(f) and ord(t) = tSolveFirst + mSettings(mSolve, 't_horizon')) = p_fProbability(f) / sum(f_$(fCentral(f_) and ord(t) = tSolveFirst + mSettings(mSolve, 't_horizon')), p_fProbability(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) ,
......
* --- Model parameters, features and switches ---------------------------------
Sets // Model related selections
mType "model types in the Backbone"
/invest, storage, schedule, realtime, building/
/invest, storage, schedule, realtime, building, schedule_dispatch/
mSetting "setting categories for models"
/t_start, t_jump, t_horizon, t_forecastStart, t_forecastLength, t_forecastJump, t_end, t_reserveLength, samples, forecasts, readForecastsInTheLoop, intervalEnd, intervalLength, IntervalInHours, t_aggregate/
counter "general counter set"
/c000*c999/
solveInfoAttributes
/
modelStat
solveStat
totalTime
iterations
nodes
numEqu
numDVar
numVar
numNZ
sumInfes
objEst
objVal
/
// Efficiency approximation related sets
op "Operating points in the efficiency curves, also functions as index for data points"
/op00*op12/ // IMPORTANT! Has to equal the same param_unit!
......
......@@ -64,6 +64,7 @@ Sets
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_fix(f,t) "ft where variable values are fixed after the solve, so that dynamic equations cannot alter them"
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"
......
......@@ -47,6 +47,7 @@ Parameters
p_sWeight(s) "Weight of sample"
p_sProbability(s) "Probability to reach sample conditioned on anchestor samples"
p_fProbability(f) "Probability of forecast"
p_sft_probability(s, f, t) "Probability of forecast"
;
Scalar p_sWeightSum "Sum of sample weights";
......@@ -63,6 +64,8 @@ Parameters
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)"
aaSolveInfo(mType, t, solveInfoAttributes) "stores information about the solve status"
p_uft_online_last(unit, f, t) "Ord of last t where unit is online"
;
* --- Stochastic data parameters ----------------------------------------------
......
......@@ -63,6 +63,7 @@ unit_minload(unit)$[p_unit(unit, 'op00') > 0 and p_unit(unit, 'op00') < 1] = yes
unit_flow(unit)$sum(flow, flowUnit(flow, unit)) = yes;
unit_fuel(unit)$sum[ (fuel, node)$sum(t, ts_fuelPriceChangenode(fuel, node, t)), uFuel(unit, 'main', fuel) ] = yes;
unit_elec(unit)$sum(gnu(grid, node, unit), p_gnu('elec', node, unit, 'maxGen')) = yes;
unit_elec(unit)$sum(gnu(grid, node, unit), p_gnu('elec', node, unit, 'maxCons')) = yes;
* Assume values for critical unit related parameters, if not provided by input data
p_unit(unit, 'eff00')$(not p_unit(unit, 'eff00')) = 1; // If the unit does not have efficiency set, it is 1
......
......@@ -40,8 +40,7 @@ q_obj ..
+ v_obj * 1000000
=E=
+ sum(msft(m, s, f, t),
p_sProbability(s) *
p_fProbability(f) *
p_sft_Probability(s,f,t) *
(
// Variable O&M costs
+ sum(gnu_output(grid, node, unit), // Calculated only for output energy
......@@ -66,8 +65,6 @@ q_obj ..
+ sum(uft_online(unit, f, t),
+ {
+ v_startup(unit, f, t) // Cost of starting up
- sum(t_${mftStart(m, f, t_) and uft_online(unit, f, t_)}, 0.5 * v_online(unit, f, t_)) // minus value of avoiding startup costs before
- sum(t_${mftLastSteps(m, f, t_) and uft_online(unit, f, t_)}, 0.5 * v_online(unit, f, t_)) // or after the model solve
} / p_unit(unit, 'unitCount')
* {
// Startup variable costs
......@@ -99,18 +96,26 @@ q_obj ..
+ p_gnu(grid, node, unit, 'rampDownCost') * v_genRampChange(grid, node, unit, 'down', f, t)
)
)
) // END * p_sProbability(s) & p_fProbability(f)
) // END * p_sft_probability(s,f,t)
) // END sum over msft(m, s, f, t)
// Value of energy storage change
- sum((mftLastSteps(m, f, t), mftStart(m, f_, t_)) $(active('storageValue')),
p_fProbability(f) *
sum(gn_state(grid, node),
- sum(mftLastSteps(m, f, t)$active('storageValue'),
sum(gn_state(grid, node),
p_storageValue(grid, node, t) *
(v_state(grid, node, f, t) - v_state(grid, node, f_, t_))
)
sum(s$p_sft_probability(s,f,t), p_sft_probability(s, f, t) * v_state(grid, node, f, t))
)
)
+ sum(mftStart(m, f, t)$active('storageValue'),
sum(gn_state(grid, node),
p_storageValue(grid, node, t) *
sum(s$p_sft_probability(s,f,t), p_sft_probability(s, f, t) * v_state(grid, node, f, t))
)
)
- sum([s, m, uft_online(unit, f, t)]$mftStart(m, f, t), p_sft_probability(s, f, t) * 0.5 * v_online(unit, f, t)) // minus value of avoiding startup costs before
- sum((s, uft_online(unit, f, t))$(ord(t) = p_uft_online_last(unit, f, t)), p_sft_probability(s, f, t) * 0.5 * v_online(unit, f, t)) // or after the model solve
// Dummy variables
+ sum(msft(m, s, f, t), p_sProbability(s) * p_fProbability(f) * (
+ sum(msft(m, s, f, t), p_sft_probability(s, f, t) * (
sum(inc_dec,
sum( gn(grid, node), vq_gen(inc_dec, grid, node, f, t) * (p_stepLength(m, f, t) + p_stepLength(m, f+pf(f,t), t+pt(t))${not p_stepLength(m, f, t)}) * PENALTY_BALANCE(grid) )
)
......@@ -125,17 +130,18 @@ q_obj ..
+ sum(gn_stateSlack(grid, node),
+ sum(msft(m, s, f, t),
+ sum(upwardSlack$p_gnBoundaryPropertiesForStates(grid, node, upwardSlack, 'slackCost'),
+ p_sProbability(s) * p_fProbability(f) * ( p_gnBoundaryPropertiesForStates(grid, node, upwardSlack, 'slackCost') * v_stateSlack(grid, node, upwardSlack, f, t) )
+ p_sft_probability(s, f, t) * ( p_gnBoundaryPropertiesForStates(grid, node, upwardSlack, 'slackCost') * v_stateSlack(grid, node, upwardSlack, f, t) )
)
+ sum(downwardSlack$p_gnBoundaryPropertiesForStates(grid, node, downwardSlack, 'slackCost'),
+ p_sProbability(s) * p_fProbability(f) * ( p_gnBoundaryPropertiesForStates(grid, node, downwardSlack, 'slackCost') * v_stateSlack(grid, node, downwardSlack, f, t) )
+ p_sft_probability(s, f, t) * ( p_gnBoundaryPropertiesForStates(grid, node, downwardSlack, 'slackCost') * v_stateSlack(grid, node, downwardSlack, f, t) )
)
)
)
;
* -----------------------------------------------------------------------------
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
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+cf(f,t), t) // The difference between current
......@@ -187,7 +193,6 @@ q_balance(gn(grid, node), m, ft_dynamic(f, t))${p_stepLength(m, f+pf(f,t), t+pt(
// 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+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+cf(f,t), t)
......@@ -198,7 +203,7 @@ q_balance(gn(grid, node), m, ft_dynamic(f, t))${p_stepLength(m, f+pf(f,t), t+pt(
$$ifi '%rampSched%' == 'yes' / 2 // Averaging all the terms on the right side of the equation over the timestep here.
;
* -----------------------------------------------------------------------------
q_resDemand(restypeDirectionNode(restype, up_down, node), ft(f, t))$(ord(t) le tSolveFirst + sum(mf(m, f), mSettings(m, 't_reserveLength'))) ..
q_resDemand(restypeDirectionNode(restype, up_down, node), ft(f, t))${ord(t) le tSolveFirst + sum[mf(m, f), mSettings(m, 't_reserveLength')] and [not (restype('tertiary') and ord(t) le tSolveFirst + tDispatchCurrent)]} ..
+ sum(nuft(node, unit, f, t)$nuRescapable(restype, up_down, node, unit), // Reserve capable units on this node
v_reserve(restype, up_down, node, unit, f, t)
)
......@@ -295,7 +300,7 @@ q_startup(uft_online(unit, ft_dynamic(f, t))) ..
;
* -----------------------------------------------------------------------------
q_genRamp(gn(grid, node), m, unit, ft_dynamic(f, t))${gnuft_ramp(grid, node, unit, f, t)} ..
+ v_genRamp(grid, node, unit, f+pf(f,t), t+pt(t))
+ v_genRamp(grid, node, unit, f, t+pt(t))
* ( p_gnu(grid, node, unit, 'maxGen') + p_gnu(grid, node, unit, 'maxCons') ) * 60 / 100 // Unit conversion from [p.u./min] to [MW/h]
* p_stepLength(m, f+pf(f,t), t+pt(t))
=E=
......@@ -313,11 +318,11 @@ q_genRamp(gn(grid, node), m, unit, ft_dynamic(f, t))${gnuft_ramp(grid, node, uni
;
* -----------------------------------------------------------------------------
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+cf(f,t), t)
- v_genRampChange(grid, node, unit, 'down', f+cf(f,t), t)
+ v_genRampChange(grid, node, unit, 'up', f, t+pt(t))
- v_genRampChange(grid, node, unit, 'down', f, t+pt(t))
=E=
+ v_genRamp(grid, node, unit, f+cf(f,t), t)
- v_genRamp(grid, node, unit, f+pf(f,t), t+pt(t))
+ v_genRamp(grid, node, unit, f, t+pt(t))
- v_genRamp(grid, node, unit, f, t+pt(t))
;
* -----------------------------------------------------------------------------
q_conversionDirectInputOutput(suft(effDirect, unit, f, t)) ..
......
* --- Variable limits ---------------------------------------------------------
// v_state absolute boundaries set according to p_gn parameters;
// When using constant values and to supplement time series with constant values (time series will override when data available)
v_state.up(gn_state(grid, node), ft_limits(f, t))$p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'useConstant') = p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'constant') * p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'multiplier');
v_state.up(gn_state(grid, node), ft_limits(f, t))${p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'useConstant') and not ft_fix(f,t)} = p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'constant') * p_gnBoundaryPropertiesForStates(grid, node, 'upwardLimit', 'multiplier');
v_state.lo(gn_state(grid, node), ft_limits(f, t))$p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'useConstant') = p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'constant') * p_gnBoundaryPropertiesForStates(grid, node, 'downwardLimit', 'multiplier');
v_state.fx(gn_state(grid, node), ft_limits(f, t))$(p_gn(grid, node, 'boundAll') and p_gnBoundaryPropertiesForStates(grid, node, 'reference', 'useConstant')) = p_gnBoundaryPropertiesForStates(grid, node, 'reference', 'constant') * p_gnBoundaryPropertiesForStates(grid, node, 'reference', 'multiplier');
// When using time series
......
if (mSolve('schedule'),
settings(mSetting) = mSettings(mSolve, mSetting);
solve schedule using mip minimizing v_obj;
tDispatchCurrent = 0;
aaSolveInfo('schedule', tSolve, 'modelStat') = schedule.modelStat;
aaSolveInfo('schedule', tSolve, 'solveStat') = schedule.solveStat;
aaSolveInfo('schedule', tSolve, 'totalTime') = schedule.etSolve;
aaSolveInfo('schedule', tSolve, 'iterations') = schedule.iterUsd;
aaSolveInfo('schedule', tSolve, 'nodes') = schedule.nodUsd;
aaSolveInfo('schedule', tSolve, 'numEqu') = schedule.numEqu;
aaSolveInfo('schedule', tSolve, 'numDVar') = schedule.numDVar;
aaSolveInfo('schedule', tSolve, 'numVar') = schedule.numVar;
aaSolveInfo('schedule', tSolve, 'numNZ') = schedule.numNZ;
aaSolveInfo('schedule', tSolve, 'sumInfes') = schedule.sumInfes;
aaSolveInfo('schedule', tSolve, 'objEst') = schedule.objEst;
aaSolveInfo('schedule', tSolve, 'objVal') = schedule.objVal;
put_utility 'gdxout' / 'output\debug_', tSolveFirst:0:0, '_', tDispatchCurrent:0:0, '.gdx';
$$ifi '%debug%' == 'yes' execute_unload; // Output debugging information
loop(restypeDirectionNode(restypeDirection(restype, up_down), node),
......@@ -9,6 +22,7 @@
v_reserve.fx(nuRescapable(restype, up_down, node, unit_elec), ft(f,t))$( nuft(node, unit_elec, f, t)
and ord(t) > p_nReserves(node, restype, 'gate_closure')
and ord(t) <= p_nReserves(node, restype, 'gate_closure') + p_nReserves(node, restype, 'gate_closure')
and not unit_flow(unit_elec) // NOTE! Units using flows can change their reserve (they might not have as much available in real time as they had bid)
)
= v_reserve.l(restype, up_down, node, unit_elec, f, t);
);
......@@ -25,15 +39,30 @@
);
solve schedule_dispatch using mip minimizing v_obj;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'modelStat') = schedule_dispatch.modelStat;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'solveStat') = schedule_dispatch.solveStat;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'totalTime') = schedule_dispatch.etSolve;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'iterations') = schedule_dispatch.iterUsd;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'nodes') = schedule_dispatch.nodUsd;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'numEqu') = schedule_dispatch.numEqu;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'numDVar') = schedule_dispatch.numDVar;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'numVar') = schedule_dispatch.numVar;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'numNZ') = schedule_dispatch.numNZ;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'sumInfes') = schedule_dispatch.sumInfes;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'objEst') = schedule_dispatch.objEst;
aaSolveInfo('schedule_dispatch', tSolveDispatch, 'objVal') = schedule_dispatch.objVal;
put_utility 'gdxout' / 'output\debug_', tSolveFirst:0:0, '_d', tDispatchCurrent:0:0, '.gdx';
$$ifi '%debug%' == 'yes' execute_unload; // Output debugging information
// Remaining solves will use bound start value for v_state once it has been established
v_state.fx(grid, node, ft_limits(f, t))$(gn_state(grid,node) and tSolveDispatch(t))
= v_state.l(grid, node, f, t);
v_state.fx(grid, node, ft_fix(f,t))$gn_state(grid,node) = v_state.l(grid, node, f, t);
// Remaining solves
v_online.fx(uft_limits_online(unit, f, t))$(ord(t) = tDispatchCurrent + 1) = round(v_online.l(unit, f, t));
v_online.fx(uft_limits_online(unit, ft_fix(f, t))) = round(v_online.l(unit, f, t));
//v_gen.fx(gnu(grid, node, unit), ft_fix(f, t)) = v_gen.l(grid, node, unit, f, t);
//v_spill.fx(gn(grid, node), ft_fix(f, t)) = v_spill.l(grid, node, f, t);
//vq_gen.fx(inc_dec, gn(grid, node), ft_fix(f, t)) = vq_gen.l(inc_dec, grid, node, f, t);
);
);
* if (mSolve('storage'), solve storage using lp minimizing v_obj; );
Supports Markdown
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