EMS-SMPC

This repository contains the scripts 📜 and data 📂 used for the EMS-SMPC based on this paper.

RANDOM FORESTS (RF) 🌳

Need to load data:

• GFA_15_min   (platform load)
• RES               (wind speed)

load('\\home.ansatt.ntnu.no\spyridoc\Documents\MATLAB\J2_PAPER\EMS-SMPC\DataFiles\dataTimeseries.mat')

1. TRAINING

   • train_RF.m   train RF models
   • Output:
     • Mdl_ld (load RF models)
     • Mdl_wp (wind power RF models)

2. OUTPUT FIGURES

   • figures_section_2_2.m    create figures of paper Section 2.2 (MSE, NRMSE, zoom forecasts k=1/k=6) (Figures 2-3)
   • forecast_plots.m            plot deterministic, probabilistic and scenario K steps ahead forecasts for selected date

MAIN SIMULATIONS (DMPC/SMPC) 📔

Need to load trained RF models:

load('\\home.ansatt.ntnu.no\spyridoc\Documents\MATLAB\J2_PAPER\EMS-SMPC\DataFiles\trainedRFmodels.mat')

1. INPUTS 📋 💾

   • user_defined_inputs.m   define the input variable

2. SYSTEM SIMULATION

   • main.m                  simulate system with DMPC/SMPC
   • optProb.m             formulate the MILP optimization for poblem for DMPC/SMPC
   • control_plots.m   plot MPC states evolution and control effort for selected method (DMPC/SMPC)

3. PAPER RESULTS

   • kpi_compare.m                     perform comparison of calculated KPIs from DMPC/SMPC (Table 8)
   • figures_section_4_1_2.m    create figures of Section 4.1.2 (forecasts on events, mean, quantiles, scenarios) (Figures 7-9)
   • figures_section_4_2.m       create figures of Section 4.2 (net load, rule-based areas, GT status and SoC (DMPC/SMPC)) (Figures 10-15)

Default input values 🔢

• input.startingDay = 100
• input.durationDays = 1
• input.giveStartingTime = 0 % {0, 1}
• inut.startingTime = 7630
• input.doAnimation = 0 % {0, 1}
• input.animationVar = 'load' % {'load', 'wind'}
• input.randomSeed = 24
• input.method = 'scn_frcst' % {'point_frcst', 'scn_frcst'}
• input.degradWeight = 'noWeight' % {'noWeight','none', 'normal', 'low', 'medium', 'high'}
• input.N_steps = 300
• input.N_prd = 6 % {MPC simulation, CRPS calculation} = {6, 12}
• input.lgdLocationDstrb = 'southwest'
• input.lgdLocationIgtOn = 'southeast'
• input.lgdLocationSoC = 'southeast'

Specialized functions used by the above main scripts ⬇️

