This program simulates a server model with 4 servers that handle requests that arrive at random times. Rewrite this program in Python (or Java), except that instead of, S = sum( -1/lambda(u) * log(rand(alpha(u),1) ) );, use library functions to generate a normal distribution.

Let us simulate one day, from 8 am till 10 pm, of a queuing system that has

four servers;

Gamma distributed service times with parameters given in the table,

Server ? ?

I 6 0.3 min?1

II 10 0.2 min?1

III 7 0.7 min?1

IV 5 1.0 min?1

a Poisson process of arrivals with the rate of 1 arrival every 4 min, independent of

service times;

random assignment of servers, when more than 1 server is available.

In addition, suppose that after 15 minutes of waiting, jobs withdraw from a queue if their

service has not started.

For an ordinary day of work of this system, we are interested to estimate the expected

values of:

the total time each server is busy with jobs;

the total number of jobs served by each server;

the average waiting time;

the longest waiting time;

the number of withdrawn jobs;

the number of times a server was available immediately (this is also the number of

jobs with no waiting time);

the number of jobs remaining in the system at 10 pm.

k = 4; % number of servers

mu = 4; % mean interarrival time

alpha = [6 10 7 5]; % parameters of service times

lambda = [0.3 0.2 0.7 1.0];

arrival = [ ]; % arrival time

start = [ ]; % service starts

finish = [ ]; % service finishes; departure time

server = [ ]; % assigned server

j = 0; % job number

T = 0; % arrival time of a new job

A = zeros(1,k); % array of times when each server

% becomes available for the next job

while T < 840; % until the end of the day

j=j+1; % next job

T = T-mu*log(rand); % arrival time of job j

arrival = [arrival T];

Nfree = sum( A < T ); % number of free servers at time T

u = 1; % u = server that will take job j

if Nfree == 0; % If all servers are busy at time T ...

for v=2:k;

if A(v)< A(u); % Find the server that gets available

u = v; % first and assign job j to it

end;

end;

if A(u)-T > 15; % If job j waits more than 15 min,

start=[start T+15]; % then it withdraws at time T+15.

finish=[finish T+15];

u = 0; % No server is assigned.

else % If job j doesnt withdraw ...

start = [start A(u)];

end;

else % If there are servers available ...

u = ceil(rand*k); % Generate a random server, from 1 to k

while A(u) > T; % Random search for an available server

u = ceil(rand*k);

end;

start = [start T]; % Service starts immediately at T

end;

server = [server u];

if u>0; % if job j doesnt withdraw ...

S = sum( -1/lambda(u) * log(rand(alpha(u),1) ) );

finish = [finish start(j)+S];

A(u) = start(j) + S; % This is the time when server u

end; % will now become available

end; % End of the while-loop

disp([(1:j) arrival start finish server])

for u=1:k; % Twork(u) is the total working

Twork(u) = sum((server==u).*(finish-start)); % time for server u

Njobs(u) = sum(server==u); % number of jobs served by u

end;

Wmean = mean(start-arrival); % average waiting time

Wmax = max(start-arrival); % the longest waiting time

Nwithdr = sum(server==0); % number of withdrawn jobs

Nav = sum(start-arrival < 0.00001); % number of jobs that did not wait

Nat10 = sum(finish > 840); % number of jobs at 10 pm