clear all
close all
clc
plot_settings

% Time intervals
T = 300;
t = 0:T;

% Heat demand and supply
mean = 200;
amplitude = 50;
Dmin = 100;
Dmax = 300;

D = mean + amplitude * sin(2*2*pi*t/T -pi/2) + normrnd(0,20,[1,T+1]);
D = min(max(D,Dmin),Dmax);

E = ones(1,length(t))*210;

figure
plot(t,D,t,E)
axis([0 T 0 400])
xlabel('Time period, $t$')
ylabel('Energy [MWh]')
legend('Heat demand, $D_t$', 'Heat supply, $E_t$')
%title('Demand model: Sinusoidal')

% Electricity price
Pmin = 0;
Pmax = 2500;

P0 = 200;
P1 = normrnd(0,50,1,T+1);
P2 = normrnd(0,500,1,T+1);

a =  unifrnd(0,1,1,T+1);
p = 0.031;
for i = 1:length(a)
    if a(i)>p
        a(i) = 0;
    else
        a(i) = 1;
    end
end

P = P0 + P1 + a.*P2;
P = min(max(P,Pmin),Pmax);

figure
plot(t,P)
xlabel('Time period, $t$')
ylabel('Price, $P_t$ [NOK/MWh]')
%title('Price model: First-order Markov chain + jump')

%%
HH = 11;
N=20;
tt = -10:N-1;

Dp = mean + amplitude * sin(2*2*pi*t/T -pi/2);
Dp = min(max(Dp,Dmin),Dmax);


P00 = P(HH);
P100 = normrnd(0,50,1,N+1);
Pp = P00 + P100;
Pp = min(max(Pp,Pmin),Pmax);

Dhp = [D(1:HH) Dp(HH+1:N+HH-1)]
Php = [P(1:HH) Pp(1:N-1)]


figure()
subplot(1,2,1)
plot(tt,Dhp)
for i = HH:size(Dhp,2)
hold on
plot(tt(i),Dhp(i),'o','Color','r','MarkerSize',5)
end
axis([-10 N 0 300])
xlabel('Time period, $t$')
ylabel('Energy Demand [MWh]')
title('(a)')


subplot(1,2,2)
plot(tt,Php)
for i = HH:size(Dhp,2)
hold on
plot(tt(i),Php(i),'o','Color','r','MarkerSize',5)
end

axis([-10 N 0 2000])
xlabel('Time period, $t$')
ylabel('Price, $P_t$ [NOK/MWh]')
title('(b)')


%%

Dhp2 = [D(1:HH) Dp(HH+1:N+HH-1)*1.5]

figure()
plot(tt,Dhp,'color','r')
for i = HH:size(Dhp,2)
hold on
plot(tt(i),Dhp(i),'x','Color','r','MarkerSize',5)
end

hold on
plot(tt,Dhp2,'color','b')
for i = HH:size(Dhp,2)
hold on
plot(tt(i),Dhp2(i),'o','Color','b','MarkerSize',5)
end

axis([-10 N 0 300])
xlabel('Time period, $t$')
ylabel('Energy Demand [MWh]')