clc
clear all
close all
warning off all

% definisce la variabile usata
syms z;

% funzione costo e valore (letterale)
cost_s = .2*z^2;
value_s = 20*(z/20)^(1/2);

% funzione costo e valore calcolata in un intervallo
int_x = 0.01;
x = [0:int_x:12];
cost = subs(cost_s,x);
value = subs(value_s,x);

% grafici
figure(1)
plot(x, cost)
hold on
plot(x, value, 'r')
legend('c', 'v', 'Interpreter', 'latex')
title('Funzioni costo e valore')

% c/v_punto_s = differenziali e funzioni calcolate in un intervallo
c_punto_s = diff(cost_s,z);
v_punto_s = diff(value_s,z);

c_punto = subs(c_punto_s,x);
v_punto = subs(v_punto_s,x);

% funzioni inverse delle derivate e calcolo in un intervallo
c_punto_inv_s = finverse(c_punto_s,z);
v_punto_inv_s = finverse(v_punto_s,z);

int_lambda = 0.01;
lambda=[0:int_lambda:300]';

c_punto_inv = subs(c_punto_inv_s,lambda);
v_punto_inv = subs(v_punto_inv_s,lambda);

figure(2)
plot(lambda, c_punto_inv)
hold on
plot(lambda, v_punto_inv, 'r')
title('$\dot{c}^{-1}$ e $\dot{v}^{-1}$', 'Interpreter', 'latex')
axis([0 300 0 500])


%% simulazione dell'algoritmo 

lambda_r0 = subs(v_punto_s, 10); % prezzo al dettaglio iniziale
durata = 100; %numero di passi

% vettori: prezzo al dettaglio, prezzo all'ingrosso, quantità consumata,
% ammanco o guadagno dell'ISO
prezzi_r = -1*ones(durata,1);
prezzi_w = -1*ones(durata,1);
quantita = -1*ones(durata,1); 
ISO_revenue = -1*ones(durata,1); 

% inizializzazioni
prezzi_r(1) = lambda_r0;
ISO_revenue(1) = 0;

% parametri
gamma = 0.5;    % coefficiente per la combinazione convessa
rho = 0.05;     % guadagno proporzionale per il controllo di ISO_revenue
sigma_2 = 0.01;    % varianza rumore sulla scelta del consumatore (media nulla)


for i=1:durata
    
    % dato il prezzo al dettaglio al tempo i
    % prezzi_r(i)
    
    % Passo 1: il consumatore decide, sulla base della sua value function,
    % la quantità di energia da consumare in [i, i+1)
    % L'eventuale aggiunta del rumore modella "indecisione" nel consumatore
    quantita(i) = subs(v_punto_inv_s, prezzi_r(i)) + sqrt(sigma_2)*randn;
    
    % Passo 2: alla fine dell'intervallo, la quantità di energia prodotta e
    % la quantità di energia consumata sono uguali (nessuna equazione)
    
    % Passo 3: l'ISO paga al produttore un prezzo per unità di energia che
    % massimizza il profitto (del produttore)
    prezzi_w(i) = subs(c_punto_s, quantita(i));
    
    % Passo 3.a: l'ISO registra la differenza fra il pagato e il ricevuto
    ISO_revenue(i+1) = quantita(i)*(prezzi_r(i) - prezzi_w(i));    
    
    % Passo 4: l'ISO, visto il prezzo (all'ingrosso) che ha dovuto pagare al produttore e
    % il prezzo (al dettaglio) che aveva proposto al consumatore all'inizio
    % dell'intervallo temporale, decide il prezzo al dettaglio per
    % l'intervallo successivo
    prezzi_r(i+1) = prezzi_r(i) + gamma*(prezzi_w(i) - prezzi_r(i)) - rho*ISO_revenue(i);
    
end   

%% La procedura seguente è quella implementata inizialmente. Se si usa
% questa bisogna ricordarsi di imporre una condizione iniziale al vettore
% quantita e non al vettore lambda_r (inoltre in questa non ci sono il rumore e il controllo sull'ISO)
% for i=1:durata
%     
%     % data la quantità di energia prodotta nell'intervallo [i-1, i]
%     % quantità(i)
%     
%     % Passo 1: l'ISO prevede quanta energia sarà utilizzata nell'intervallo
%     % successivo. La previsione è semplicemente
%     % s(t+1) = d(t)
%     
%     % Passo 2: l'ISO paga al produttore un prezzo per unità di energia che
%     % massimizza il profitto (del produttore)
%     prezzi_w(i) = subs(c_punto_s, quantita(i));
%     
%     % Passo 3: l'ISO, visto il prezzo (all'ingrosso) che ha dovuto pagare al produttore e
%     % il prezzo (al dettaglio) che aveva proposto al consumatore all'inizio
%     % dell'intervallo temporale, decide il prezzo al dettaglio per
%     % l'intervallo successivo
%     prezzi_r(i+1) = prezzi_r(i) + gamma*(prezzi_w(i) - prezzi_r(i));
%     
%     % Passo 1: il consumatore decide, sulla base della sua value function,
%     % la quantità di energia da consumare in [i, i+1)
%     quantita(i+1) = subs(v_punto_inv_s, prezzi_r(i+1));
%    
% end   


tempi=[1:durata]';

figure(3)
plot(tempi,prezzi_r(1:durata))
hold on
plot(tempi,prezzi_w, 'r')
legend('\lambda_r', '\lambda_w', 'Interpreter', 'latex')
title('Prezzi', 'Interpreter', 'latex')

figure(4)
plot(tempi,prezzi_w-prezzi_r(1:durata))
title('Differenza prezzi')

figure(5)
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 12 0 5])

figure(6)
plot(tempi, ISO_revenue(1:durata))
title('Ammanco/guadagno dell''ISO')

