clear all
close all
clc
tic

%% Load
addpath('../Exogenous_information')
load('sample_path.mat')
plot_settings

%% setting

% number of scenario & time horizon
n_scenario = size(sample_path.D,1);
time = size(sample_path.D,2);

% Initial state
S0 = S0();
S = S0.S;
rate_d = S0.rate_d*ones(n_scenario,1);
rate_c = S0.rate_c*ones(n_scenario,1);
R_max = S0.R_max;

    %initial states
    R = S.R*ones(n_scenario,1);
    D = S.D*ones(n_scenario,1);
    E = S.E*ones(n_scenario,1);
    P = S.P*ones(n_scenario,1);

%Price_parameters
P_high = 190;
P_low = 120;

%% PFA policy based on specific parameters for low and high prices.
%data storage platform
uSim.WD = [];
uSim.GD = [];
uSim.RD = [];
uSim.WR = [];
uSim.GR = [];
xSim.R = [];
xSim.D = [];
xSim.E = [];
xSim.P = [];
xSim.C = [];
for t = 1:time
    
    %policy => control variable
    WD = min(D,E);
    
    comp_h = P>P_high;
    RD = comp_h.*min(min(D-WD,rate_d),R)*S0.eff_d;
    
    GD = D-WD-RD;
    
    WR = min(min(E-WD,rate_c),R_max-R);
    
    comp_l = P<P_low;
    GR = comp_l.*min(rate_c-WR,R_max-R-WR);
    
    
    C = P.*(GD+GR);
    
    %Data storage
    uSim.WD = [uSim.WD, WD];
    uSim.GD = [uSim.GD, GD];
    uSim.RD = [uSim.RD, RD];
    uSim.WR = [uSim.WR, WR];
    uSim.GR = [uSim.GR, GR];
    xSim.R = [xSim.R, R];
    xSim.D = [xSim.D, D];
    xSim.E = [xSim.E, E];
    xSim.P = [xSim.P, P];
    xSim.C = [xSim.C, C];

    %update states
    R = R+S0.eff_c*(WR+GR)-RD;
    D = sample_path.D(1:n_scenario,t);
    E = sample_path.E(1,t)*ones(n_scenario,1);
    P = sample_path.P(1:n_scenario,t);
end
% Cost
Cost = sum(xSim.P.*(uSim.GD+uSim.GR),2);
for k = 1:n_scenario
    Avg_Cost_cum(k) = mean(Cost(1:k));
end
Avg_Cost = mean(Cost)

save('P01_simple')

%% Plot
Line = 1.1;
Font = 14;
N=300;
time = linspace(1,301,301);
figure()
subplot(5,1,1)
stairs(time,uSim.WD(n_scenario,:),'LineWidth',1.1)
%legend('$Q_{WD}$')
axis([0 N 0 300])
%xlabel('Time period, $t$')
ylabel('$Q^{WD}$ [MWh]','FontSize',Font)
set(gca,'FontSize',Font)

subplot(5,1,2)
stairs(time,uSim.GD(n_scenario,:),'LineWidth',1.1)
%legend('$Q_{GD}$')
axis([0 N 0 300])
%xlabel('Time period, $t$')
ylabel('$Q^{GD}$ [MWh]','FontSize',Font)
set(gca,'FontSize',Font)

subplot(5,1,3)
stairs(time,uSim.RD(n_scenario,:),'LineWidth',1.1)
%legend('$Q_{GD}$')
axis([0 N 0 300])
%xlabel('Time period, $t$')
ylabel('$Q^{SD}$ [MWh]','FontSize',Font)
set(gca,'FontSize',Font)

subplot(5,1,4)
stairs(time,uSim.WR(n_scenario,:),'LineWidth',1.1)
%legend('$Q_{WS}$')
axis([0 N 0 300])
%xlabel('Time period, $t$')
ylabel('$Q^{WS}$ [MWh]','FontSize',Font)
set(gca,'FontSize',Font)

subplot(5,1,5)
stairs(time,uSim.GR(n_scenario,:),'LineWidth',1.1)
%legend('$Q_{GS}$')
axis([0 N 0 300])
xlabel('Time period, $t$','FontSize',Font)
ylabel('$Q^{GS}$ [MWh]','FontSize',Font)
set(gca,'FontSize',Font)

Font = 30;
figure()
%subplot(2,1,1)
stairs(time,xSim.R(n_scenario,:),'LineWidth',1.5)
%legend('$Q_{GS}$')
axis([0 N 0 S0.R_max])
%xlabel('Time period, $t$')
ylabel('Energy amount in Storage, $S_t$ [MWh]','FontSize',Font)
xlabel('Time period, $t$','FontSize',Font)
set(gca,'FontSize',Font)

