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+1); %funzione valore (letterale)

cost=subs(cost_s,x); %funzione costo (valori)
value=subs(value_s,x);%funzione valore (valori)

figure(1);
plot(x,cost,'r',x,value,'b');
title('Funzione costo (rosso), funzione valore (blu)');

%% 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

figure(2);
plot(x,c_punto,'r',x,v_punto,'b');
title('Derivate costo e valore');
axis([0 15 0 25]);

%%

c_punto2_s =diff(c_punto_s,z);
c_punto2 =subs(c_punto2_s,x); % derivata seconda di c(x) = c
v_punto2_s =diff(v_punto_s,z);
v_punto2 =subs(v_punto2_s,x); % derivata seconda di v(x) = derivata di a(x)

figure(3);
plot(x,c_punto2,'r',x,v_punto2,'b');
title('Derivate seconde costo e valore');
axis([0 15 -40 25]);

%%
funz_s = 2*v_punto2_s/(v_punto2_s - c_punto2_s);
funz = subs(funz_s,x);
gamma_lim = subs(funz_s,4.3733);

figure(4);
plot(x,funz);
title('Funzione 2*a''(x)/[a''(x) - c]');
axis([0 15 0 4]);

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

% figure(3);
% plot(lambda,c_punto_inv,lambda,v_punto_inv);
% title('$\dot{c}^{-1}$','$\dot{v}^{-1}$', 'Interpreter', 'latex');

%%
x0=10; %quantita iniziale di energia
durata=100; %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;
gamma = gamma_lim - 0.02*gamma_lim;

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

for i=1:durata
    lambda_w = subs(c_punto_s, xt);
    lambda_r = subs(v_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;
    
   if(i==1)
        integrale(1,1) = (prezzi_r(1,1)-prezzi_w(1,1))*xt;
   end
   if (i>1)
        temp1 = (prezzi_r(i,1)-prezzi_w(i,1))*xt;
        integrale(i,1) = integrale(i-1,1)+temp1;
   end
   
   xt = subs(v_punto_inv_s, lambda_ott);
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,':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    
title('$\dot{c},\dot{v}$ e andamento delle quantita', 'Interpreter', 'latex');
axis([0 15 0 15]);



