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

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

cost_s =0.5*z^(2); %funzione costo (letterale)
value_s = 20*log(z*5); %funzione valore (letterale)

cost=subs(cost_s,x); %funzione costo (valori)
value=subs(value_s,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=diff(value_s,z);

c_punto=subs(c_punto_s,x); % lambda_w
v_punto=subs(v_punto_s,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 =finverse(v_punto_s,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 =subs(v_punto_inv_s,lambda);



%%
x0=10; %quantita iniziale
durata=100; %numero di passi
prezzi_r=-1*ones(durata,1);
prezzi_w=-1*ones(durata,1);
xt=x0;
quantita=-1*ones(durata,1); %vettore delle quantita di energia
gamma = 0.2;
rho = 0.01;

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

for i=1:durata
    lambda_w = subs(c_punto_s, xt);
    lambda_r = subs(v_punto_s, xt);
    lambda_ott = lambda_r + gamma*(lambda_w - lambda_r) - rho*temp;
    
    prezzi_r(i,1)=lambda_r;
    prezzi_w(i,1)=lambda_w;
    quantita(i,1)=xt;
    xt = subs(v_punto_inv_s, lambda_ott);

    y(i) = temp + (prezzi_r(i,1)-prezzi_w(i,1))*xt;
    temp = y(i);
end   

%%
t=[1:durata]';

% figure(1);
% plot(x,cost,'b',x,value,'r');
% title('Funzione costo (blu), funzione valore (rosso)');
% 
% figure(2);
% plot(x,c_punto,'b',x,v_punto,'r');
% title('Derivate costo e valore');

% 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');

figure(7);
plot(t,y, 'k', 'LineWidth', 2);
grid on
set(gca, 'FontSize', 16, 'FontName', 'Times New Roman')
xlabel('Time')
ylabel('Money')
title('Total revenue of the ISO', 'FontSize', 16, 'FontName', 'Times New Roman', 'Interpreter', 'latex');



figure(8);
plot(x,c_punto, 'k');
hold on;
plot(x,v_punto,':k', 'LineWidth', 2);
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', 'Interpreter', 'latex');
axis([0 15 0 15]);


figure(9);
plot(x,c_punto, 'k');
hold on;
plot(x,v_punto,':k', 'LineWidth', 2);
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 - Zoom', 'Interpreter', 'latex');
axis([3 5.5 3 6]);



