急！请教matlab 排队论

M/M/m排队系统性能仿真代码matlab
function r=randexp(lambda) r = -lambda*log(rand);

说明同上

%mmm simulation in matlab
clc;clear;

ST_Idle=0;
ST_Busy=1;

EV_NULL=0;
EV_Arrive=1;
EV_Depart=2;
EV_LEN=3;

Q_LIMIT=1000;

% next_event_type=[];
% next_depart=[]
% num_custs_delayed=[];
% num_delays_required=[];
% num_events=[];
% num_in_q=[];
% server_status=[];
% area_num_in_q=[];
% area_server_status=[];
% mean_interarrival=[];
% mean_service=[];
% sim_time=[];
% time_last_event=[];
% total_of_delays=[];
%
time_arrival=[]; %到达时刻

time_next_event=zeros(1,EV_LEN);
%仿真参数
num_events=EV_LEN-1;
num_server=3; %M/M/m m=2
mean_interarrival=1;
mean_service=.5;
num_delays_required=2000; %

outfile=fopen("MM m.txt","w");
fprintf(outfile, "Multiple-server queueing system ");
fprintf(outfile, "Mean interarrival time%11.3f minutes ",mean_interarrival);
fprintf(outfile, "Mean service time%16.3f minutes ", mean_service);
fprintf(outfile, "Number of servers%20d ", num_server);
fprintf(outfile, "Number of customers%14d ", num_delays_required);

%part1 initialize
sim_time=0.0;
%/* Initialize the state variables. */

server_status =zeros(1,num_server); %idle
num_in_q = 0;
time_last_event = 0.0;

%/* Initialize the statistical counters. */

num_custs_delayed = 0;
total_of_delays = 0.0;
total_of_time = 0.0;

area_num_in_q = 0.0;
area_server_status = 0.0;

%/* Initialize event list. Since no customers are present, the departure
%(service completion) event is eliminated from consideration. */
time_next_event(EV_Arrive) = sim_time + randexp(mean_interarrival);
time_next_event(EV_Depart) = 1.0e+230;

time_depart=zeros(1,num_server);
next_depart=0;

%part2
while (num_custs_delayed < num_delays_required)
%Run the simulation while more delays are still needed.
%/* Determine the next event. */

min_time_next_event = 1.0e+290;
next_event_type = 0;
% Determine m depart event
min_time_depart=1e290;
next_depart=-1;
for i=1:num_server
if(server_status(i)==1 & time_depart(i)<min_time_depart)
min_time_depart=time_depart(i);
next_depart=i;
end
end
time_next_event(2)=min_time_depart;
%

%/* Determine the event type of the next event to occur. */

for i = 1: num_events
if (time_next_event(i) < min_time_next_event)
min_time_next_event = time_next_event(i);
next_event_type = i;
end
end

%/* Check to see whether the event list is empty. */

if (next_event_type == 0)

%/* The event list is empty, so stop the simulation. */

fprintf(outfile, " Event list empty at time %f", sim_time);
exit(1);
end
%/* The event list is not empty, so advance the simulation clock. */

sim_time = min_time_next_event;

%/* Update time-average statistical accumulators. */

double time_since_last_event;

%/* Compute time since last event, and update last-event-time marker. */

time_since_last_event = sim_time - time_last_event;
time_last_event = sim_time;

%/* Update area under number-in-queue function. */

area_num_in_q=area_num_in_q + num_in_q * time_since_last_event;

%/* Update area under server-busy indicator function. */

for i=1:num_server
area_server_status =area_server_status + server_status(i) * time_since_last_event;
end

%/* Invoke the appropriate event function. */
%arrival

if(next_event_type==EV_Arrive)
double delay;

%/* Schedule next arrival. */

time_next_event(1) = sim_time + randexp(mean_interarrival);

%/* Check to see whether server is busy. */

s_idle=-1;
for i=1:num_server
if (server_status(i) == ST_Idle)
s_idle=i;
break;
end
end

%/* all Server is busy, so increment number of customers in queue. */
if(s_idle== -1 )
num_in_q=1+num_in_q;

%/* Check to see whether an overflow condition exists. */

if (num_in_q > Q_LIMIT)
%/* The queue has overflowed, so stop the simulation. */

fprintf(outfile, " Overflow of the array time_arrival at");
fprintf(outfile, " time %f", sim_time);
exit(2);
end
%/* There is still room in the queue, so store the time of arrival of the arriving customer at the (new) end of time_arrival. */

time_arrival(length(time_arrival)+1)=sim_time;

else
%/* Server is idle, so arriving customer has a delay of zero. (The following two statements are for program clarity
%and do not affect the results of the simulation.) */

delay = 0.0;
total_of_delays =total_of_delays + delay;

%/* Increment the number of customers delayed, and make server busy. */

num_custs_delayed = 1 + num_custs_delayed;
server_status(s_idle) = ST_Busy;

%/* Schedule a departure (service completion). */

time_depart(s_idle) = sim_time + randexp(mean_service);
end % if (server_status == ST_Busy)
%depart
else
double delay;

%/* Check to see whether the queue is empty. */

if (num_in_q == 0)

% /* The queue is empty so make the server idle and eliminate the departure (service completion) event from consideration. */

server_status(next_depart) = ST_Idle;
time_depart(next_depart) = 1.0e+230;
%check_depart()
min_time_dapart=1e290;
next_depart=-1;
for i=1:num_server
if(server_status(i)==1 & time_depart(i)<min_time_depart)
min_time_depart=time_depart(i);
next_depart=i;
end
end
time_next_event(2)=min_time_depart;
else

%/* The queue is nonempty, so decrement the number of customers in queue. */

num_in_q=num_in_q-1;

%/* Compute the delay of the customer who is beginning service and update the total delay accumulator. */

delay = sim_time - time_arrival(1);
total_of_delays =total_of_delays + delay;

%/* Increment the number of customers delayed, and schedule departure. */

num_custs_delayed = 1 + num_custs_delayed;
serv_time=randexp(mean_service);
time_depart(next_depart) = sim_time + serv_time;
total_of_time = total_of_time + delay + serv_time;
%check_depart()
min_time_dapart=1e290;
next_depart=-1;
for i=1:num_server
if(server_status(i)==1 & time_depart(i)<min_time_depart)
min_time_depart=time_depart(i);
next_depart=i;
end
end
time_next_event(2)=min_time_depart;

%/* Move each customer in queue (if any) up one place. */

tempForPop=time_arrival(2:length(time_arrival));
time_arrival=tempForPop;
end %if (num_in_q == 0)

end %if(next_event_type==EV_Arrive)

end %while

%%%%%%%%%% part 3
%/* Invoke the report generator and end the simulation. */

fprintf(outfile, " Average delay in queue%11.3f minutes ",total_of_delays / num_custs_delayed);
fprintf(outfile, " Average delay in system%11.3f minutes ",total_of_time / num_custs_delayed);

fprintf(outfile, "Average number in queue%10.3f ",area_num_in_q / sim_time);
fprintf(outfile, "Server utilization%15.3f ",area_server_status / sim_time);
fprintf(outfile, "Time simulation ended%12.3f minutes", sim_time);
fclose(outfile);