close all;
clear all;

% Nominal energy capacity in kWh
E_bar = 100;

% Sampling time in hours (ie, you have one sample each Delta_T)
Delta_T = 600/3600;

% Initial SOC
SOC_0 = 0.5;

% Lenght of the "simulation" in number of steps.
T = 100;

% Generate random power profile with a standard deviation of 10 kW per
% sample.
rng(2);
P = randn(T, 1)*10;



% Efficiency (assumed equal for charging and discharging)
eta = 0.9;



% In the following, I present a few options on how to code the exercise.


% -- 1. Some explicit way.

% Create a vector (or array) where to store the result of the SOC
% computation. Its first element is the initial SOC.
SOC = nan(T, 1);
SOC(1) = SOC_0;

for t=1:T-1
    if P(t) >= 0
        SOC(t+1) = SOC(t) + eta * P(t)/E_bar * Delta_T;
    else
        SOC(t+1) = SOC(t) + 1/eta * P(t)/E_bar * Delta_T;
    end
end

% Vector time for plotting
time = ((1:T) - 1) * Delta_T;

figure
plot(time, SOC)
hold on

plot(time, time * 0 + 0, 'r--')
plot(time, time * 0 + 1, 'r--')

hold off

xlabel('Time [hours]');
ylabel('SOC (per unit)');
legend('SOC estimation', 'Range');
ylim([-0.1, 1.1]);

title('Method 1')

% Note that the physical boundaries are plotted just for reference. 
% We don't no any effort in enforcing them so far, namely the P profile
% (random), could drive the SOC to unfeasible levels; we will see this
% later.


% -- 2. We avoid the if by introducing the negative and positive part
% operators

% Positive and negative part [x]+ and [x]- operators code as inline functions.
pos_part = @(x) max(x, 0);
neg_part = @(x) -min(x, 0);


SOC = nan(T, 1);
SOC(1) = SOC_0;
for t=1:T-1
    SOC(t+1) = SOC(t) + (eta * pos_part(P(t)) - 1/eta * neg_part(P(t))) * Delta_T/E_bar;    
end


figure
plot(time, SOC)
hold on

plot(time, time * 0 + 0, 'r--')
plot(time, time * 0 + 1, 'r--')

hold off

xlabel('Time [hours]');
ylabel('SOC (per unit)');
legend('SOC estimation', 'Range');
ylim([-0.1, 1.1]);

title('Method 2')



% -- 3. One liner using vector notation

SOC = SOC_0 + Delta_T/E_bar * cumsum(eta * pos_part(P) - 1/eta * neg_part(P));


% Cumsum (cumulative sum) is a linear operation and could be rendered as a
% matrix product. We will see this when formulating optimization problems.

figure
plot(time, SOC)
hold on

plot(time, time * 0 + 0, 'r--')
plot(time, time * 0 + 1, 'r--')

hold off

xlabel('Time [hours]');
ylabel('SOC (per unit)');
legend('SOC estimation', 'Range');
ylim([-0.1, 1.1]);

title('Method 3')