Skip to content

GitLab

Commit 6708c8f7 authored by Stefan Englberger's avatar Stefan Englberger 🤠
Browse files

initial commit

parents
Show whitespace changes
Inline Side-by-side
.gitignore 0 → 100644
View file @ 6708c8f7
# Ignore everything
*
# But not these files...
!.gitignore
!LICENSE
!README.md
!logo.png
!script_00_main.m
!script_10_initialization.m
!script_11_load_profiles.m
!script_12_create_profiles.m
!script_13_fcr_bounds.m
!script_20_optimization.m
!script_21_constraints.m
!script_22_soh_calculation.m
!script_30_post_calculation.m
!script_31_plotting.m
\ No newline at end of file
LICENSE 0 → 100644
View file @ 6708c8f7
BSD 3-Clause License
Copyright (c) 2020, Stefan Englberger, Institute for Electrical
Energy Storage Technology, Technical University of Munich.
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
* Neither the name of the copyright holder nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
README.md 0 → 100644
View file @ 6708c8f7
# Multi-use optimization tool for energy storage systems (mu_opt)
This stand-alone tool generates the optimized operating strategy for energy storage systems. Due to its flexible structure, this optimization framework enables the simultaneous use of multiple storage applications (dynamic multi-use).
The tool was created by [Stefan Englberger](https://www.ei.tum.de/en/ees/staff/mitarbeiter-liste/englberger-stefan/) at the Institute for Electrical Energy Storage Technology at the Technical University of Munich. For efficient optimization, the tool was designed in Mathworks Matlab and also supports the Gurobi solver.
### **Highlights & Features:**
#### Maximization of profit from multi-use operation
- Dynamic allocation of energy and power of physical energy storage system
- Segmentation into virtual energy and power partitions
- Behind-the-meter and front-the-meter applications
#### Optimization
- Mixed-integer linear programming
- Problem-based optimization framework
- Linearization of non-linear internal processes
#### Supporting state-of-the-art applications
- Frequency containment reserve
- Peak shaving
- Self-consumption improvement
- Spot market trading
#### Technical considerations include
- Component efficiencies and system peripherals
- Degradation model for well-established cell chemistry
- Self-discharge of storage system
#### Regulatory constraints include
- System topology
- Application's degrees of freedom
- Spot market regulations
### **How to cite:**
[Stefan Englberger, Andreas Jossen, Holger Hesse. (2020). Unlocking the Potential of Battery Storage With the Dynamic Stacking of Multiple Applications. Cell Reports Physical Science, 1(11).](https://doi.org/10.1016/j.xcrp.2020.100238)
![alt text](https://gitlab.lrz.de/open-ees-ses/mu_opt/raw/master/logo.png "mu_opt / Englberger")
\ No newline at end of file
logo.png 0 → 100644
View file @ 6708c8f7
logo.png

133 KB

script_00_main.m 0 → 100644
View file @ 6708c8f7
%% Script Description
% Coded and provided by:
% Stefan Englberger (stefan.englberger@tum.de)
% Institute for Electrical Energy Storage Technology
% Technical University of Munich
%
% 2020-10-10 Stefan Englberger (author)
%% How to cite:
% Stefan Englberger, Andreas Jossen, Holger Hesse. (2020). Unlocking the Potential of Battery
% Storage With the Dynamic Stacking of Multiple Applications. Cell Reports Physical Science, 1(11).
% https://doi.org/10.1016/j.xcrp.2020.100238
%% clear workspace
clc; clear; close all;
%% run optimization
try
Input.scenario_folder = '';
Input.scenario_name = datestr(now,'yyyymmddHHMMSS');
disp(['Simulation ',Input.scenario_folder,Input.scenario_name,' started.']);
run('script_10_initialization.m'); % Data initialization
run('script_20_optimization.m'); % Run optimization
run('script_30_post_calculation.m'); % Concatenate and save files
disp(['Simulation ',Input.scenario_folder,Input.scenario_name,' was successful.']);
catch
try config_step = uint32(optimization_loop); catch; config_step = 0; end
try config_steps = uint32(Input.simulation_time/Input.rolling_horizon); catch; config_steps = 0; end
disp(['WARNING: Simulation ',Input.scenario_folder,Input.scenario_name,' was aborted at [',num2str(config_step),' / ',num2str(config_steps),'].']);
end
%% clear workspace
clearvars -except Input* Optimization*;
\ No newline at end of file
script_10_initialization.m 0 → 100644
View file @ 6708c8f7
%% Script Description
% Coded and provided by:
% Stefan Englberger (stefan.englberger@tum.de)
% Institute for Electrical Energy Storage Technology
% Technical University of Munich
%
% 2020-10-10 Stefan Englberger (author)
%% How to cite:
% Stefan Englberger, Andreas Jossen, Holger Hesse. (2020). Unlocking the Potential of Battery
% Storage With the Dynamic Stacking of Multiple Applications. Cell Reports Physical Science, 1(11).
% https://doi.org/10.1016/j.xcrp.2020.100238
% start timer
tic;
%% general settings
Input.sample_time = 5 / 60; % [h] length of sample time
Input.optimization_time = 24; % [h] length of optimization period
Input.rolling_horizon = 2; % [h] length of rolling horizon period
Input.simulation_time = 365 * 24; % [h] length of simulation period
Input.simulation_start = 0 * 24; % [h] start time of simulation (effects profiles)
path_profile = 'profiles.mat'; % path of profiles
path_results = 'results\'; % path of results folder
Input.intermediate_data = 0; % save intermediate results all x optimizatinon loops
Input.optimization_full = false; % save full optimization data
Input.soh_caluclation = true; % state of health calculation
% rolling horizon has to be a multiple of the sample time
if mod(Input.rolling_horizon,Input.sample_time) ~= 0
error('Rolling horizon must be a multiple of the sample time.');
end
% optimization time has to be a multiple of the sample time
if mod(Input.optimization_time,Input.sample_time) ~= 0
error('Optimization time must be a multiple of the sample time.');
end
% simulation time has to be a multiple of the rolling horizon
if mod(Input.simulation_time,Input.rolling_horizon) ~= 0
Input.simulation_time = Input.rolling_horizon * ceil(Input.simulation_time/Input.rolling_horizon);
warning(['Simulation time must be a multiple of the rolling horizon period.',newline,' Simulation time corrected to ',num2str(Input.simulation_time/24),' days.']);
end
% enable/disable behind-the-meter (btm) and front-the-meter (ftm) applications
Input.fcr = true; % conduct frequency containment reserve (fcr, ftm)
Input.ps = true; % conduct peak shaving (ps, btm)
Input.qc = true; % conduct reactive power compensation
Input.sci = true; % conduct self-consumption increase (sci, btm)
Input.smt = true; % conduct spot market trading (smt, ftm)
% define the input folder and profiles
Input.number_price_fcr = 13; % number of profile for fcr price
Input.type_price_smt = 'idc'; % type of spot market ('dam' or 'idc')
Input.number_price_smt_low = 37; % number of profile for spot market price (low price signal)
Input.number_price_smt_high = 38; % number of profile for spot market price (high price signal)
Input.number_power_generator = 5; % number of profile for active power from generator
Input.number_power_load_active = 212; % number of profile for active load power
Input.number_power_load_reactive = 1; % number of profile for reactive load power
Input.number_frequency = 8; % number of profile for frequency
% optimization options
Input.max_optimization_time = 900; % [s] maximum time per optimization
Input.max_relative_optimization_gap = 5e-3; % [1] maximum relative gap per optimization
%% energy storage system (ess)
Input.ESS.state_of_health_initial = 1.0; % [1] initial state of health of storage system
Input.ESS.energy_nominal = 1340; % [kWh] nominal energy content of storage system
Input.ESS.inverter_amount = 24; % integer amount of inverters
Input.ESS.power_apparent = 65.5; % [kVA] rated apparent power per inverter
Input.ESS.power_active = 0.8 * Input.ESS.power_apparent; % [kW] rated active power per inverter
Input.ESS.power_reactive = 0.3 * Input.ESS.power_apparent; % [kvar] rated reactive power per inverter
Input.ESS.state_of_charge_min = 0.00; % [1] minimum state of charge (lower bound)
Input.ESS.state_of_charge_max = 1.00; % [1] maximum state of charge (upper bound)
Input.ESS.state_of_charge_initial = 0.50; % [1] state of charge at the begin of simulation
Input.ESS.state_of_charge_end = 0.50; % [1] state of charge at the end of simulation
Input.ESS.efficiency_battery = 0.97; % [1] efficiency of battery cells
Input.ESS.efficiency_inverter = 0.96; % [1] efficiency of power electronics
Input.ESS.efficiency_transformer = 1.00; % [1] efficiency of transformer
Input.ESS.efficiency_charge = prod([Input.ESS.efficiency_battery,Input.ESS.efficiency_inverter,Input.ESS.efficiency_transformer]); % [1] efficiency during charging
Input.ESS.efficiency_discharge = prod([Input.ESS.efficiency_battery,Input.ESS.efficiency_inverter,Input.ESS.efficiency_transformer]); % [1] efficiency during discharging
Input.ESS.power_active_charge = min(1.00*Input.ESS.energy_nominal/1,Input.ESS.inverter_amount*Input.ESS.power_active); % [kW] maximum active power during charging of ESS
Input.ESS.power_active_discharge = min(1.00*Input.ESS.energy_nominal/1,Input.ESS.inverter_amount*Input.ESS.power_active); % [kW] maximum active power during discharging of ESS
Input.ESS.allocation_switching_time = 1.0; % [h] switching time of power and energy allocation
Input.ESS.temperature_cell = 20 + 273.15; % [K] cell temperature
Input.ESS.capacity_cell_nominal = 2.05; % [Ah] nominal cell capacity
Input.ESS.energy_self_discharge = 0.015/(30.5*24/Input.sample_time) * Input.ESS.energy_nominal; % self-discharge depending on nominal energy content
Input.ESS.energy_shift = [-0.0 , 0.0] * Input.ESS.energy_nominal; % maximum energy exchange between ftm and btm partition
Input.ESS.support_power = 0.00 * Input.ESS.energy_nominal; % [kW] additional power consumption for support energy of storage system
Input.ESS.invest_system = 380 * Input.ESS.energy_nominal; % [EUR] invest costs for storage system
Input.ESS.efc_capability = 3000; % [1] equivalent full cycles capability until system's end of life
%% price signals
% https://www.bdew.de/service/daten-und-grafiken/strompreis-fuer-die-industrie/
Input.price_base_distribution = 0.04300; % [EUR/kWh] base price for electricity distribution
Input.price_power_charges = 100.0; % [EUR/kW] power related grid charges
Input.price_grid_charges = 0.00800; % [EUR/kWh] energy related grid charges
Input.price_concession_fee = 0.00110; % [EUR/kWh] concession fee in Germany (http://www.gesetze-im-internet.de/kav/__2.html)
Input.price_eeg_surcharge = 0.06405; % [EUR/kWh] eeg surcharge in Germany (https://www.bundesnetzagentur.de/DE/Sachgebiete/ElektrizitaetundGas/Verbraucher/PreiseRechnTarife/preiseundRechnungen-node.html)
Input.price_kwkg_surcharge = 0.00280; % [EUR/kWh] kwkg surcharge in Germany (http://www.gesetze-im-internet.de/kwkg_2016/)
Input.price_strom_nev_surcharge = 0.00200; % [EUR/kWh] strom nev surcharge (paragraph 19 Strom NEV) in Germany (https://www.gesetze-im-internet.de/stromnev/__19.html)
Input.price_offsu_surcharge = 0.00416; % [EUR/kWh] offshore liability levy in Germany (https://www.gesetze-im-internet.de/enwg_2005/__17f.html)
Input.price_ablav_surcharge = 0.00005; % [EUR/kWh] surcharge for interruptible loads in Germany (Abschaltverordnung AbLaV) (https://www.gesetze-im-internet.de/ablav_2016/__18.html)
Input.price_electricity_tax = 0.01537; % [EUR/kWh] electricity tax in Germany
Input.price_trading_fee = 0.0001425; % [EUR/kWh] electricity trading fee
%% import profiles and application-specific data
% import profiles
Input.peak_generation = 0; % [kW] peak power of generation profile
Input.consumption = 20e6; % [kWh] annual energy consumption
Input.curtailment_btm_max = Inf; % [kW] maximum curtailment power at btm partition
% front-the-meter (ftm) settings
Input.curtailment_ftm_max = 0; % [kW] maximum curtailment power at ftm partition
run('script_11_load_profiles.m'); % load existing profiles
run('script_12_create_profiles.m'); % create necessary profiles
% frequency containment reserve (fcr) settings
Input.FCR.objective_weight = 1; % [1] relative weighting of the application in the objective function
Input.FCR.period = 4; % [h] length of minimum period for fcr provision (must be multiple of sample time)
Input.FCR.minimum_power = 100; % [kW] minimum power for fcr provision
Input.FCR.reserve_time = 15/60; % [h] reserved time for additional energy criteria
Input.FCR.reserve_power = 0.25; % [1] reserved power for additional power criteria
Input.FCR.frequency_nominal = 50.0; % [Hz] nominal frequency
Input.FCR.frequency_deadband = 0.01; % [Hz] frequency dead band of fcr provision
Input.FCR.frequency_droop = 0.4/100; % [1] frequency droop
Input.FCR.overfulfillment = 1.2; % [1] maximum overfulfillment
Input.FCR.energy_difference = 1.0 * Input.ESS.energy_nominal; % [kWh] maximum energy difference at start and end of provision period
Input.FCR.revenue_share_min = 0.60; % [1] minimum share of revenue from fcr provision
Input.FCR.revenue_share_max = 0.85; % [1] maximum share of revenue from fcr provision
Input.FCR.revenue_share_max_power = 1000; % [kW] power threshold to reach maximum share of revenue
run('script_12_create_profiles.m'); % create necessary profiles
run('script_13_fcr_bounds.m'); % create lower and upper bound of fcr power considering the specified degrees of freedom
% peak shaving (ps) settings
Input.PS.objective_weight = 1; % [1] relative weighting of the application in the objective function
Input.PS.price = Input.price_power_charges; % [EUR/kW] electricity tariff for peak power
Input.PS.billing_period = 365 * 24; % [h] billing period for peak power threshold
Input.PS.power_threshold_initial = Input.consumption / 8760; % [kW] initial peak shave threshold (for each billing period)
Input.PS.reset_threshold = true; % reset power threshold to initial value for each billing period
% self-consumption increase (sci) settings
Input.SCI.objective_weight = 1; % [1] relative weighting of the application in the objective function
Input.SCI.grid_tariff = sum([Input.price_base_distribution,Input.price_grid_charges,Input.price_electricity_tax,...
Input.price_eeg_surcharge,Input.price_kwkg_surcharge,Input.price_concession_fee,...
Input.price_strom_nev_surcharge,Input.price_offsu_surcharge,Input.price_ablav_surcharge]); % [EUR/kWh] electricity tariff for energy purchase
Input.SCI.grid_remuneration = 0.10610; % [EUR/kWh] remuneration price for energy sale
Input.SCI.price_selfconsumed = 0.4 * Input.price_eeg_surcharge; % [EUR/kWh] costs for self-consumed electricity
Input.SCI.feedin_limit = 1e3; % [kW] feed-in limitation for power flow to grid
% spot market trading (smt) settings
Input.SMT.objective_weight = 1; % [1] relative weighting of the application in the objective function
Input.SMT.period = 15/60; % [h] length of minimum period (must be multiple of sample time)
Input.SMT.minimum_power = 100; % [kW] minimum power for transaction
Input.SMT.charges_purchase = Input.price_trading_fee + Input.price_electricity_tax + Input.price_concession_fee; % [EUR/kWh] charges for purchased energy
Input.SMT.charges_sell = Input.price_trading_fee; % [EUR/kWh] charges for sold energy
run('script_12_create_profiles.m'); % create necessary profiles
%% check application-related errors
% rolling horizon has to be a multiple of the spot market provision period
if mod(Input.rolling_horizon,Input.SMT.period) ~= 0
error('Rolling horizon must be a multiple of the SMT provision period.');
end
% optimization time has to be a multiple of the allocation switching time
if mod(Input.optimization_time,Input.ESS.allocation_switching_time) ~= 0
error('Optimization time must be a multiple of the allocation switching time.');
end
% optimization time has to be a multiple of the fcr provision period
if mod(Input.optimization_time,Input.FCR.period) ~= 0
error('Optimization time must be a multiple of the FCR provision period.');
end
% optimization time has to be a multiple of the spot market provision period
if mod(Input.optimization_time,Input.SMT.period) ~= 0
error('Optimization time must be a multiple of the SMT provision period.');
end
% allocation switching time has to be a multiple of the sample time
if mod(Input.ESS.allocation_switching_time,Input.sample_time) ~= 0
error('Allocation switching time must be a multiple of the sample time.');
end
% fcr provision period has to be a multiple of the sample time
if mod(Input.FCR.period,Input.sample_time) ~= 0
error('FCR provision period must be a multiple of the sample time.');
end
% smt provision period has to be a multiple of the sample time
if mod(Input.FCR.period,Input.sample_time) ~= 0
error('Spot market provision period must be a multiple of the sample time.');
end
\ No newline at end of file
script_11_load_profiles.m 0 → 100644
View file @ 6708c8f7
%% Script Description
% Coded and provided by:
% Stefan Englberger (stefan.englberger@tum.de)
% Institute for Electrical Energy Storage Technology
% Technical University of Munich
%
% 2020-10-10 Stefan Englberger (author)
%% How to cite:
% Stefan Englberger, Andreas Jossen, Holger Hesse. (2020). Unlocking the Potential of Battery
% Storage With the Dynamic Stacking of Multiple Applications. Cell Reports Physical Science, 1(11).
% https://doi.org/10.1016/j.xcrp.2020.100238
% define the length of profiles
profile_length = ( Input.simulation_time + Input.optimization_time - Input.rolling_horizon ) / Input.sample_time;
% load profiles file if available
if exist(path_profile,'file') == 2
load(path_profile, ...
['profile_price_fcr_',sprintf('%03d',Input.number_price_fcr)], ...
['profile_price_',Input.type_price_smt,'_',sprintf('%03d',Input.number_price_smt_low)], ...
['profile_price_',Input.type_price_smt,'_',sprintf('%03d',Input.number_price_smt_high)], ...
['profile_power_gen_',sprintf('%03d',Input.number_power_generator)], ...
['profile_power_active_load_',sprintf('%03d',Input.number_power_load_active)], ...
['profile_power_reactive_load_',sprintf('%03d',Input.number_power_load_reactive)], ...
['profile_frequency_',sprintf('%03d',Input.number_frequency)] ...
);
% fcr price signal [EUR/kWh]
profile_price_fcr = eval(['profile_price_fcr_',sprintf('%03d',Input.number_price_fcr),';']);
Input.price_fcr = repelem(double(profile_price_fcr),60,1); % [EUR/kWh] % call fcr price profile
Input.price_fcr = repmat(Input.price_fcr,ceil((Input.simulation_start+Input.simulation_time+Input.optimization_time-Input.rolling_horizon)/24/365),1);
dummy = zeros(profile_length,1);
for i = uint32(1:profile_length); dummy(i) = mean(Input.price_fcr(Input.simulation_start*(60*60)+(i-1)*(Input.sample_time*(60*60))+1:Input.simulation_start*(60*60)+i*(Input.sample_time*(60*60)))); end; clear i;
Input.price_fcr = dummy; clear dummy;
% low spot market price signal [EUR/kWh]
profile_price_smt_low = eval(['profile_price_',Input.type_price_smt,'_',sprintf('%03d',Input.number_price_smt_low),';']);
Input.price_smt_low = repelem(double(profile_price_smt_low),60,1); % [EUR/kWh] % call spot market price profile
Input.price_smt_low = repmat(Input.price_smt_low,ceil((Input.simulation_start+Input.simulation_time+Input.optimization_time-Input.rolling_horizon)/24/365),1);
dummy = zeros(profile_length,1);
for i = uint32(1:profile_length); dummy(i) = mean(Input.price_smt_low(Input.simulation_start*(60*60)+(i-1)*(Input.sample_time*(60*60))+1:Input.simulation_start*(60*60)+i*(Input.sample_time*(60*60)))); end; clear i;
Input.price_smt_low = dummy; clear dummy;
% high spot market price signal [EUR/kWh]
profile_price_smt_high = eval(['profile_price_',Input.type_price_smt,'_',sprintf('%03d',Input.number_price_smt_high),';']);
Input.price_smt_high = repelem(double(profile_price_smt_high),60,1); % [EUR/kWh] % call spot market price profile
Input.price_smt_high = repmat(Input.price_smt_high,ceil((Input.simulation_start+Input.simulation_time+Input.optimization_time-Input.rolling_horizon)/24/365),1);
dummy = zeros(profile_length,1);
for i = uint32(1:profile_length); dummy(i) = mean(Input.price_smt_high(Input.simulation_start*(60*60)+(i-1)*(Input.sample_time*(60*60))+1:Input.simulation_start*(60*60)+i*(Input.sample_time*(60*60)))); end; clear i;
Input.price_smt_high = dummy; clear dummy;
% correct spot market price signal if smt is disabled (but scheduled transactions during fcr are enabled)
if ~Input.smt
Input.price_smt_low(:) = mean(mean([Input.price_smt_low,Input.price_smt_high],2));
Input.price_smt_high = Input.price_smt_low;
end
% active power supply [kW]
profile_power_gen = eval(['profile_power_gen_',sprintf('%03d',Input.number_power_generator),';']);
Input.power_generation = repelem(double(profile_power_gen),60,1); % call generation profile
Input.power_generation = Input.power_generation / max(Input.power_generation) * Input.peak_generation; % scaling of profile
Input.power_generation = repmat(Input.power_generation,ceil((Input.simulation_start+Input.simulation_time+Input.optimization_time-Input.rolling_horizon)/24/365),1);
dummy = zeros(profile_length,1);
for i = uint32(1:profile_length); dummy(i) = mean(Input.power_generation(Input.simulation_start*(60*60)+(i-1)*(Input.sample_time*(60*60))+1:Input.simulation_start*(60*60)+i*(Input.sample_time*(60*60)))); end; clear i;
Input.power_generation = dummy; clear dummy;
% active power demand [kW]
profile_power_active_load = eval(['profile_power_active_load_',sprintf('%03d',Input.number_power_load_active),';']);
Input.power_active_load = repelem(double(profile_power_active_load),60,1); % call load power profile
Input.power_active_load = Input.power_active_load / sum(Input.power_active_load) * Input.consumption * 3.6e3; % scaling of profile
Input.power_active_load = repmat(Input.power_active_load,ceil((Input.simulation_start+Input.simulation_time+Input.optimization_time-Input.rolling_horizon)/24/365),1);
dummy = zeros(profile_length,1);
for i = uint32(1:profile_length); dummy(i) = mean(Input.power_active_load(Input.simulation_start*(60*60)+(i-1)*(Input.sample_time*(60*60))+1:Input.simulation_start*(60*60)+i*(Input.sample_time*(60*60)))); end; clear i;
Input.power_active_load = dummy; clear dummy;
% reactive power [kvar]
profile_power_reactive_load = eval(['profile_power_reactive_load_',sprintf('%03d',Input.number_power_load_reactive),';']);
Input.power_reactive_load = repelem(double(profile_power_reactive_load),60,1); % call load power profile
Input.power_reactive_load = Input.power_reactive_load / sum(Input.power_reactive_load) * Input.consumption / 6; % scaling of profile; 6 = ratio between active and reactive profile (in this particular case)
Input.power_reactive_load = repmat(Input.power_reactive_load,ceil((Input.simulation_start+Input.simulation_time+Input.optimization_time-Input.rolling_horizon)/24/365),1);
dummy = zeros(profile_length,1);
for i = uint32(1:profile_length); dummy(i) = mean(Input.power_reactive_load(Input.simulation_start*(60*60)+(i-1)*(Input.sample_time*(60*60))+1:Input.simulation_start*(60*60)+i*(Input.sample_time*(60*60)))); end; clear i;
Input.power_reactive_load = dummy; clear dummy;
% frequency signal [Hz]
profile_frequency = eval(['profile_frequency_',sprintf('%03d',Input.number_frequency),';']);
Input.frequency = double(profile_frequency); % [Hz] % call frequency profile
Input.frequency = repmat(Input.frequency,ceil((Input.simulation_start+Input.simulation_time+Input.optimization_time-Input.rolling_horizon)/24/365),1); % [Hz] % call frequency profile
dummy = zeros(profile_length,1);
for i = uint32(1:profile_length); dummy(i) = mean(Input.frequency(Input.simulation_start*(60*60)+(i-1)*(Input.sample_time*(60*60))+1:Input.simulation_start*(60*60)+i*(Input.sample_time*(60*60)))); end; clear i;
Input.frequency = dummy; clear dummy;
end
% clear unused variables
clear profile_* path_profile;
\ No newline at end of file
script_12_create_profiles.m 0 → 100644
View file @ 6708c8f7
%% Script Description
% Coded and provided by:
% Stefan Englberger (stefan.englberger@tum.de)
% Institute for Electrical Energy Storage Technology
% Technical University of Munich
%
% 2020-10-10 Stefan Englberger (author)
%% How to cite:
% Stefan Englberger, Andreas Jossen, Holger Hesse. (2020). Unlocking the Potential of Battery
% Storage With the Dynamic Stacking of Multiple Applications. Cell Reports Physical Science, 1(11).
% https://doi.org/10.1016/j.xcrp.2020.100238
% define the length of profiles
profile_length = ( Input.simulation_time + Input.optimization_time - Input.rolling_horizon ) / Input.sample_time;
% fcr price signal [EUR/kWh]
if isfield(Input,'price_fcr') == 0 && isfield(Input,'FCR') == 1 && isfield(Input.FCR,'period') == 1
dummy = repelem(smoothdata(normrnd(0.012,0.003,[ceil(profile_length*Input.sample_time/Input.FCR.period),1]),'gaussian',1/Input.FCR.period*12),Input.FCR.period/Input.sample_time,1);
Input.price_fcr = dummy(mod(Input.simulation_start,Input.FCR.period)/Input.sample_time+1:mod(Input.simulation_start,Input.FCR.period)/Input.sample_time+profile_length);
clear dummy;
end
% spot market price signals [EUR/kWh]
if (isfield(Input,'price_smt_low') == 0 || isfield(Input,'price_smt_high') == 0) && isfield(Input,'SMT') == 1 && isfield(Input.SMT,'period') == 1
dummy = repelem(smoothdata(normrnd(0.045,0.03,[ceil(profile_length*Input.sample_time/Input.SMT.period),1]),'gaussian',Input.SMT.period*12),Input.SMT.period/Input.sample_time,1);
Input.price_smt_low = dummy(mod(Input.simulation_start,Input.SMT.period)/Input.sample_time+1:mod(Input.simulation_start,Input.SMT.period)/Input.sample_time+profile_length);
if Input.number_price_smt_low == Input.number_price_smt_high
Input.price_smt_high = Input.price_smt_low;
else
dummy = repelem(smoothdata(normrnd(0.045,0.03,[ceil(profile_length*Input.sample_time/Input.SMT.period),1]),'gaussian',Input.SMT.period*12),Input.SMT.period/Input.sample_time,1);
Input.price_smt_high = dummy(mod(Input.simulation_start,Input.SMT.period)/Input.sample_time+1:mod(Input.simulation_start,Input.SMT.period)/Input.sample_time+profile_length);
clear dummy;
dummy = [Input.price_smt_low,Input.price_smt_high];
Input.price_smt_low = min(dummy,[],2);
Input.price_smt_high = max(dummy,[],2);
clear dummy;
end
% correct spot market price signal if smt is disabled (but scheduled transactions during fcr are enabled)
if ~Input.smt
Input.price_smt_low(:) = mean(mean([Input.price_smt_low,Input.price_smt_high],2));
Input.price_smt_high = Input.price_smt_low;
end
end
% active power supply [kW]
if isfield(Input,'power_generation') == 0 && isfield(Input,'peak_generation') == 1
dummy = ceil(profile_length/24*Input.sample_time)+1;
dummy = sin(linspace(0,2*pi*dummy,dummy*24/Input.sample_time)).*abs(sin(linspace(0,2*pi*dummy,dummy*24/Input.sample_time)));
dummy(dummy<0) = 0; dummy = circshift(dummy,6/Input.sample_time); dummy = transpose(dummy);
dummy = (0.7+0.3*rand(size(dummy))).*dummy; dummy(dummy<0) = 0;
dummy = dummy/max(dummy)*Input.peak_generation;
Input.power_generation = dummy(mod(Input.simulation_start,24)/Input.sample_time+1:mod(Input.simulation_start,24)/Input.sample_time+profile_length);
clear dummy;
end
% active power demand [kW]
if isfield(Input,'power_active_load') == 0 && isfield(Input,'consumption') == 1
dummy = 0.3+0.7*rand(profile_length,1);
dummy = movmean(dummy,6/Input.sample_time,'omitnan').^5;
Input.power_active_load = dummy/sum(dummy)*Input.consumption*(profile_length*Input.sample_time)/(365*24)/Input.sample_time;
clear dummy;
end
% reactive power [kvar]
if isfield(Input,'power_reactive_load') == 0 && isfield(Input,'consumption') == 1
dummy = 0.3+0.7*rand(profile_length,1);
dummy = movmean(dummy,6/Input.sample_time,'omitnan').^5;
Input.power_reactive_load = dummy/sum(dummy)*Input.consumption/6*(profile_length*Input.sample_time)/(365*24)/Input.sample_time; % scaling of profile; 6 = ratio between active and reactive profile (in this particular case)
clear dummy;
end
% frequency signal [Hz]
if isfield(Input,'frequency') == 0 && isfield(Input,'FCR') == 1 && isfield(Input.FCR,'frequency_nominal') == 1
dummy = normrnd(Input.FCR.frequency_nominal,0.02,[profile_length,1]);
dummy = movmean(dummy,0.25/Input.sample_time,'omitnan');
Input.frequency = dummy;
clear dummy;
end
% clear unused variables
clear profile_*;
\ No newline at end of file
script_13_fcr_bounds.m 0 → 100644
View file @ 6708c8f7
%% Script Description
% Coded and provided by:
% Stefan Englberger (stefan.englberger@tum.de)
% Institute for Electrical Energy Storage Technology
% Technical University of Munich
%
% 2020-10-10 Stefan Englberger (author)
%% How to cite:
% Stefan Englberger, Andreas Jossen, Holger Hesse. (2020). Unlocking the Potential of Battery
% Storage With the Dynamic Stacking of Multiple Applications. Cell Reports Physical Science, 1(11).
% https://doi.org/10.1016/j.xcrp.2020.100238
%% calculation of fcr bounds (minimum and maximum)
% define input variables
frequency_delta = Input.frequency - Input.FCR.frequency_nominal; % [Hz]
fcr_power_min = zeros(size(Input.frequency));
fcr_power_max = zeros(size(Input.frequency));
% calculate power signal of fcr considering the frequency
for i = uint32(1:length(Input.frequency))
if frequency_delta(i) >= 0
if abs(frequency_delta(i)) <= Input.FCR.frequency_deadband
fcr_power_min(i) = 0;
else
fcr_power_min(i) = max(-1, -(frequency_delta(i) * 1) / (Input.FCR.frequency_nominal * Input.FCR.frequency_droop));
end
fcr_power_max(i) = max(-1, -(frequency_delta(i) * 1) / (Input.FCR.frequency_nominal * Input.FCR.frequency_droop)) * Input.FCR.overfulfillment;
else
if abs(frequency_delta(i)) <= Input.FCR.frequency_deadband
fcr_power_min(i) = 0;
else
fcr_power_min(i) = min(1, -(frequency_delta(i) * 1) / (Input.FCR.frequency_nominal * Input.FCR.frequency_droop));
end
fcr_power_max(i) = min(1, -(frequency_delta(i) * 1) / (Input.FCR.frequency_nominal * Input.FCR.frequency_droop)) * Input.FCR.overfulfillment;
end
end
% create dummy variables to transform the relative extremes to the absolute minimum and absolute maximum
dummy_min = min(fcr_power_min,fcr_power_max);
dummy_max = max(fcr_power_min,fcr_power_max);
% transformation of relative extremes to the absolute minimum and absolute maximum values
fcr_power_min = dummy_min;
fcr_power_max = dummy_max;
% clear unused variables
clear dummy_* frequency_delta i;
\ No newline at end of file
script_20_optimization.m 0 → 100644
View file @ 6708c8f7
%% Script Description
% Coded and provided by:
% Stefan Englberger (stefan.englberger@tum.de)
% Institute for Electrical Energy Storage Technology
% Technical University of Munich
%
% 2020-10-10 Stefan Englberger (author)
%% How to cite:
% Stefan Englberger, Andreas Jossen, Holger Hesse. (2020). Unlocking the Potential of Battery
% Storage With the Dynamic Stacking of Multiple Applications. Cell Reports Physical Science, 1(11).
% https://doi.org/10.1016/j.xcrp.2020.100238
% add the path of the gurobi solver, if available
if exist('gurobi','file') == 3
gurobipath = fileparts(which('gurobi'));
addpath([gurobipath(1:end-6),'examples\matlab']);
else
warning('It is recommended to use the Gurobi solver.')
promptMessage = sprintf('The Gurobi solver was not found in the MATLAB environment.\nDo you want to continue anyway?');
button = questdlg(promptMessage, 'Continue without Gurobi solver', 'Continue', 'Cancel', 'Continue');
if strcmpi(button, 'Cancel')
error('Simulation was aborted by user.');
end
end
% clear unused variables