1. PARAMETERS 📊

   preamble.m assigns 📎 parameter values to variable par.

   Default par values 🔢

   
   

   Geenric

   • par.Ts = 15 % Timestep (minutes)
   • par.dol2eur = 0.89 % dollars to euros conversion
   • par.rhoGas = 0.717 % Natural Gas density [kg/m^3]

   
   

   Sets

   • par.N_pwl = 11 % # of discretization points for PieceWise Linear approx.
   • par.N_gt = 4 % # of Gas Turbines
   • par.N_scn = 10 % # of scenarios

   
   

   Random forests

   • par.leafSizeIdx = 1
   • par.lamda = 0.5
   • par.tau = linspace(0,1,21)
   • par.lagsNum = 6

   
   

   Cost coeeficicents

   • par.c_dump = 10*100 % artificial cost (per unit of dumped power per period)
   • par.c_soc_dev = 0*10*100*100 % artificial cost (per unit of absolute SoC deviation in the end)
   • par.c_fuel = 0.24/par.rhoGas * par.dol2eur % [euros/kgGas]
   • par.c_gt_srt = 1217 % [euros/GTstart]
   • par.c_gt_ON = 5000 % [euros/GT_ON sattus]
   • par.c_Bat_rpl = par.degradWeight * 500000 * par.dol2eur % replacement cost [euros/MWh]
   • par.c_Bat_res = par.degradWeight * 50000 * par.dol2eur % residual value [euros/MWh]

   
   

   Gas Turbines

   • par.P_gt_nom = 20.2 % Nominal GT power rating
   • par.P_gt_min = 0.20 * par.P_gt_nom
   • par.P_gt_max = 1.09 * par.P_gt_nom
   • par.gt_RR = par.P_gt_max % Ramping Rate
   • par.spinRes = 1.05
   • par.idleFuel = 172*0.2*20.2+984 % [kg/h] coming from min GT fuel consumption (linear curve) - intercept @ no load: 1720.220.2+984
   • par.P_gt_data = linspace(par.P_gt_min,par.P_gt_max,par.N_pwl)
   • par.fuel_data = (0.5109 * par.P_gt_data.^2 -20.933 .* par.P_gt_data + 433.83) % [kg/MWh]

   
   

   BESS

   • par.eta_ch = 0.95 % charging efficieny
   • par.eta_dis = 0.95 % discharging efficieny
   • par.P_bat_max = 5 % power rating [MW] nominal: 5
   • par.E_bat_max = 10 % capacity rating [MWh] nominal: 10
   • par.socUPlim = 0.8 % up SoC limit [-]
   • par.socDOWNlim = 0.2 % down SoC limit [-]
   • par.SoC_ref = 0.5 % reference SoC

   
   

   Degradation

   • par.batLifetime = 10 % Lifetime expectancy of battery
   • par.a = 1591.1 % Proportional constant of cycling curve
   • par.k = 2.089 % Exponent of cycling curve
   • ar.DoD_data = linspace(0,1,par.N_pwl)' % Depth-Of-Discharge [0-1]
   • par.Ncyc = par.a*par.DoD_data.^(-par.k) % Cycle lifetime (# of cycles)
   • par.rho_data = 100*100./par.Ncyc % Percentage degradation [%] - (times 100 for scaling purposes)

2. SCENARIO FORECASTING 🔮

• funScenGenQRF1.m             issue K steps ahead scenarios forecasts at time t
• funScenGenQRF.m               issue K steps ahead scenarios forecasts itteratively for t>1 (used with funAnimateQRF.m)
• funScenFrcstFig1step.m   plot K steps ahead scenario forecasts for selected time t
• funProbFrcstFig1step.m   plot K steps ahead quantile forecasts for selected time t
• funFrcstFig1step.m          plot both scenario and quantile K steps ahead forecasts for selected time t
• funAnimateQRF.m               produce forecasting animation (.gif) for itterative forecasts (t>1)

3. GET/PLOT RESULTS 💡

• funGetCtrlRslt.m             calculated KPIs from simulation in variable kpi
• funPltCtrlRslt.m             plot the control variables (states and inputs) results for a method (DMPC/SMPC)
• funPltCtrlRsltPretty.m   plot the control variables (states and inputs) results for a method (DMPC/SMPC) in a pretty way

4. CRPS 📈

• funCalcCRPS.m            calculate CRPS, skill score and plot CRPS for all lead times, for selected period of time with QRF forecasts and benchmark method (paper: Figure 6, Table 1)
• funCalcCRPS1step.m   calculate CRPS for single predictions and plot cdf to compare QRF and benchmark method forecasts for each lead time
• funCovCorrGenQRF.m   Estimate covariance and correlation matrices from inverse transformed data based on probabilistic forecasts from QRF models

SIMULATION CASE STUDIES ⌛

Cases studies (Section 4.2) 📅

input.durationDays = 1 and input.giveStartingTime = 0

• input.startingDay=100 (10 April)
• input.startingDay=118 (27 April)
• input.startingDay=226 (14 August)
• input.startingDay=61 (02 March)
• input.startingDay=166 (15 June)
• input.startingDay=160 (09 June)

Irregular events (Section 4.1.2) 📅

input.durationDays = 0 and input.giveStartingTime = 1 and par.N_scn = 25 (for scenarios visualization)

Load

• inut.startingTime= 7630
• inut.startingTime= 7635
• inut.startingTime= 7636
• inut.startingTime= 7709

Wind

• inut.startingTime= 7646

• inut.startingTime= 7647

• inut.startingTime= 7648

• inut.startingTime= 7649

• inut.startingTime= 7760

• inut.startingTime= 7761

• inut.startingTime= 7762

• inut.startingTime= 7763

DATA FILES 📂

• To upload the data