Hydro-Thermal Scheduling
The source codes can be downloaded from here and is based on the course notes of LINMA2491 on SDDP at UCLouvain.Problem statement
We apply SDDP on the hydro-thermal scheduling problem. Suppose you need to take actions on two devices:- a thermal production unit is available at a marginal cost of 25 $/MWh and a capacity of 500 MW
- a hydroelectric dam that can store up to 1000 MWh of energy per hour
- independent uniformly distributed rainfall between time stages,
- an additive auto-regressive model of rainfall,
- a multiplicative auto-regressive model of rainfall.
Rainfall model
Consider a model of rainfall over 24 stages, where each stage represents a month. The three alternative models of rainfall (uniform, auto-regressive additive, and auto-regressive multiplicative) are developed by using a lattice which obeys serial independence. This lattice is built by defining a random variable \(w_t\) which is uniformly distributed in the set \(\{1, 2,\dots, 20\}\). Suppose that the rainfall follows a uniform distribution over the interval \([0, 1000]\) MW. Then this model of rainfall can be approximated by assigning to each \( (t, k) \)of the lattice the following rainfall value: $$ R_{t,k} = 50 w_{t,k}.$$ Now consider the alternative case where rainfall is assumed to follow an additive auto-regressive process, according to the following model: $$ R_t = \phi R_{t-1} + \sigma\epsilon_t. $$ where \(\phi\) is an auto-regressive parameter, \(\epsilon_t\) are independent identically distributed random variables that obey a standard normal distribution, and \(\sigma\) corresponds to the standard deviation of the random input. Consider the following choice of parameters for the additive auto-regressive model: \( c = 0 \) , \( \phi = 1 \), \(\sigma = 100 MW\), and an initial rainfall level of 500 MW. Then the 24-stage, 20-node lattice can be used for approximating additive auto-regressive rainfall, according to the following assignment of rainfall values in each node \((t, k)\) of the lattice: $$ R_{t,k} = R_{t-1} + 100 \Phi^{-1}(\frac{w_{t,k}}{20} - 0.025),$$ where \(\Phi\) is the cumulative distribution function of a standard normal variable and \(w_t\) is uniformly distributed in the interval \(\{1, 2,\dots, 20\}\). Finally, consider using the same lattice for approximating multiplicative auto-regressive rainfall of the following form: $$R_t = (c + \phi R_{t-1})\eta_t,$$ where \(\phi\) is an auto-regressive parameter, and \(\eta_t\) are independent identically distributed random variables that obey a normal distribution of mean 1 and standard deviation \(\sigma\). Consider the following choice of parameters for the multiplicative auto-regressive model: \(c = 0, \phi = 1, \sigma = 0.1 MW \), and an initial rainfall level of 500 MW. In order to approximate the multiplicative auto-regressive model on the 24-stage, 20-node lattice, the following assignment of rainfall values in each node \((t, k)\) of the lattice can be used: $$ R_{t,k} = R_{t-1} (1 + 0.01 \Phi^{-1}(\frac{w_{t,k}}{20} - 0.025)), $$Version 1, with basic rain model
% Hydro thermal planning problem , chapter 7.4
clc ; close all ; clear all ;
% Create lattice
H = 24 ;
nNodes = 20 ;
lattice = Lattice.latticeEasy(H, nNodes, @rain) ;
% Run SDDP
params = sddpSettings('algo.McCount',10, ...
'stop.iterationMax',10,...
'stop.pereiraCoef',0.01,...
'precise.computeEnd',true,...
'precise.count',200,...
'solver','linprog') ;
var.x = sddpVar(H) ;
var.p = sddpVar(H) ;
var.q = sddpVar(H) ;
var.l = sddpVar(H) ;
var.r = sddpVar(H) ;
lattice = compileLattice(lattice, @(scenario)nlds(scenario, var), params) ;
output = sddp(lattice, params) ;
plotOutput(output) ;
function [cntr, obj] = nlds(scenario, var)
t = scenario.getTime() ;
C = 25 ;
V = 1000 ;
D = 1000 ;
E = 1000;
PMax = 500;
InitialWater = 500;
InitialRain = 500;
x = var.x;
p = var.p;
q = var.q;
l = var.l;
r = var.r;
% Objective
obj = V * l(t) + C * p(t) ;
SupplyDemandBalance = p(t) + q(t) + l(t) >= D ;
Positivity = [x(t) >= 0, l(t) >= 0, p(t) >= 0, q(t) >= 0] ;
ThermalCapacityLimit = p(t) <= PMax ;
ReservoirLimit = x(t) <= E ;
Rain = r(t) == scenario.data;
if(scenario.time == 1)
HydroUtilizationLimit0 = q(t) <= InitialWater + scenario.data ;
HydroLevelDynamics0 = x(t) == InitialRain + scenario.data - q(t) ;
cntr = [SupplyDemandBalance , Positivity , ...
ThermalCapacityLimit , ReservoirLimit , ...
HydroUtilizationLimit0 , HydroLevelDynamics0 , Rain];
else
HydroUtilizationLimit = q(t) <= x(t-1) + scenario.data ;
HydroLevelDynamics = x(t) == x(t-1) + scenario.data - q(t) ;
cntr = [SupplyDemandBalance , Positivity , ...
ThermalCapacityLimit , ReservoirLimit , ...
HydroUtilizationLimit , HydroLevelDynamics , Rain];
end
end
function out = rain(t,i)
if t == 1
out = 500;
else
out = 50*i;
end
end
Version 2, with AR rain model
% Hydro thermal planning problem , chapter 7.4
clc ; close all ; clear all ;
% Create lattice
H = 24 ;
nNodes = 20 ;
lattice = Lattice.latticeEasy(H, nNodes, @rainDisturbanceG) ;
% Run SDDP
params = sddpSettings('algo.McCount', 150, ...
'stop.iterationMax', 10,...
'stop.pereiraCoef', 2,...
'precise.computeEnd', true,...
'precise.count', 200,...
'solver','gurobi') ;
var.x = sddpVar(H) ;
var.p = sddpVar(H) ;
var.q = sddpVar(H) ;
var.l = sddpVar(H) ;
var.r = sddpVar(H) ;
lattice = compileLattice(lattice, @(scenario)nldsAR(scenario, var), params) ;
output = sddp(lattice, params) ;
plotOutput(output) ;
function [cntr, obj] = nldsAR(scenario, var)
t = scenario.getTime() ;
C = 25 ;
V = 1000 ;
D = 1000 ;
E = 1000;
PMax = 500;
InitialWater = 500;
InitialRain = 700;
x = var.x;
p = var.p;
q = var.q;
l = var.l;
r = var.r;
% Objective
obj = V * l(t) + C * p(t);
SupplyDemandBalance = p(t) + q(t) + l(t) >= D ;
Positivity = [l(t) >= 0, p(t) >= 0];
%Posr = r(t) >= 0;
ThermalCapacityLimit = p(t) <= PMax ;
ReservoirLimit = x(t) <= E ;
if(scenario.time == 1)
HydroUtilizationLimit0 = q(t) <= InitialWater + r(t) ;
HydroLevelDynamics0 = x(t) <= InitialWater + r(t) - q(t) ;
RainDynamics0 = r(t) == max(InitialRain + scenario.data, 0);
cntr = [SupplyDemandBalance , Positivity, ...
ThermalCapacityLimit, ReservoirLimit, ...
HydroUtilizationLimit0, RainDynamics0, HydroLevelDynamics0];
else
HydroUtilizationLimit = q(t) <= x(t-1) + r(t) ;
HydroLevelDynamics = x(t) == x(t-1) + r(t) - q(t) ;
RainDynamics = r(t) == r(t-1) + scenario.data;
cntr = [SupplyDemandBalance , Positivity, ...
ThermalCapacityLimit, ReservoirLimit, ...
HydroUtilizationLimit, RainDynamics, HydroLevelDynamics];
end
end
function out = rainDisturbanceG(t,i)
if t == 1
out = 100*noiseG(11);
else
out = 100*noiseG(i);
end
end
Version 3, with ARM rain model
% Hydro thermal planning problem , chapter 7.4
clc ; close all ; clear all ;
% Create lattice
H = 24 ;
nNodes = 20 ;
lattice = Lattice.latticeEasy(H, nNodes, @rainDisturbanceARM) ;
%Run SDDP
params = sddpSettings('algo.McCount', 150, ...
'stop.iterationMax', 10,...
'stop.pereiraCoef', 2,...
'precise.computeEnd', true,...
'precise.count', 200,...
'solver','gurobi') ;
var.x = sddpVar(H) ;
var.p = sddpVar(H) ;
var.q = sddpVar(H) ;
var.l = sddpVar(H) ;
var.r = sddpVar(H) ;
lattice = compileLattice(lattice, @(scenario)nldsARM(scenario, var), params) ;
output = sddp(lattice, params) ;
plotOutput(output) ;
function [cntr, obj] = nldsARM(scenario, var)
t = scenario.getTime() ;
C = 25 ;
V = 1000 ;
D = 1000 ;
E = 1000;
PMax = 500;
InitialWater = 500;
InitialRain = 700;
x = var.x;
p = var.p;
q = var.q;
l = var.l;
r = var.r;
% Objective
obj = V * l(t) + C * p(t);
SupplyDemandBalance = p(t) + q(t) + l(t) >= D ;
Positivity = [l(t) >= 0, p(t) >= 0, x(t) >= 0, r(t) >= 0, q(t) >= 0];
%Posr = r(t) >= 0;
ThermalCapacityLimit = p(t) <= PMax ;
ReservoirLimit = x(t) <= E ;
if(scenario.time == 1)
HydroUtilizationLimit0 = q(t) <= InitialWater + r(t) ;
HydroLevelDynamics0 = x(t) <= InitialWater + r(t) - q(t) ;
RainDynamics0 = r(t) == max(InitialRain*scenario.data, 0);
cntr = [SupplyDemandBalance , Positivity, ...
ThermalCapacityLimit, ReservoirLimit, ...
HydroUtilizationLimit0, RainDynamics0, HydroLevelDynamics0];
else
HydroUtilizationLimit = q(t) <= x(t-1) + r(t) ;
HydroLevelDynamics = x(t) == x(t-1) + r(t) - q(t) ;
RainDynamics = r(t) == r(t-1)*scenario.data;
cntr = [SupplyDemandBalance , Positivity, ...
ThermalCapacityLimit, ReservoirLimit, ...
HydroUtilizationLimit, RainDynamics, HydroLevelDynamics];
end
end
function out = rainDisturbanceARM(t,i)
if t == 1
out = 1 + 0.1*noiseG(11);
else
out = 1 + 0.1*noiseG(i);
end
end
Optimal Management of Storage in Burkina Faso
The source codes can be downloaded from here and is based on the following publication: T. Kaneda, B. Losseau, A. Papavasiliou, L. Cambier, D. Scieur, P. Henneaux, N. Leemput, Optimal Management of Storage for Offsetting Solar Power Uncertainty using Multistage Stochastic Programming, 20th Power Systems Computation Conference, Dublin, Ireland, June 11-15, 2018.Data can be downloaded from load_1000_weekdays.mat and pv_1000_sample.mat. All come from this repository.
Problem statement
Africa has recently engaged in implementing an aggressive renewable energy integration plan. The major challenge is the management of excess energy, and introduction of large-scale battery storage has been considered as an attractive solution. Also we need to deal with the unpredictability of renewable energy supply and load. We here apply SDDP to this storage management problem under uncertainty. Suppose we have the following system: A large Photo-voltaic (PV) power plant is available, and its production is stochastic Load is stochastic, and might be different at each time step There exists a thermal generator with marginal cost \(MC = 200\) $/MWh and its lower/upper limit \(Pmin/Pmax\) is \(0/300\) MW Importation from other country \(pi\) is available with a cost \(CI = 100\) $/MWh and its capacity \(PI\) is 200 MW (but exportation is not allowed) There are five identical batteries with the charging/discharging rate \(BC/BD\) of 200 MW and their capacities \(ST\) are 1000 /MWh For the balancing, we are able to choose shedding production (\(ps\)) or load (\(ls\)) The value of lost load (\(VOLL\)) is 1000 \$/MWh Since the PV production and the load are both stochastic, we introduce net load (\(NL\)) which is defined as the difference between the load and the PV power. Net load is thus the stochastic parameter of the problem. Suppose we solve the problem on a 4-day horizon with hourly time step (\(H = 96\)), our objective is to minimize the expected cost over time.NLDS
We here present the \(NLDS_{t,k}\) at stage \(t\) and scenario \(k\) of the problem. The cost is caused by the thermal generator, import and load shedding. $$\min \left[ MC \cdot p_g + CI \cdot pi + VOLL \cdot ls \right]$$ where \(p_g\) is the production of generator. The power balance constraint can be written as follows: $$NL_{t,k} + \sum_{j = 1}^5 pd_j + ps = \sum_{j = 1}^5 pb_j + p_g + pi + ls$$ where \(pb_j\) is the power supply by discharging, \(pd_j\) is the charging power of battery \(j\). Net load \(NL_{t,k}\) can be different from one node of the lattice to another. The dynamics of the batteries is expressed as: $$s_j = s_{j,t-1} + \left(\eta_j \cdot pd_j - \frac{pb_j}{\mu_j}\right),\ j \in 1,...,5$$ where \(s_j\) represents the storage of battery \(j\) at the current stage \(t\), while \(s_{j,t-1}\) is a parameter which refers to the storage of battery \(j\) at the previous stage \(t-1\). \(\eta_j \ (<1)\) is the efficiency of charging and \(\mu_j \ (<1)\) is the efficiency of discharging battery \(j\). The following physical constraints are additionally introduced: $$\begin{aligned} &s_j \leq ST_j, pd_j \leq PD_j, pb_j \leq PB_j,\ j \in 1,...,5 \\ & pi \leq PI \\ &PMin_g \leq p_g \leq PMax_g \end{aligned} \label{eq:capacity}$$ Finally, we impose non-negativity constraints: $$ls,ps, pi, s_j,pb_j,pd_j,p_g \geq 0,\ g \in G,j \in 1,...,5 \label{eq:nonnegativity}$$Lattice
Let \(|\Omega|\) be the number of scenarios at each stage \(t = 2,...,H\). Since PV power and load can vary in time and day, at a given stage \(t\) we need to represent for each scenario the transition probability of going to every possible scenario at time \(t+1\), making a lattice of size \(|\Omega| \times |\Omega|\). We need to create a lattice of net load with the transition probability of \(|\Omega| \times |\Omega| \times H-1\) . This can be done using latticeEasyMarkovNonConst(H,P), a method included in Lattice of the toolbox. The argument P is a cell array of size \(H-1\), where \(P\{t\}(i, j)\) represents the probability of transition from realization \(i\) at time \(t\) to realization \(j\) at time \(t + 1\). In this case study, we set \(|\Omega| = 10\).Codes
close all;
clear all;
clc;
% Set parameters
n_days = 4; % horizon of the model (number of days)
n_nodes = 10; % number of nodes per a stage on the lattice
n_batteries = 5; % number of batteries
%% Import the data and build the lattice
% a) import the load data
load('load_1000_weekdays.mat');
% b) import the pv generation data
load('pv_1000_sample.mat');
% c) refine the data
[~,col]=find(pv);
start_time = min(col); % find the eariest sunrise time in the data
n_hours = size(pv,2); % number of hours to the horizon
data_size = size(pv,1);
pv = [pv(:,start_time:n_hours) ,pv(:,1:start_time-1)]; % redefine PV data
load_weekdays = [load_weekdays(:,start_time:n_hours),...
load_weekdays(:,1:start_time-1)]; % redefine load data
% c) build net load scenario
net_load = zeros(data_size,size(pv,2)*n_days);
for i= 1:n_days
perm1 = randperm(data_size);
perm2 = randperm(data_size);
for k = 1:data_size
net_load(k,(n_hours*(i-1)+1):n_hours*i) = ...
load_weekdays(perm1(k),:) - pv(perm2(k),:);
end
end
% d) build the lattice
H = size(net_load,2); % set horizon H
[transitionProba, data] = transprob(n_nodes,net_load);
% e) visualize the lattice
lattice = Lattice.latticeEasyMarkovNonConst(H, transitionProba, data);
figure ;
lattice.plotLattice(@(data) num2str(data));
%% Set the variables
var.ls = sddpVar(H); % load shedding
var.ps = sddpVar(H); % production shedding
var.s = sddpVar(n_batteries,H); % battery storage
var.pd = sddpVar(n_batteries,H); % pumping demand for batteries
var.pb = sddpVar(n_batteries,H); % supply of battery
var.pg = sddpVar(H); % thermal energy production
var.pi = sddpVar(H); % importing energy
var.cg = sddpVar(H); % cost of generator
params = sddpSettings('algo.McCount', 25,...
'stop.iterationMax',10,...
'stop.stopWhen','never',...
'verbose',1,...
'solver','gurobi');
lattice = compileLattice(lattice,@(scenario)nlds(scenario,var),params) ;
%% Solve the problem
output = sddp(lattice,params) ;
plotOutput(output);
function [ cntr, obj ] = nlds(scenario, var)
% Build NLDS for a given node
t = scenario.getTime() ;
% SETTING PARAMETERS AND VARIABLES
% parameters
VOLL = 1000; % cost of load shedding (in $/unit)
VOLP = 1e-5; % cost of production shedding (in $/unit)
COST_BT = 0.6e-5; % cost of charging battery ($/MWh)
MC = 200; % Marginal cost of diesel $/MWh
CI = 100; % Marginal cost of imports $/MWh
eta = 0.95; % fraction of converted energy to batteries
mu = 0.97; % fraction of converted energy from batteries
s_0 = 0; % level of energy in the batteries at the beginning
ST = 1000; % capacity of battery storage
PD = 200; % capacity of the pumping demand to store
PB = 200; % capacity of the battery power to extract
PG = 300; % total capacity of generators
PI = 200; % capacity of importation
% Variables
ls = var.ls; % load shedding
ps = var.ps; % production shedding
s = var.s; % battery storage
pd = var.pd; % pumping demand for batteries
pb = var.pb; % supply of battery
pg = var.pg; % thermal energy production
pi = var.pi; % importing energy
cg = var.cg; % cost of generator
data = scenario.data; % extract the stocastic parameter
net_load = data;
pd_total = sum(pd(:,t)); % total demand for batteries
pb_total = sum(pb(:,t)); % total supply of batteries
% OBJECTIVE FUNCTION
total_cost = ls(t)*VOLL + cg(t) + pi(t) * CI + ps(t) * VOLP + ...
(pd_total +pb_total)* COST_BT;
% CONSTRAINTS
% generator cost function
generator_cost = cg(t) >= pg(t) * MC;
% Power balance
power_balance = net_load + pd_total + ps(t) ==...
pb_total + pg(t) + pi(t) + ls(t) ;
% Storage balance in batteries
for i = 1:5
if t>1
storage_balance(i) = s(i,t) == ...
s(i,t-1) + (eta*pd(i,t) - 1/mu*pb(i,t));
else
storage_balance(i) = s(i,t) == ...
s_0 + (eta*pd(i,t) - 1/mu*pb(i,t));
end
end
% Maximum capacity
for i = 1:5
storage_capacity(i) = s(i,t) <= ST;
pumping_demand_capacity(i) = pd(i,t) <= PD;
extract_power_capacity(i) = pb(i,t) <= PB;
end
power_generated_capacity = pg(t) <= PG;
power_imported_capacity = pi(t) <= PI;
% Nonnegativity
for i = 1:5
s_non_negativity(i) = s(i,t) >= 0 ;
pd_non_negativity(i) = pd(i,t) >= 0 ;
pb_non_negativity(i) = pb(i,t) >= 0 ;
end
non_negativity = [s_non_negativity, pd_non_negativity, pb_non_negativity,...
ls(t) >= 0, pg(t) >= 0, pi(t) >= 0, cg(t) >= 0 , ps(t) >= 0 ] ;
% Inserting the model in the FAST format
obj = total_cost;
cntr = [generator_cost,...
power_balance, ...
storage_balance, ...
power_imported_capacity, ...
power_generated_capacity, ...
storage_capacity, ...
pumping_demand_capacity, ...
extract_power_capacity, ...
non_negativity
] ;
end
function [transitionProba_cell,nodeValue] = transprob(n_nodes,net_load)
% Generate transition probability (cell) and
% corresponding node values (matrix)
n_scenario = size(net_load,1);
n_hours = size(net_load,2);
transitionNode = zeros(n_nodes,n_nodes,n_hours-1);
transitionProb = zeros(n_nodes,n_nodes,n_hours-1);
transitionProba_cell = cell(n_hours-1,1);
% make ranges
ranges = zeros(n_nodes+1,n_hours); % ranges of net loads
nodeValue = zeros(n_hours,n_nodes); % node values
sort_nl = sort(net_load,1); % sorted net load
ranges(1,:) = sort_nl(1,:);
for i = 1:n_nodes
for t = 1:n_hours
ranges(i+1,t) = sort_nl(round(n_scenario/n_nodes*i),t);
if i == 1
nodeValue(t,i) = mean(sort_nl(1:round(n_scenario/n_nodes*i),t)); % take the average within the range
else
nodeValue(t,i) = mean(sort_nl(round(n_scenario/n_nodes*(i-1)):round(n_scenario/n_nodes*i),t));
end
end
end
% Building transition probability
for i=1:n_scenario
transisionCount = zeros(n_nodes,n_hours); % count the pv_values that in the range
for t=1:n_hours
for j=1:n_nodes
if j ~= 1
if ranges(j,t) < net_load(i,t) && net_load(i,t) <= ranges(j+1,t)
transisionCount(j,t) = transisionCount(j,t)+1;
end
else
if ranges(j,t) <= net_load(i,t) && net_load(i,t) <= ranges(j+1,t)
transisionCount(j,t) = transisionCount(j,t)+1;
end
end
end
end
for t = 2:n_hours
for l = 1:n_nodes
for q = 1:n_nodes
if transisionCount(l,t-1) == 1 && transisionCount(q,t) == 1
transitionNode(l,q,t-1) = transitionNode(l,q,t-1) + 1;
end
end
end
end
end
for t = 1:n_hours-1
for j = 1:n_nodes
for k = 1:n_nodes
if sum(transitionNode(j,:,t)) == 0
transitionProb(j,k,t) = 0;
else
transitionProb(j,k,t) = transitionNode(j,k,t)/sum(transitionNode(j,:,t));
end
end
end
end
% For stage 2 to H-1
for t = 2:n_hours-1
transitionProba_cell{t} = transitionProb(:,:,t);
end
% For stage 1: Choose the middile value and prob
nodeValue(1,1) = nodeValue(1,floor(n_nodes/2));
nodeValue(1,2:end) = zeros(1,n_nodes-1);
A = transitionProb(:,:,1);
transitionProba_cell{1} = A(floor(n_nodes/2),:);
end