nIterazioni = 10000;  		% Il  numero di iterazioni, cioè il numero di volte che viene scelta la nuova parte da lavorare
delta = 1; % Il tempo di setup delta
tau = [1/35;1/7;1/20];		% I tempi di lavorazione dei vari tipi di parti: tau = [tau1, tau2, ...]
% tau = [1; 1; 1]/40;
d = [18;3;1];				% I tassi di domanda: d = [d1, d2, ...]
% d = [10; 10; 10];
P = numel(d);			% Il numero P di tipi di parti
rhoV = d.*tau;			% Il vettore dei rhoj=dj tauj, cioè rhoV = [rho1, rho2, rho3, ...]
rho = sum(rhoV);		% Il parametro rho somma dei rhoj
xOld = 0*ones(P,1);		% xOld è il vettore x = [x1, x2, ...] del contenuto dei buffer all'inizio del nuovo ciclo. Qui
						% rappresenta lo stato iniziale x0
xNew = xOld;  			% xNew invece è il valore di x alla fine del ciclo
tOld = 0;				% tempo all'inizio del ciclo (Tn)
wipOld = sum(xOld);	   	% wipOld è il WIP corrispondente a xOld
integraleWIP = 0;		% integraleWIP è l'integrale del WIP da 0 al tempo attuale
mediaWIP = zeros(nIterazioni,1);   % mediaWIP è un vettore la cui i-esima componente dà il valore medio del WIP da 0 fino alla fine del ciclo i
tempi = zeros(2*nIterazioni,1);	   % vettore contenente tutti gli istanti Tn in cui sono state scelte le parti da lavorare e i Tn+delta
tempiTn = zeros(nIterazioni,1);	   % vettore contenente tutti gli istanti Tn in cui sono state scelte le parti da lavorare
xBuffer = zeros(P,2*nIterazioni);  % vettore che memorizza x ad ogni tempo Tn, n>=1
scelte = zeros(nIterazioni,1);     % vettore che memorizza tutte le scelte fatte dalla politica adottata (in questo caso CLB)

parteDaLavorare = 0; % Serve solo per simulare la RR

tic

for i=1:nIterazioni,
    
    % Scelta politica: scommentare quella da simulare
     [massimo,iMax] = max(xOld); %CLB
%      [massimo,iMax] = max(xOld+d.*delta); %CLAB
%      [massimo,iMax] = max(tau.*xOld); %CLW
%      [massimo,iMax] = max(tau.*xOld+rhoV.*delta); %CLAW
%      [massimo,iMax] = max((tau.*xOld+delta)./(1-rhoV));	% CLAT
%      [massimo,iMax] = max((xOld+delta*d)./(sqrt((1-rhoV).*d)));  % CLAR (con ci = 1)
%      parteDaLavorare = mod(parteDaLavorare,P)+1; iMax = parteDaLavorare;  %RR
     
   tNew = tOld + (delta+tau(iMax)*xOld(iMax))/(1-rhoV(iMax)); % è la formula (1.13) del testo (Tn+1 in funzione di Tn)
   xDelta = xOld + d*delta;  % il contenuto dei buffer x al tempo Tn + delta
   xNew = xDelta + d*(tNew-tOld-delta);		% il contenuto dei buffer x al tempo Tn+1 .... continua con la riga successiva
   xNew(iMax) = xDelta(iMax) + (d(iMax)-1/tau(iMax))*(tNew-tOld-delta); %deve essere = 0, controllato nell'istruzione seguente
   if xNew(iMax)>0.001, 
       disp('Errore');
       return
   end
   tempi(2*i-1) = tOld + delta;	     % si memorizza in tempi il tempo Tn + delta
   tempi(2*i) = tNew;				 % si memorizza in tempi il tempo Tn+1
   tempiTn(i) = tNew;                % si memorizza in tempiTn il tempo Tn+1 
   xBuffer(:,2*i-1) = xDelta; % si memorizza in xBuffer il contenuto dei buffer x al tempo Tn + delta   
   xBuffer(:,2*i)=xNew;		  % si memorizza in xBuffer il contenuto dei buffer x al tempo Tn+1
   wipDelta = sum(xDelta);		      % si calcola il WIP corrispondente a Tn + delta
   wipNew = sum(xNew);		      % si calcola il WIP corrispondente a Tn+1
   integraleWIP = integraleWIP + (wipOld+wipDelta)*delta/2 + (wipDelta+wipNew)*(tNew-tOld-delta)/2;	% si aggiorna l'integrale del WIP
   wipOld = wipNew;		  	   % si aggiorna wipold
   tOld = tNew;					% si aggiorna told
   xOld = xNew;				   % si aggiorna xOld
   mediaWIP(i) = integraleWIP/tNew; % si calcola il WIP medio fino a questo punto e lo si memorizza in mediaWIP
   scelte(i) = iMax;				% si memorizza la scelta fatta nell'ultimo ciclo
end

toc

WIPmedio = mediaWIP(nIterazioni)		% si restituisce la media finale del WIP
LB = delta*sum(sqrt(d.*(1-rhoV)))^2/(2*(1-rho)) % si fornisce il lower bound LB del costo J dato riportato nel teorema 1.4 del testo (se ci=1)
WIPmedioRR = 0.5*(P*delta/(1-rho))*sum(d.*(1-rhoV)) % WIP medio di una RR (formula (1.26) testo con ci=1)
fTeorica = (1-rho)/delta  % frequenza a regime di una politica stabile non-idling (teorema 1.2 testo)
fSimRegime = nIterazioni/tempi(2*nIterazioni)  % frequenza setup a regime relativa alla simulazione

figure
colori = ['r','g','b','m','c','k','y'];
nColori = numel(colori);
legenda = cell(1,P);
legenda(1) = cellstr('x1');
plot(tempi,xBuffer(1,:),'r');
hold on
for indParte = 2:P,
    colore = colori(mod(indParte-1,nColori)+1);
    plot(tempi,xBuffer(indParte,:),colore)
    legenda(indParte) = cellstr(['x' num2str(indParte)]);
end
legend(legenda)
xlabel('Tempi')

figure
plot(tempi,tau'*xBuffer)
xlabel('Tempi [s]')
ylabel('Lavoro w(t)')

figure
plot(tempi,sum(xBuffer))
xlabel('Tempi [s]')
ylabel('Work In Process WIP(t)')