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

// State
volume_max=100;
volume_min=0;

// Control
control_max=floor(0.4/7*volume_max);
control_max=control_max+1;
control_min=0;

// Time
t0=1;
horizon=52;
T=t0:(horizon-1);
bT=t0:(horizon);

// Costs
cost=-66*2.7;
fcost=0;
cost=cost*(1+0.5*(rand(T)-1/2));

// Uncertainties
uncertainty_max=floor(0.5/7*volume_max) ;
uncertainty=[0:uncertainty_max];

// Probabilities
unnormalized_proba=cumsum(ones(1,uncertainty_max+1))-1;
proba1=unnormalized_proba/sum(unnormalized_proba);
// more rain in winter
proba182=proba1($:-1:1);
// less rain in summer

// simulation of K independent sequences of water inflows
//between 1 and uncertainty_max+1 

// ----------------------------------------------------
//  Dynamical System Functions
// ----------------------------------------------------


// Dynamics
function sdot=dynamics(s,u,d,a)
sdot=s-u-d+a;
endfunction

// Instantaneous cost function
function c=instant_cost(t,s,u,d,a)
c=cost(t)*u;
endfunction

// Final cost function
function c=final_cost(s)
c=fcost*s;
endfunction  // here the final value of water is 0


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

// inflow probability modelization
function prob = gen_proba()
  // the vector 'proba' of probabilities differs according to time t
  // more rain in winter, less rain in summer
prob=zeros(uncertainty_max+1,horizon-t0+1);
for t=bT
  prob(:,t)=((1-sin(%pi*t/(horizon-t0+1)))*proba1 + sin(%pi*t/(horizon-t0+1))*proba182)';
end
endfunction

// Scenario generation
function Scenarios= gen_scenario(nb,proba)
  // The function generate a number of scenarios of external inflows
Scenarios=zeros(nb,horizon-t0+1);
for k=1:nb
    for t=bT
    Pcum=cumsum((proba(:,t))');
    a=rand(1);
    j=1
    while a>Pcum(1,j)
      j=j+1;
    end
    Scenarios(k,t)=j-1; 
  end
end
endfunction

// Generating the probability distribution of inflows
proba = gen_proba();

// Trajectories simulations 

function [S,U,D,C]=trajectories(s0,policy,scenarios)
// Compute the trajectories of the system over a set of scenarios
// according to a given trajectories

S=[];// Stock trajectories
U=[];// Control trajectories
D=[];// Control (spillage) trajectories
C=[];// total costs


nb_simulations=size(scenarios,'r');

for ksimu=1:nb_simulations

  s=s0;
  u=[];
  d=[];
  c=0;
  
  for t=T
    [turb,spil]=policy(t,s($),scenarios(ksimu,t));
    u=[u,turb];
    d=[d,spil];
    s=[s,dynamics(s($),u($),d($),scenarios(ksimu,t))] ;
    c=c+ instant_cost(t,s($),u($),d($),scenarios(ksimu,t)) ; 
  end
c=c+final_cost(s($)) ; 

  S=[S;s];
  U=[U;u];
  D=[D;d];
  C=[C;c];
end
endfunction