% Simulazione_Attianese_et_al.m
%
% Simulazione di Monte Carlo per il caso di studio  "5 MIXTURES"
%
%
% $Revision: 1                  $Date: Sept. 24, 2004
% $Revision: 2.1 by M.Miccio    $Date: May 31, 2007

N=input ('Qual è la durata della simulazione in giorni? ');
disp (' ')
disp ('Qual è il N. di miscelatori per la simulazione ? ')
mixers=input('(per default è N. mixers = 13)      N = ');
if isempty(mixers) mixers=13; end;

rand('state',0)    % ---Resets the random No. generator to its initial state----- 

for k=1:N
   %%%%%%%%%%%%%%%%%%%%%%%%%%MISCELA A%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
   uni=rand*1000;
   if (uni>=0 & uni<=20), y(1,k)=1; end
   if (uni>20 & uni<=223), y(1,k)=2; end
   if (uni>223 & uni<=569), y(1,k)=3; end
   if (uni>569 & uni<=793), y(1,k)=4; end
   if (uni>793 & uni<=924), y(1,k)=5; end
   if (uni>924 & uni<=968), y(1,k)=6; end
   if (uni>968 & uni<=995), y(1,k)=7; end
   if (uni>995 & uni<=999), y(1,k)=8; end
   %%%%%%%%%%%%%%%%%%%%%%%%%%MISCELA B%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
   uni=rand*1000;
   if (uni>=0 & uni<=142), y(2,k)=1; end
   if (uni>142 & uni<=632), y(2,k)=2; end
   if (uni>632 & uni<=898), y(2,k)=3; end
   if (uni>898 & uni<=985), y(2,k)=4; end
   if (uni>985 & uni<=999), y(2,k)=5; end
   %%%%%%%%%%%%%%%%%%%%%%%%%%MISCELA C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
   uni=rand*1000;
   if (uni>=0 & uni<=3), y(3,k)=1; end
   if (uni>3 & uni<=90), y(3,k)=2; end
   if (uni>90 & uni<=527), y(3,k)=3; end
   if (uni>527 & uni<=901), y(3,k)=4; end
   if (uni>901 & uni<=988), y(3,k)=5; end
   if (uni>988 & uni<=999), y(3,k)=6; end
   %%%%%%%%%%%%%%%%%%%%%%%%%%MISCELA D%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
   uni=rand*1000;
   if (uni>=0 & uni<=20), y(4,k)=1; end
   if (uni>20 & uni<=607), y(4,k)=2; end
   if (uni>607 & uni<=978), y(4,k)=3; end
   if (uni>978 & uni<=999), y(4,k)=4; end
   %%%%%%%%%%%%%%%%%%%%%%%%%%MISCELA E%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
   uni=rand*1000;
   if (uni>=0 & uni<=366), y(5,k)=1; end
   if (uni>366 & uni<=905), y(5,k)=2; end
   if (uni>905 & uni<=992), y(5,k)=3; end
   if (uni>992 & uni<=999), y(5,k)=4; end

   KK=k;
   disp ('giorno della simulazione: ')
   disp (KK)
   disp (' ')
   
        for i=1:5
        disp (y(i,KK))
        end
   continua=input('Premi INVIO per continuare','s');
end
%_-------------Somma degli ordini giornalieri--------------
% For matrices, SUM(X) is a row vector with the sum over each column. 
sum_day=sum(y);

sum_orders=sum(sum_day);
NB=(sum_day>mixers).*(sum_day-mixers);    % N. ordini inevasi (No. of orders - No. of mixers) each day of simulation

%--------------Calcolo del N. di giorni in cui si verifica la stessa richiesta di ordini  
%              (trascurando il limite al numero dei miscelatori)
for k=0:25
   h(k+1)=length(find(y==k));
end

%-------------Calcolo del N. di giorni in cui si verifica la stesso N. di ordini inevasi 
%             (considerando il numero di miscelatori)
Nmax=30;
for k=0:Nmax
   p(k+1)=length(find(NB==k));
end
%--------------Calcolo della Somma degli ordini inevasi -------------
sum_back_orders=0;
for k=0:Nmax
   sum_back_orders=k*p(k+1)+sum_back_orders;
end

sum_orders   %---------------------------------Somma totale degli ordini ---------------------
sum_back_orders        
percentuale_back_orders=(sum_back_orders/sum_orders)*100   %---Percentuale degli ordini inevasi nello stesso giorno
percentuale_orders_no_delay=(1-sum_back_orders/sum_orders)*100   %---Percentuale degli ordini evasi senza ritardo-----------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Grafico "numero ordini non evasi vs. giorni" %%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
stem([0:Nmax],p,'o')
grid
xlabel('No. of back orders')
ylabel('Days')

disp('... The end ...')
disp('  ')


