function results=simulace_n(sim,lambda,mu,n)

%funkce simulace_n simuluje cinnost systemu M/M/n/Inf/FIFO

%syntaxe:

%results=simulace_n(sim,lambda,mu,n)

%vstupni parametry:

 %sim ... celkovy pocet simulaci

 %lambda ... intenzita vstupniho proudu zakazniku (tj. pocet zakazniku za

 %jednotku casu)

 %mu ... intenzita obsluhy (tj. pocet obslouzenych zakazniku za jednotku

 %casu)

 %n ... pocet linek v systemu

 %Pozn.: lambda musi byt mensi nez mu*n, jinak se system nestabilizuje

%vystupni parametry:

 %results ... obsahuje dalsi struktury dosazitelne pmoci teckove notace

  %.cas ... delky intervalu mezi prichody zakazniku

  %.obsluha ... delka obsluhy

  %.prichod ... cas prichodu zakaznika

  %.odchod ... cas odchodu zakaznika

  %.cekani ... 0 pokud zakaznik neceka, 1 pokud zakaznik ceka

  %.linka ... dulezite hlavne pro vypocet, v prvnim radku je cas, ve kterem

  %bude n-ta linka opet volna a v druhem radku je celkovy cast obsluhy v

  %n-te lince

  %.citac ... serazeny cas prichodu a odchodu

  %.citac ... pocita, kolik bylo v danem okamziku pritomno zakazniku

  %Pozn.: pomoci plot(results.poradi,results.citac,'*') lze vykreslit graf

  %celkoveho poctu zakazniku v systemu v prubehu casu 

  %.pst ... v prvnim radku obsahuje vsechny varianty poctu zakazniku a ve

  %druhem radku jejich pravdepodobnost

% Autor: Jakub Buček, 2014

%   **08/04/2015 Daniela Kuruczova: oprava vypoctu strednej doby cekani a obsluhy 

 

if lambda/mu/n>=1 error('System se nemuze stabilizovat. Intenzita provozu je vetsi nez 1.') 

end



% inicializace promennych

cas=zeros(1,n);

obsluha=zeros(1,n);

prichod=zeros(1,n);

odchod=zeros(1,n);

cekani=zeros(1,n);

linka=zeros(2,n);



for k=1:sim

    % vygenerovani casu mezi prichody a doby obsluhy

    cas(k)=exprnd(1/lambda);

    obsluha(k)=exprnd(1/mu);

    

    % pocitadlo casu, tj. okamziky prichodu zakazniku

    prichod(k)=sum(cas(1:k));

    

    % nejprve se podivame do 1. radku promenne linka, jestli je volna linka

    % pokud ano, tak ji obsadime zakaznikem

    % pokud ne, tak pricteme k lince, ktera bude nejdrive volna cas na

    % obsluhu posledniho zakaznika a zaznamename do promenne cekani, ze

    % zakaznik musel cekat

    [hodnota,pozice]=min(linka(1,:));

    if (hodnota<=prichod(k))

        linka(1,pozice)=prichod(k)+obsluha(k);

        linka(2,pozice)=linka(2,pozice)+obsluha(k);

        odchod(k)=prichod(k)+obsluha(k);

        cekani(k)=0;

    else

        linka(1,pozice)=linka(1,pozice)+obsluha(k);

        linka(2,pozice)=linka(2,pozice)+obsluha(k);

        odchod(k)=hodnota+obsluha(k);

        cekani(k)=1;

    end   

end



% hledani pravdepodobnosti

[poradi,index]=sort([0,prichod,odchod]);

m=length(poradi);

citac=zeros(1,m);



for k=1:m

    for h=1:sim

        if (prichod(h)<=poradi(k)) && (odchod(h)>poradi(k))

        citac(k)=citac(k)+1;

        end

    end

end



citac=citac(1:(end-1));



varianty=unique(citac);

intervaly=diff(poradi);

m=length(varianty);

pst=zeros(1,m);



for k=1:m

    pst(k)=sum(intervaly(citac==varianty(k)));

end



pst=[varianty;pst/odchod(sim)];



results.cas=cas;

results.obsluha=obsluha;

results.prichod=prichod;

results.odchod=odchod;

results.cekani=cekani;

results.linka=linka;

results.poradi=poradi;

results.citac=[citac,0];

results.pst=pst;



disp('**************************************************************************')

fprintf('Prumerna doba mezi prichody zakazniku: %6.3f \n',mean(cas));

fprintf('Prumerna doba obsluhy: %6.3f \n',mean(obsluha));

fprintf('Prumerna doba cekani ve fronte: %6.3f \n',mean(odchod-prichod)-mean(obsluha));

fprintf('Prumerna doba v systemu: %6.3f \n',mean(odchod-prichod));

fprintf('Prumerny pocet zakazniku v systemu: %6.3f \n',pst(1,2:m)*pst(2,2:m)');

fprintf('Prumerny pocet zakazniku ve fronte: %6.3f \n',(pst(1,(n+2):m) - n)*pst(2,(n+2):m)');

fprintf('Prumerny pocet obsluhovanych zakazniku: %6.3f \n',(pst(1,2:(n+1))*pst(2,2:(n+1))') + n * sum(pst(2,(n+2):m)));

fprintf('Intenzita provozu: %6.3f \n',sum(linka(2,:))/(n*odchod(sim)));

fprintf('Celkova doba cinnosti linek: %6.3f \n',sum(linka(2,:)));

fprintf('Celkova doba prostoje linek: %6.3f \n',n*odchod(sim)-sum(linka(2,:)));

fprintf('Pravdepodobnost, ze prichazejici zakaznik bude cekat ve fronte: %6.3f \n',mean(cekani));

disp('**************************************************************************')

for k=1:m

fprintf('Pravdepodobnost, ze v systemu bude %2u zakazniku je %6.3f \n', [pst(1,k) pst(2,k)]);

end

disp('**************************************************************************')