///////////////////////////////////////////////////////////////////////
//    STOCHASTIC DUAL DYNAMIC PROGRAMMING EQUATION
///////////////////////////////////////////////////////////////////////

clear ;
exec damdata.sce ;
exec CutPerAlea.sce ;

// ----------------------------------------------------
// DATA 
// ----------------------------------------------------

// MODEL PARAMETERS
// initial state
x0=50;
// dynamics
Q=[1];
R=[-1,-1];
// constraints (see 3.5)
b=[-control_min;control_max;0;0;volume_max];
A=[-1,0;1,0;0,-1;1,1;-1,-1];
E=[0;0;0;-1;1];


// ALGORITHM PARAMETERS

// Initialization
L=3; // Forced forward loops, the cut is computed where we want it to
// Number of cuts computed before stopping the algorithm
optimal_cuts=50;
// Number of scenarios to evaluate final policy
nb_eval = 1000


// ----------------------------------------------------
// MACROS 
// ----------------------------------------------------


// APPROXIMATE OF THE BELLMAN FUNCTION STOCHASTIC CASE

function V=SDDP(A,E,b,Q,R,L,oc)
  // V is the approximation of the Bellman function obtained. 
  // It is a collection of hyperplane, the slope being the first column,
  // and the ordinate at 0 the second column.

  // A is the inequality matrix
  V=cell(1,horizon);
  for t=t0:horizon
    V{1,t}=zeros(1,2);
    V{1,t}(1,2)=-1/%eps;
  end
  V{1,horizon}=[fcost,0];

  // Stock in the dam
  // The initial stock is x0
  x=zeros(oc,horizon);
  x(:,1)=x0;


// We first compute a cut for Vt at L points
  delta=(volume_max-volume_min)/(L+1);
  for l=1:L
     for t=horizon-1:-1:t0  
       for i=1:uncertainty_max+1
 	b(4:5,1)=[uncertainty(1,i);volume_max-uncertainty(1,i)];
 	[LAMBDA(i,t),BETA(i,t),u_opt]=CutPerAlea([cost(1,t),0],A,b,E,...
			uncertainty(1,i),Q,R,V{1,t+1},volume_min+l*delta);
       end
       //computation of the new plane
       V{1,t}=[V{1,t};[proba(:,t)'*LAMBDA(:,t),proba(:,t)'*BETA(:,t)]];
     end
  end

  // THE MAIN LOOP is computed oc times
  // oc is chosen in advance

  for k=1:oc
    // FORWARD PHASE
    // The optimal trajectory is computed thanks to
    // the current approximate of the Bellman function
    
    // Computing noise scenario for forward phase
    inflow_scenar = gen_scenario(1,proba)

    for t=1:horizon-1
      //the constraint vector is actualized
      b(4:5,1)=[inflow_scenar(1,t);volume_max-inflow_scenar(1,t)];
      //Use of an external linear solver through CutPerAlea, to obtain the optimal control
      [lambda,beta,u_opt]=CutPerAlea([cost(1,t),0],A,b,E, inflow_scenar(1,t),Q,R,V{1,t+1},x(k,t));
      //Computation of the new stock
      x(k,t+1)=dynamics(x(k,t),u_opt(1,1),u_opt(2,1), inflow_scenar(1,t));
    end

    // BACKWARD PHASE
    // a new cut is built near the optimal solution at each step

    for t=horizon-1:-1:1  
      for i=1:length(uncertainty) // Computing the parameters per uncertainty
	//The constraint vector is actualized
	b(4:5,1)=[uncertainty(1,i);volume_max-uncertainty(1,i)];
	//Use of an external linear solver through CutPerAlea
	[LAMBDA(i,t),BETA(i,t),u_opt]=CutPerAlea([cost(1,t),0],A,b,E,...
				uncertainty(1,i),Q,R,V{1,t+1},x(k,t));
      end
      
      //computation of the new hyperplane parameters by taking into account
      //all the uncertainties
      V{1,t}=[V{1,t};[proba(:,t)'*LAMBDA(:,t), proba(:,t)'*BETA(:,t)]];
    end
  end

endfunction


//V=SDDP(A,E,b,Q,R,L,optimal_cuts);

// POLITIQUE SDDP
// Warning : use the result of the computation of SDDP. 
function [turbined,spilled]=sddp_policy(t,s,a)
  if t==horizon-1 then
    b=[-control_min;control_max;0;a+s;volume_max-a-s];
    C=[cost(1,t)-fcost,-fcost];
    [u_opt,c_opt,flag,d]=linprog(C,A,b,[],[],lb=[-%inf,-%inf]);
  else
    b=[-control_min;control_max;0;a;volume_max-a];
    [lambda,beta,u_opt]=CutPerAlea([cost(1,t),0],A,b,E,a,Q,R,V{1,t+1},s);
  end
  turbined = u_opt(1,1);
  spilled = u_opt(2,1);
endfunction

function res = V_function(t,x) 
// Tranform the collection of hyperplane in a function. 
// Depends on the last V computed
   res = max(V{t}(:,1)*x + V{t}(:,2))
endfunction

// ----------------------------------------------------
//  SIMULATIONS
// ----------------------------------------------------
 
// Trajectories simulations and visualization

//Scenarios= gen_scenario(nb_eval,proba);

//[S,U,D,C]=trajectories(x0,sddp_policy,Scenarios);
//xset("window",31); xbasc();
//plot2d(bT,S');
//xtitle('Stock volumes in a dam following an optimal policy with a zero final value of water','(time)','(volume)');

//xset("window",32); xbasc();
//plot2d2(T,U');
//xtitle('turbined volumes in a dam following an optimal policy with a zero final value of water','(time)','(volume)');



//xset("window",33);xbasc();
//plot2d2(T,D');
//xtitle('spilled volumes in a dam following an optimal policy with a zero final value of water','(time)','(volume)');

print('SDDP algorithm with '+string(optimal_cuts)+' iterations evaluated on '+string(nb_eval)+'scenarios.') ;
print('The optimal payoff given by simulation is '+string(mean(C)));