%=================== CDA6530, Fall 2010 =============================
%========= discrete event simulation of  =========== 
%========programming project 3 ======================================
%========Use realistic simulation considering the actions of all ===
%======== 3 machines                                             ===
%====== corresponding to slides Project3-solution.ppt ==============

% first, comput the theoretical value
lambda = 0.2;  mu = 0.5; c = 3;
Pi_theory = zeros(c+1, 1);
k=0:c;
Pi_theory(1) = 1/sum( power(mu/lambda, k)./factorial(k) );  %Pi_0
k=1:c;
Pi_theory(k+1) = power(mu/lambda, k)./factorial(k)* Pi_theory(1);
eta_theory = 1 - Pi_theory(c+1)
avgWorkingNum_theory = mu/lambda * eta_theory;

%-------------- discrete event simulation --------------------
Queue = zeros(c,1);  Status = zeros(c,1); 
EventList = zeros(c,1);
finishRepairNum = 100;  
Tran = zeros(10*finishRepairNum,1); 
SystemState = zeros(10*finishRepairNum,1); 
% initialization
queueNum = 0;  currentTime = 0; repairNum = 0; transitionNum = 0;
workingNum = c;
for k=1:c, 
    Status(k) = 'w';   %all machines are working at beginning
    EventList(k) = -log(rand)/lambda + currentTime;
end

while repairNum < finishRepairNum,
    [currentTime, k] = min(EventList);
    switch Status(k),
        case 'w'  % k-th machine breaks
            % update queue, put k-th machine at the end of queue
            queueNum = queueNum + 1;
            Queue(queueNum) = k;
            workingNum = workingNum - 1;
            if queueNum == 1,  % immediately repaired by repairman
                Status(k) = 'r';
                EventList(k) = currentTime - log(rand)/mu;
            else   % the machine waits in the queue
                Status(k) = 'q';
                EventList(k) = inf;
            end
        case 'r'   %  finishing repair, Queue(1) must be the k-th machine
            % update queue - move all waiting jobs one step ahead
            repairNum = repairNum + 1;
            Status(k) = 'w';
            EventList(k) =  -log(rand)/lambda + currentTime;
            workingNum = workingNum + 1;
            if queueNum >= 2,  % still have waiting job in queue
                for i=2:queueNum,
                    Queue(i-1) = Queue(i);
                end
                queueNum = queueNum - 1;
                % the job of Queue(1) should change from 'q' to 'r'
                j = Queue(1);
                Status(j) = 'r'; 
                EventList(j) = currentTime - log(rand)/mu;
            else   %  Queue becomes empty after k-th job finishes repair
                queueNum = queueNum - 1;
            end
        case 'q'  % should never have this case
            error('wrong status appeared!');
    end
    transitionNum = transitionNum + 1;
    Tran(transitionNum) = currentTime;
    SystemState(transitionNum) = workingNum + 1 ;
end

%  comput simulation results
Time = zeros(c+1,1); Time(c+1) = Tran(1); 
for k=1:c+1,
    for i=1:transitionNum-1,
        if SystemState(i) == k,
            Time(k) = Time(k) + Tran(i+1) - Tran(i);
        end
    end
end
Pi_simul = zeros(c+1,1);
Pi_simul = Time./Tran(transitionNum);
eta_simul = 1 - Pi_simul(c+1);
avgWorkingNum_simul = sum(Pi_simul.* 0:c);



figure; hold on; 
plot(Pi_theory, '-ob'); plot(Pi_simul, '--^k');
xlabel('Steady state probability \pi'); 
legend('Theoretical', 'Simulation');





