clc
clear all
close all
warning off all

% definisce la variabile usata
syms z;

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

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

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

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


%% 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.6;    % coefficiente per la combinazione convessa
rho = 0.03;     % guadagno proporzionale per il controllo di ISO_revenue
sigma_2 = 0.05;    % varianza rumore sulla scelta del consumatore (media nulla)

Sigma = zeros(2,2);
var_camp_lambda = 0;
var_camp_I = 0;

numero_prove = 1;

for num_prove = 1:numero_prove
    
    num_prove
    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) = ISO_revenue(i) + 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

    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', '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 12 0 12]);

    figure(6)
    plot(tempi, ISO_revenue(1:durata), 'k', 'LineWidth', 2)
    grid on
    set(gca, 'FontSize', 16, 'FontName', 'Times New Roman')
    xlabel('Time')
    ylabel('ISO revenue')
    title('ISO revenue')

    
    temp_r = prezzi_r(50:end, 1) - sqrt(20);
    temp_I = ISO_revenue(50:end,1);
        
    var_camp_lambda = var_camp_lambda + var(temp_r);
    var_camp_I = var_camp_I + var(temp_I);
     
end

