Commit 8f4ac29e authored by unknown's avatar unknown
Browse files

First commit as Backbone. A working version for scheduling that probably contains several bugs.

parents
Backbone.gpr
Backbone.lst
Backbone.lxi
Backbone.txt
equations.lst
equations.lxi
*.~gm
input/
output/
225a/
$title Backbone
$ontext
Backbone - chronological energy systems model
===================================
Created by:
Juha Kiviluoma <juha.kiviluoma@vtt.fi>
Erkka Rinne <erkka.rinne@vtt.fi>
Based on Stochastic Model Predictive Control method [1]. Improved by
generalising the idea of storages and adding the ability to load storages.
Can handle multiple stochastic parameters.
GAMS command line arguments
---------------------------
--debug=[yes|no]
Switch on/off debugging mode. In debug mode, writes debug.gdx
with all symbols as well as a gdx file for each solution containing
model parameters, variables and equations.
--dummy=[yes|no]
Do not solve the model, just do preliminary calculations.
For testing purposes.
--<name of model parameter>=<value>
Set model parameter value. See file settings.inc for available
parameters.
--<name of model feature>=[yes|no]
Switch model features on/off. See file settings.inc for available
features.
References
----------
[1] K. Nolde, M. Uhr, and M. Morari, Medium term scheduling of a hydro-thermal
system using stochastic model predictive control, Automatica, vol. 44,
no. 6, pp. 15851594, Jun. 2008.
$offtext
* === Settings ================================================================
* Activate end of line comments and set comment character to '//'
$oneolcom
$eolcom //
* Set log file (output to IDE process window)
file log /''/;
* Allow empty data definitions
$onempty
option profile = 3;
* Load model settings
$include 'settings.inc'
* === Sets ====================================================================
$include 'inc/sets.gms'
* === Parameters ==============================================================
$include 'inc/parameters.gms'
$include 'inc/results.gms'
Scalars
errorcount /0/
elapsed "Model time elapsed since simulation start (t)" /0/
t_solveOrd "ord of t_solve"
;
* Debug arrays
Parameters
x_stoContent(storage, f, t) "Storage content at the end of the time period in a sample (ratio of max)"
x_storageControl(storage, f, t) "Storage control value during a time period in a sample (MWh)"
;
* === Macros ==================================================================
$include 'inc/macros.gms'
options
solvelink = %Solvelink.Loadlibrary%
$ifi not '%debug%' == 'yes'
solprint = Silent
;
* Load updates made for BackBone
$gdxin 'input/inputData.gdx'
$loadm param
$loaddc geo
$loaddc flow
$loaddc bus
$loaddc unit
$loaddc fuel
$loaddc storage
*$loaddc eg
*$loaddc unitVG
$loaddc hydroBus
*$loaddc gu
$loaddc egu
$loaddc egu_input
$loaddc ggu
$loaddc egs
$loaddc flow_unit
$loaddc unit_fuel
$loaddc unit_storage
$loaddc gu_fixed_output_ratio
$loaddc gu_constrained_output_ratio
*$loaddc resCapable
$loaddc emission
$loaddc ts_energyDemand
$loaddc ts_import
*$load ts_reserveDemand
$loaddc ts_cf
$loaddc ts_stoContent
$loaddc ts_fuelPriceChange
$loaddc ts_inflow
$loaddc p_data2d
$loaddc uData
$loaddc usData
$loaddc uReserveData
$loaddc p_transferCap
$loaddc p_transferLoss
*$loaddc etype_storage
$gdxin
* Generate sets based on parameter data
eg(etype, geo)$sum(unit, egu(etype, geo, unit)) = yes;
gu(geo, unit)$sum(etype, egu(etype, geo, unit)) = yes;
*egu(etype,geo,unit)$(gu(geo,unit) and eg(etype,geo)) = YES;
ggu(geo, geo_, unit)$(gu(geo, unit) and ord(geo) = ord(geo_)) = yes;
eg2g(etype, from_geo, to_geo)$p_transferCap(etype, from_geo, to_geo) = yes;
bus_to_bus(from_geo, to_geo)$p_transferCap('elec', from_geo, to_geo) = yes;
ts_fuelPriceChangeGeo(fuel, geo, t) = ts_fuelPriceChange(fuel, t);
unitOnline(unit)$[ sum(egu(etype, geo, unit), udata(etype, geo, unit, 'startup_cost') or udata(etype, geo, unit, 'startup_fuelcons') or udata(etype, geo, unit, 'leadtime') ) ] = yes;
unitVG(unit)$sum(flow, flow_unit(flow, unit)) = yes;
*unitConversion(unit)$sum(eg(etype, geo), egu_input(etype, geo, unit)) = yes;
unitElec(unit)$sum(egu(etype, geo, unit), udata('elec', geo, unit, 'max_cap')) = yes;
unitHeat(unit)$sum(egu(etype, geo, unit), udata('heat', geo, unit, 'max_cap')) = yes;
unitFuel(unit)$sum[ (fuel, geo)$sum(t, ts_fuelPriceChangeGeo(fuel, geo, t)), unit_fuel(unit, fuel, 'main') ] = yes;
unitVG(unit)$sum(flow, flow_unit(flow, unit)) = yes;
unitMinload(unit)$sum(egu(etype, geo, unit), p_data2d(etype, unit, 'min_load')) = yes;
unitHydro(unit)$sum(unit_fuel(unit,'WATER','main'), 1) = yes;
unitHydro(unit)$sum(unit_fuel(unit,'WATER_RES','main'), 1) = yes;
storageHydro(storage)$sum(unitHydro, unit_storage(unitHydro, storage)) = yes;
resCapable(resType, resDirection, geo, unit)$uReserveData(geo, unit, resType, resDirection) = yes;
* === Variables ===============================================================
$include 'inc/variables.gms'
* === equations ===============================================================
$include 'inc/equations.gms'
* Load stochastic scenarios
$batinclude 'inc/gdxload_fluctuation.inc' wind
$batinclude 'inc/gdxload_fluctuation.inc' solar
$ifthen exist 'input/scenarios_hydro.gdx'
$$gdxin 'input/scenarios_hydro.gdx'
$endif
$gdxin
* --- Data corrections etc. ---------------------------------------------------
$if exist 'extra_data.gms' $include 'extra_data.gms'
$include 'inc/schedule.gms'
* === Calculations ============================================================
* Define long-term storages
loop(egs(etype, geo, storage) $(usData(etype, geo, storage, 'max_content') > 0
and sum(unit_storage(unit, storage),
uData('elec', geo, unit, 'max_cap')
+ uData('heat', geo, unit, 'max_cap')) > 0),
storageLong(storage) = yes$(usData(etype, geo, storage, 'max_content')
/ sum(unit_storage(unit, storage),
uData('elec', geo, unit, 'max_cap')
+ uData('heat', geo, unit, 'max_cap'))
> 0.5 * modelSolveRules('schedule', 't_horizon') );
);
* Link units to genTypes
$iftheni '%genTypes%' == 'yes'
loop(gu_fuel(geo, unit, fuel, 'main'),
genType_g('pumped storage', unit) = yes$(sameas(fuel, 'water')
and p_data(unit, 'max_loading') > 0);
genType_g('hydropower', unit) = yes$(sameas(fuel, 'water')
and not genType_g('pumped storage', unit));
genType_g('nuclear', unit) = yes$sameas(fuel, 'nuclear');
genType_g('coal', unit) = yes$sameas(fuel, 'coal');
genType_g('OCGT', unit) = yes$(sameas(g, 'OCGT') or sameas(unit, 'DoE_Peaker'));
genType_g('CCGT', unit) = yes$(sameas(fuel, 'nat_gas')
and not genType_g('OCGT', unit));
genType_g('solar', unit) = yes$sameas(fuel, 'solar');
genType_g('wind', unit) = yes$sameas(fuel, 'wind');
genType_g('imports', unit) = yes$(not sameas(area, 'SA'));
genType_g('dummy', unit) = yes$sameas(unit, 'dummy');
);
$endif
* Calculate average hourly loads from demand
ts_energyDemand(etype, geo, f, t)$ts_energyDemand(etype, geo, f, t)
= ts_energyDemand(etype, geo, f, t) * p_data2d(etype, geo, 'annualDemand') / 1;
* Calculate power based time series for ramp scheduling
*if(active('rampSched'),
* $$include 'inc/rampSched/rampSchedTimeSeries_rampSearch.gms'
*);
*$include 'inc/rampSched/killStuff.gms'
* --- Use input data to select which samples are included in the model run
*$include 'input/samples_in_model.inc';
* --- Calculate trees ---------------------------------------------------------
*$include 'inc/treeCalc.inc'
* --- Extra calculations ------------------------------------------------------
$if exist 'extra_calculations.gms' $include 'extra_calculations.gms'
* === Files ===================================================================
files gdx, cmd;
file f_info /'output/info.txt'/;
* === Simulation ==============================================================
* --- Generate model rules ----------------------------------------------------
* Check the modelSolves for preset patterns for model solve timings
* If not found, then use modelSolveRules to set the model solve timings
loop(modelType,
if(sum(t$modelSolves(modelType, t), 1) = 0,
t_skip_counter = 0;
loop(t$( ord(t) = modelSolveRules(modelType, 't_start') + modelSolveRules(modelType, 't_jump') * t_skip_counter and ord(t) <= modelSolveRules(modelType, 't_end') ),
modelSolves(modelType, t)=yes;
t_skip_counter = t_skip_counter + 1;
);
);
);
* Select models to be used
m(modelType)=no;
loop(modelType,
// If there are solves for a model, then the model is turned on
if(sum(t$modelSolves(modelType, t), 1) > 0,
m(modelType) = yes;
);
);
* Select samples for the model
loop(m,
// Set scenarios in use for the models, if they haven't been provided in the data
if (not sum(s, ms(m, s)),
ms(m, s)$(ord(s) <= modelSolveRules(m, 'samples')) = yes;
// Use all scenarios if modelSolveRules/scenarios is 0
if (modelSolveRules(m, 'samples') = 0,
ms(m, s) = yes;
);
);
);
// If the model does not have preset step lengths...
loop(m,
if(sum[(ms(m,s), f, t)$fRealization(f), p_stepLength(m, f, t)] = 0,
// ...and if there is a t_interval for the model, then use it to set constant step lengths
if(modelSolveRules(m, 't_interval'),
p_stepLength(mf(m, f), t) = modelSolveRules(m, 't_interval');
);
);
);
loop(modelSolves(m_solve, t_solve),
t_solveOrd = ord(t_solve);
elapsed = t_solveOrd - modelSolveRules(m_solve, 't_start');
// Set mft for the modelling period and model forecasts
mft(m_solve,f,t) = no;
mft(m_solve, f, t)${ [ord(t) >= ord(t_solve)]
and [ord(t) <= ord(t_solve) + modelSolveRules(m_solve, 't_forecastLength')]
and mf(m_solve, f)
} = yes;
mftStart(m_solve,f,t) = no;
mftStart(m_solve,fRealization,t)$[ord(t) = ord(t_solve)] = yes;
mftLastForecast(m_solve,f,t) = no;
mftLastForecast(m_solve,f,t)$[ord(f)-1 <= modelSolveRules(m_solve, 'forecasts') and ord(t) = ord(t_solve) + modelSolveRules(m_solve, 't_forecastLength')] = yes;
mftLastSteps(m_solve,f,t) = no;
mftLastSteps(m_solve,f,t)$[ord(f)-1 <= modelSolveRules(m_solve, 'forecasts') and ord(t) = ord(t_solve) + max(modelSolveRules(m_solve, 't_forecastLength'), modelSolveRules(m_solve, 't_horizon'))] = yes;
mftBind(m_solve,f,t) = no;
mft_bind(m_solve,f,t) = no;
mt_bind(m_solve,t) = no;
* mftBind(mft(m_solve,f,t))$[ord(t) = ord(t_solve) + modelSolveRules(m_solve, 't_forecastLength')] = yes;
* mft_bind(mft(m_solve,f,t))$[ord(t) = ord(t_solve) + modelSolveRules(m_solve, 't_forecastLength')] = 1 - ord(f);
* mt_bind(m_solve,t)$[ord(t) = ord(t_solve) + modelSolveRules(m_solve, 't_forecastLength')] = -1;
msft(m_solve, s, f, t) = no;
msft(m_solve, 's000', f, t) = mft(m_solve,f,t);
msft(m_solve, 's000', fRealization(f), t)${ [ord(t) >= ord(t_solve) + modelSolveRules(m_solve, 't_forecastLength')]
and [ord(t) <= ord(t_solve) + modelSolveRules(m_solve, 't_horizon')]
and mf(m_solve, f)
} = yes;
ft(f,t) = no;
ft(f,t) = mft(m_solve, f, t);
ft_realized(f,t) = no;
ft_realized(f,t)$[fRealization(f) and ord(t) = ord(t_solve)] = yes;
pf(ft(f,t))$(ord(t) eq ord(t_solve) + 1) = 1 - ord(f);
pt(t)$[ ord(t) > ord(t_solve) and ord(t) <= ord(t_solve) + modelSolveRules(m_solve, 't_horizon') ] = -1;
// Arbitrary value for energy in storage
p_storageValue(egs(etype, geo, storage), t)$sum(fRealization(f), ft(f,t)) = 50;
// PSEUDO DATA
ts_reserveDemand(resType, resDirection, bus, fRealization(f), t) = 50;
* --- Set variable limits (.lo and .up) ---------------------------------------
$$include 'inc/set_variable_limits.gms'
* --- Solve Model -------------------------------------------------------------
if (m_solve('schedule'),
solve schedule using lp minimizing v_obj;
);
* --- Output debugging inetypeation --------------------------------------------
$$ifi '%debug%' == 'yes'
execute_unload 'output/debug.gdx';
*--- Store results ------------------------------------------------------------
// Deterministic stage
loop(ft(fRealization(f), t),
* p_stoContent(f, t)$(p_data(storage, 'max_content') > 0)
* = v_stoContent(storage, f, t) / p_data(f, 'max_content');
r_gen(etype, unit, t) = sum(geo$egu(etype, geo, unit), v_gen.l(etype, geo, unit, f, t));
$iftheni.genTypes '%genTypes%' == 'yes'
r_elec_type(genType, t) = sum(g $genType_g(genType, unit),
v_gen.l('elec', unit, f, t));
$endif.genTypes
r_demand(bus, t)
= sum(eg('elec', bus), ts_energyDemand('elec', bus, f, t));
$ontext
r_transmission(h, from_bus, to_bus, t)
= v_transmission(from_bus, to_bus, t);
r_elecConsumption(h, consuming(elec))
= v_elecConsumption.l(elec, t);
r_elecPrice(h, bus)
= q_elecDemand.m(bus, f, t) / p_blockLength(t);
r_heat(h, heat) = v_heat.l(heat, t);
r_storageControl(h, storage)
= v_stoCharge.l(storage, f, t) / p_blockLength(t);
r_onlineCap(h, elec) = v_gen.l('elec', elec, fCentral(f), t);
r_onlineCap(h, unitMinLoad(elec)) = v_online.l(elec, f, t);
r_onlineCap(h, unitVG(elec)) = p_unitVG(elec, f, t);
r_onlineCap(h, unitVG(elec))
= v_elec.l(unitVG, t);
loop(step_hour(h, t),
r_stoContent(h, f)$p_data(f, 'max_content')
= r_stoContent(h - 1, f)
+ (r_storageControl(h, f)
+ ts_inflow(h, f)
+ sum(unit_storage(unitVG, f),
ts_inflow(h, unitVG)
)
) / p_data(f, 'max_content');
r_storageValue(h, f) = p_storageValue(f, t);
r_elecLoad(h, bus)
= sum(load_in_hub(load, bus), ts_elecLoad(h, load));
);
$offtext
);
r_totalCost = r_totalCost + v_obj.l;
// Debugging results
$iftheni.debug '%debug%' == 'yes'
* putclose gdx;
* put_utility 'gdxout' / 'output\'m_solve.tl:0, '-', t_solve.tl:0, '.gdx';
* execute_unload
* $$include inc/debug_symbols.inc
* ;
$endif.debug
);
* Time independent results
$iftheni '%genTypes%' == 'yes'
r_capacity_type(genType)
= sum(g$genType_g(genType, g), p_data(g, 'max_power'));
$endif
* === Output ==================================================================
* Get model version
$echon "'StoSSch version' " > 'stosschver'
$call 'git describe --dirty=+ --always >> stosschver'
* Create metadata
set metadata(*) /
'User' '%sysenv.username%'
'Date' '%system.date%'
'Time' '%system.time%'
'GAMS version' '%system.gamsrelease%'
'GAMS system' '%system.gstring%'
$include 'stosschver'
/;
if(errorcount > 0, metadata('FAILED') = yes);
put f_info
put "***********************************************************************"/;
put "* MODEL RUN DETAILS *"/;
put "***********************************************************************"/;
loop(metadata,
put metadata.tl:20, metadata.te(metadata) /;
);
put /;
put "time (s)":> 21 /;
put "---------------------"/;
put "Compilation", system.tcomp:> 10 /;
put "Execution ", system.texec:> 10 /;
put "Total ", system.elapsed:> 10 /;
put /;
put "***********************************************************************"/;
put "* MODEL FEATURES *"/;
put "***********************************************************************"/;
loop(feature $active(feature),
put feature.tl:20, feature.te(feature):0 /;
);
put /;
f_info.nd = 0; // Set number of decimals to zero
put "Start time: ", modelSolveRules('schedule', 't_start')/;
put "Length of forecasts: ", modelSolveRules('schedule', 't_forecastLength')/;
put "Model horizon: ", modelSolveRules('schedule', 't_horizon')/;
put "Model jumps after solve: ", modelSolveRules('schedule', 't_jump')/;
put "Last time period to solve: ", modelSolveRules('schedule', 't_end')/;
put "Length of each time period: ", modelSolveRules('schedule', 't_interval')/;
put "Number of samples: ", modelSolveRules('schedule', 'samples')/;
putclose;
* -----------------------------------------------------------------------------
* Post-process results
$if exist 'postprocess.gms' $include 'postprocess.gms'
execute_unload 'output/results.gdx',
$$include 'inc/result_symbols.inc'
;
$ifi '%debug%' == 'yes'
execute_unload 'output/debug.gdx';
if(errorcount > 0, abort errorcount);
* === THE END =================================================================
* Parameters
p_stepLength
p_data
p_data2d
ts_cf
ts_inflow
ts_stoContent
ts_reserveDemand
* Variables
v_obj
v_stoCharge
v_gen
v_fuelUse
v_spill
v_stoContent
v_online
v_startup
v_transfer
v_resTransCapacity
v_reserve
* equations
q_obj
q_balance
q_resDemand
q_maxDownward
q_maxUpward
q_storageControl
q_storageDynamics
q_bindStorage
q_startup
q_bindOnline
q_fuelUse
q_conversion
q_outputRatioFixed
q_outputRatioConstrained
q_stoMinContent
q_stoMaxContent
q_maxHydropower
q_transferLimit
* Dummy variables
vq_gen
vq_resDemand
vq_stoCharge
equations
q_obj "Objective function"
q_balance(etype, geo, f, t) "Energy demand must be satisfied at each geographical location"
q_resDemand(resType, resDirection, geo, f, t) "Demand for each reserve type is greater than demand"
q_maxDownward(etype, geo, unit, f, t) "Downward commitments will not undercut minimum load or maximum elec. consumption"
q_maxUpward(etype, geo, unit, f, t) "Upward commitments will not exceed maximum available capacity or consumed power"
q_storageControl(etype, geo, storage, f, t) "Storage energy control"
q_storageDynamics(etype, geo, storage, f, t) "Dynamic equation for storages"
q_bindStorage(etype, geo, storage, modelType, f, t) "Couple storage contents for joining forecasts or for joining sample time periods"
q_startup(geo, unit, f, t) "Capacity started up is greater than the difference of online cap. now and in previous time step"
q_bindOnline(geo, unit, modelType, f, t) "Couple online variable for joining forecasts or for joining sample time periods"
q_fuelUse(geo, unit, fuel, f, t) "Use of fuels in units equals generation and losses"
* q_storageEnd(geo, storage, f, t) "Expected storage end content minus procured reserve energy is greater than start content"
q_conversion(geo, unit, f, t) "Conversion of energy between etypes of energy presented in the model, e.g. electricity consumption equals unitHeat generation times efficiency"
q_outputRatioFixed(etype, etype, geo, unit, f, t) "Fixed ratio between two etypes of energy output"
q_outputRatioConstrained(etype, etype, geo, unit, f, t) "Constrained ratio between two etypes of energy output; e.g. electricity generation is greater than c_V times unitHeat generation in extraction plants"
q_stoMinContent(etype, geo, storage, f, t) "Storage should have enough content at end of time step to deliver committed upward reserves"
q_stoMaxContent(etype, geo, storage, f, t) "Storage should have enough room to fit committed downward reserves"
q_maxHydropower(etype, geo, storage, f, t) "Sum of unitHydro generation in storage is limited by total installed capacity"
q_transferLimit(etype, geo, geo, f, t) "Transfer of energy and capacity reservations are less than the transfer capacity"
;
$setlocal def_penalty 10e6
Scalars
PENALTY "Default equation violation penalty" / %def_penalty% /
;
Parameters
PENALTY_BALANCE(etype) "Penalty on violating energy balance eq. (/MWh)"
PENALTY_RES(resType, resDirection) "Penalty on violating a reserve (/MW)"
;
PENALTY_BALANCE(etype) = %def_penalty%;
PENALTY_RES(resType, resDirection) = %def_penalty%;
* -----------------------------------------------------------------------------
q_obj ..
+ v_obj
=E=
+ sum(msft(m, s, f, t),
p_sProbability(s) *
p_fProbability(f) *
(
// Variable O&M costs
sum(egu(etype, geo, unit),
uData(etype, geo, unit, 'OaM_costs') *
$$ifi not '%rampSched%' == 'yes' p_stepLength(m, f, t) *
$$ifi '%rampSched%' == 'yes' (p_stepLength(m, f, t) + p_stepLength(m, f, t+1))/2 *
(
v_gen(etype, geo, unit, f, t)
)
)
// Fuel and emission costs
+ sum((geo, unitFuel, fuel)$(gu(geo, unitFuel) and unit_fuel(unitFuel, fuel, 'main')),
v_fuelUse(geo, unitFuel, fuel, f, t)
* ( sum{t_fuel$[ord(t_fuel) <= ord(t)],
ts_fuelPriceChangeGeo(fuel, geo, t_fuel) } // Fuel costs, sum initial fuel price plus all subsequent changes to the fuelprice
+ sum(emission, // Emission taxes
p_data2d(fuel, emission, 'emission_intensity') / 1e3
* p_data2d(emission, geo, 'emission_tax')
)
)
)
// Start-up costs
+ sum(egu(etype, geo, unitOnline),
+ {
+ v_startup(geo, unitOnline, f, t) // Cost of starting up
- sum(t_$mftStart(m, f, t_), 0.5 * v_online(geo, unitOnline, f, t_)) // minus value of avoiding startup costs before
- sum(t_$mftLastSteps(m, f, t_), 0.5 * v_online(geo, unitOnline, f, t_)) // or after the model solve
}
* {
// Startup variable costs
+ uData(etype, geo, unitOnline, 'startup_cost')
* uData(etype, geo, unitOnline, 'max_cap')
// Start-up fuel and emission costs
+ sum(unit_fuel(unitFuel, fuel, 'startup'),
+ uData(etype, geo, unitOnline, 'max_cap')
* uData(etype, geo, unitOnline, 'startup_fuelcons')
// Fuel costs for start-up fuel use
* ( + sum{t_fuel$[ord(t_fuel) <= ord(t)],
ts_fuelPriceChangeGeo(fuel, geo, t_fuel) }
// Emission taxes of startup fuel use
+ sum(emission,
p_data2d(emission, fuel, 'emission_intensity') / 1e3
* p_data2d(emission, geo, 'emission_tax')
)
)
)
}
)
) // p_sProbability(s) & p_fProbability(f)
) // msft(m, s, f, t)
// Value of energy storage change
- sum((mftLastForecast(m, f, t), mftStart(m, f_, t_)) $(active('storageValue')),
p_fProbability(f) *
sum(egs(etype, geo, storage),
p_storageValue(etype, geo, storage, t) *
(v_stoContent(etype, geo, storage, f, t) - v_stoContent(etype, geo, storage, f_, t_))
)
)
// Dummy variables
+ sum(msft(m, s, f, t), p_sProbability(s) * p_fProbability(f) * (
sum(inc_dec,
sum( eg(etype, geo), vq_gen(inc_dec, etype, geo, f, t) * p_stepLength(m, f, t) * PENALTY_BALANCE(etype) )
)
+ sum((resType, resDirection, bus),
vq_resDemand(resType, resDirection, bus, f, t)