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

int_x = 0.01;
x=[0:int_x:12]';
cas_1=2; %1*rand(1,1)+1
cas_2=2; %1*rand(1,1)+1;
syms z;

cost_s =.5*z^(cas_1); %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);
%title('z^{cas_1}');
title('Funzione costo');

%% 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,x,v_punto,'r');
title('Derivate costo e valore');

%% 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);
title('$\dot{c}^{-1}$', 'Interpreter', 'latex');

%%
x0=10; %quantita iniziale di energia
durata=30; %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
[xottimo,lambdaottimo]= polyxpoly(x,v_punto,x,c_punto);             % punto di ottimo
gamma = 0.4;
flag=false;
regime=0; % quanto ci mette a portarsi circa al prezzo ottimo
y = zeros(durata,1); %vettore contenente i "guadagni" dell'ISO
temp = 0;


% primo passo
lambda_r = subs(v_punto_s, xt);
lambda_w = subs(c_punto_s, xt);

prezzi_r(1,1)=lambda_r;
prezzi_w(1,1)=lambda_w;

lambda_ott = gamma*lambda_w + (1 - gamma)*lambda_r;
quantita(1,1)=xt;
y(1) = temp + (prezzi_r(1,1)-prezzi_w(1,1))*xt;
temp = y(1);
xt = subs(v_punto_inv_s, lambda_ott);


for i=2:durata
    quantita(i,1)=xt;
    lambda_r = subs(v_punto_s, xt);
    lambda_w = subs(c_punto_s, xt);

    prezzi_r(i,1)=lambda_r;
    prezzi_w(i,1)=lambda_w; 
    
    x_temp = quantita(1:i);
    y_temp = prezzi_r(1:i);
    
    if(i==2) % approssimo con una retta
        p = polyfit(x_temp,y_temp,1);
        xi=linspace(0,12,1/int_x+1);
        yi=polyval(p,xi);
        figure(12)
        plot(x,c_punto,xi,yi,x_temp,y_temp,'o');
    end  
    if (i>2)  % approssimo con una parabola
        p = polyfit(x_temp,y_temp,2);
        xi=linspace(0,12,1/int_x+1);
        yi=polyval(p,xi);
        figure(12) %guardo come evolve
        plot(x,c_punto,xi,yi,x_temp,y_temp,'o',x,v_punto,'r')
     end
        
  
    [xt,lambda_ott]= polyxpoly(xi,yi,x,c_punto);
      if (length(xt)>1) % ci sono + intersezioni e va presa la + piccola   
        xt=min(xt);
        lambda_ott=min(lambda_ott);
      end
   
%    gamma_new = (lambda_ott-lambda_r)/(lambda_w-lambda_r);
%     if (gamma_new >=1 || gamma_new<=0) % controlla se è un valore di gamma tra 0 e 1 altrimenti old style
%         lambda_ott=gamma*lambda_w + (1 - gamma)*lambda_r;
%         prezzi_r(i,1)=lambda_r;
%         prezzi_w(i,1)=lambda_w;
%         temp=round(lambda_ott/int_lambda+1);
%         xt= round(c_punto_inv(temp)/int_x)*int_x;
%     end
%     xt= round(xt/int_x)*int_x;
    %temp=round(lambda_ott/int_lambda+1);

       if (abs(lambdaottimo-lambda_ott)<0.05*(lambdaottimo)&& flag==false);
            regime=i
            flag=true;
       end
    
    
    % aggiorna l'integrale
    y(i) = temp + (prezzi_r(i,1)-prezzi_w(i,1))*xt;
    temp = y(i);
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_w-prezzi_r);
title('Differenza');

figure(7);
plot(t,y);
title('Somme delle differenze di prezzo');

somma_persa=y(durata,1)   % valore a regime che perde/guadagna 

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



