clc;
clear all;
close all;
warning off all;

int_x = 0.01;
x=[0:int_x:15]';
syms z;

cost_s =.5*z^(2); %funzione costo (letterale)
value_s_1 =20*log(z*5+1); %funzione valore (letterale)
value_s_2 =30*log(z*5+1); %funzione valore (letterale)
value_s_3 =25*log(z*5+1); %funzione valore (letterale)
value_s_4 =15*log(z*5+1); %funzione valore (letterale)

cost=subs(cost_s,x); %funzione costo (valori)
value_1=subs(value_s_1,x);%funzione valore (valori)
value_2=subs(value_s_2,x);%funzione valore (valori)
value_3=subs(value_s_3,x);%funzione valore (valori)
value_4=subs(value_s_4,x);%funzione valore (valori)

%% Calcolo delle derivate delle funzioni costo e value
% c/v_punto_s = differenziali

c_punto_s=diff(cost_s,z);
v_punto_s_1=diff(value_s_1,z);
v_punto_s_2=diff(value_s_2,z);
v_punto_s_3=diff(value_s_3,z);
v_punto_s_4=diff(value_s_4,z);

c_punto=subs(c_punto_s,x); % lambda_w
v_punto_1=subs(v_punto_s_1,x); % lambda_r
v_punto_2=subs(v_punto_s_2,x); % lambda_r
v_punto_3=subs(v_punto_s_3,x); % lambda_r
v_punto_4=subs(v_punto_s_4,x); % lambda_r


%% Calcola le inverse della c punto e v punto
c_punto_inv_s =finverse(c_punto_s,z);
v_punto_inv_s_1 =finverse(v_punto_s_1,z);
v_punto_inv_s_2 =finverse(v_punto_s_2,z);
v_punto_inv_s_3 =finverse(v_punto_s_3,z);
v_punto_inv_s_4 =finverse(v_punto_s_4,z);

int_lambda = 0.001; %intervallo dei valori dei prezzi
lambda=[0:int_lambda:300]';

% calcola inverse 
c_punto_inv =subs(c_punto_inv_s,lambda);
v_punto_inv_1 =subs(v_punto_inv_s_1,lambda);
v_punto_inv_2 =subs(v_punto_inv_s_2,lambda);
v_punto_inv_3 =subs(v_punto_inv_s_3,lambda);
v_punto_inv_4 =subs(v_punto_inv_s_4,lambda);

%%

x0=4.4; %quantita iniziale di energia
durata=24; %numero di passi
prezzi_r = -1*ones(durata,1); %vettore contenente i prezzi al dettaglio
prezzi_w = -1*ones(durata,1); %vettore contenente i prezzi all'ingrosso
xt=x0;
quantita = -1*ones(durata,1); %vettore delle quantita di energia
gamma = 0.5;

% somme delle differenze dei prezzi al dettaglio e all'ingrosso
integrale = zeros(durata,1); %vettore contenente i "guadagni" dell'ISO
temp = 0;

for i=1:6
    lambda_r = subs(v_punto_s_1, xt);
    lambda_w = subs(c_punto_s, xt);
    lambda_ott = gamma*lambda_w + (1 - gamma)*lambda_r;
    
    prezzi_r(i,1)=lambda_r;
    prezzi_w(i,1)=lambda_w;
    quantita(i,1)=xt;
    
    %lambda_ott=round(lambda_ott/int_lambda+1);    
    xt= subs(v_punto_inv_s_1, lambda_ott);
    
    integrale(i,1) = (prezzi_r(i,1)-prezzi_w(i,1))*xt + temp;
    temp = integrale(i,1);
end
for i=7:12
    lambda_r = subs(v_punto_s_2, xt);
    lambda_w = subs(c_punto_s, xt);
    lambda_ott = gamma*lambda_w + (1 - gamma)*lambda_r;
    
    prezzi_r(i,1)=lambda_r;
    prezzi_w(i,1)=lambda_w;
    quantita(i,1)=xt;
    
    %lambda_ott=round(lambda_ott/int_lambda+1);    
    xt= subs(v_punto_inv_s_2, lambda_ott);
    
    integrale(i,1) = (prezzi_r(i,1)-prezzi_w(i,1))*xt + temp;
    temp = integrale(i,1);
end
for i=13:18
    lambda_r = subs(v_punto_s_3, xt);
    lambda_w = subs(c_punto_s, xt);
    lambda_ott = gamma*lambda_w + (1 - gamma)*lambda_r;
    
    prezzi_r(i,1)=lambda_r;
    prezzi_w(i,1)=lambda_w;
    quantita(i,1)=xt;
    
    %lambda_ott=round(lambda_ott/int_lambda+1);    
    xt= subs(v_punto_inv_s_3, lambda_ott);
    
    integrale(i,1) = (prezzi_r(i,1)-prezzi_w(i,1))*xt + temp;
    temp = integrale(i,1);
end
for i=19:24
    lambda_r = subs(v_punto_s_4, xt);
    lambda_w = subs(c_punto_s, xt);
    lambda_ott = gamma*lambda_w + (1 - gamma)*lambda_r;
    
    prezzi_r(i,1)=lambda_r;
    prezzi_w(i,1)=lambda_w;
    quantita(i,1)=xt;
    
    %lambda_ott=round(lambda_ott/int_lambda+1);    
    xt= subs(v_punto_inv_s_4, lambda_ott);
    
    integrale(i,1) = (prezzi_r(i,1)-prezzi_w(i,1))*xt + temp;
    temp = integrale(i,1);
end

%%

t=[1:durata]';
% 
% figure(4);
% plot(t,prezzi_r);
% title('$\lambda_r$', 'Interpreter', 'latex');
% 
% figure(5);
% plot(t,prezzi_w);
% title('$\lambda_w$', 'Interpreter', 'latex');
% 
% figure(6);
% plot(t,prezzi_r-prezzi_w);
% title('Differenza dei prezzi al dettaglio e all''ingrosso');
% 
% figure(7);
% plot(t,integrale);
% title('Somme delle differenze di prezzo');

integrale(durata,1)

figure(8);
plot(x,c_punto, 'k');
hold on;
plot(x,v_punto_1,':m');
hold on;
plot(x,v_punto_2,':c');
hold on;
plot(x,v_punto_3,':b');
hold on;
plot(x,v_punto_4,':k');
for t=1:durata-1
    plot([quantita(t) quantita(t)],[prezzi_r(t) prezzi_r(t+1)],'r')
    plot([quantita(t) quantita(t+1)],[prezzi_r(t+1) prezzi_r(t+1)],'r')
    plot([quantita(t) quantita(t)], [prezzi_r(t+1) subs(c_punto_s,quantita(t))],'k')
end    
set(gca, 'FontSize', 16, 'FontName', 'Times New Roman')
xlabel('Energy unit')
ylabel('\lambda^r, \lambda^w')
title('Energy vs price -- Trend', 'Interpreter', 'latex');
axis([2 8 2 8]);

figure(12);
plot(x,c_punto,'b',quantita,prezzi_r,'o',x,v_punto_1,'r',x,v_punto_2,'r',x,v_punto_3,'r',x,v_punto_4,'r');
set(gca, 'FontSize', 16, 'FontName', 'Times New Roman')
xlabel('Energy unit')
ylabel('\lambda^r, \lambda^w')
title('Energy vs price -- Settling', 'Interpreter', 'latex');
axis([2 8 2 8]);

save file_dati quantita prezzi_r;

