PHP Code:
/*
C:\Fall2003Henderson\WorkFolder\IpaVar3.exe 1.5 2 150 5 4 11 17 20 20 17 0.5 1 1.3 1.5 1.3 1.1
C:\Fall2003Henderson\WorkFolder\IpaVar5.exe 1.5 999 1440 24 4 11 20 15 15 26 20 22 22 23 27 27 27 27 27 33 34 51 38 43 43 32 38 30 29 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.3 1.3 1.3 1.3 1.45 1.6 1.75 1.9 2.05 2.05 2.05 2.05 2.05 1.95 1.85 1.75 1.65
*/
/*
This version implements version 2 of the IPA algorithm
*/
#include<stdio.h>
#include "RNG.h"
#include "RngStream.h"
#include<math.h>
#include<fstream.h>
#include<vector>
#include<list>
#include "Customer.hpp"
using namespace std;
// 3 types of events- ARRIVAL = 0, COMPLETION = 1, ........
enum Event {ARRIVAL, COMPLETION, PERIODEND};
/**********************Functions***********************************************/
void ReadCommandLine(char *pArgs[]);
void CalculateStats();
void WriteStatsFile();
void Simulate();
double RateFunct(double time);
void Initialize();
Event Transition();
void Arrival();
void Completion();
void PeriodEnd();
void Termination();
void UpdateIPA();
double Integral(double a, double b);
double NHPP(RngStream &generator, double time);
double Exponential(RngStream &generator, double lambda);
void ShowParameters();
void MakeTrace(Event event);
void PeriodEndTrace();
void ErrorTrace(fstream &DataFile);
void CheckIPACounters();
typedef unsigned short int USHORT;
typedef unsigned long int ULONG;
bool traceOn = false; //boolean to turn on and off trace. TRUE to turn trace on.
fstream TraceFile("Trace.txt",ios::in | ios::out |ios::trunc);
double maxArrivalRate; // maximum of arrival rate function over entire simulation
vector<double> arrivalRate; //holds piecewise arrival rate, arrivalRate[i] = arrival rate at END of period i
vector<int> numServers; //number of servers in each period
USHORT maxServers; //max number of servers in any period
double serviceRate; // mu, the rate of service (services per minute)
double simTime = 0; // simulation time in minutes
double maxWait; //maximum waiting time in minutes to count as success
double endTime = 0; // time that simulation terminates (in minutes from start of sim)
ULONG numCallsTotal = 0; // number of incoming calls in total
ULONG numCallsServiced = 0; // number of calls that have started service
ULONG numCallsInSystem = 0; // number of calls in system now
double timeToPeriodEnd; //time to end of period
double periodLength; // length of each period in minutes
USHORT currentPeriod = 1; //current period
USHORT totalPeriods; // total number of periods
USHORT totalReplications; // number of replications
USHORT currentReplication =1; // current replication
double timeToNextArrival = 0; // time to next arrival
USHORT serverLeastWork = 0; //index of server with least amount of work among all busy servers,
//=0 if all servers free
vector<double> workRemain; //holds remaining work time of each server
vector < vector<double> > delta; //holds accumulated infite delay, delta[server][period]
vector<double> T; //time that server i last completed a service, even if server i busy
vector< vector<double> > delKs; //array of nperiods partial derivatives for call Ks. delKs[Ks][period]
vector<double> xiKs; //holds time when a server became available to serve call Ks
vector<double> Ak; //holds time of kth arrival
vector< vector<double> > meanIPAGrad; //the mean of the IPA gradients obtained over all replications
vector< vector<double> > varIPAGrad;
vector<ULONG> successCounter; //number of successful services during period during each replication
vector<ULONG> numCalls; //number of arrivals during each period during each replication
vector<double> meanSuccess; //average number of succeses during each period
vector<double> varSuccess; //variance of the number of successes during each period
vector<double> meanNumCalls; //average number of arrivals during each period
list<Customer*> queue; //the arrival queue, list of pointers to Customers in queue
// Random number Streams
RngStream arrivalStream; // stream used for generating interarrival times in NHPP
RngStream workStream; // stream used for generating exp(1) r.v for workload
int main(int argc, char *pArgs[])
{
//Read in values from command line arguments
//This version takes in command line arguments
// 1st Arg: maximum wait time
// 2nd Arg: Number of Replications
// 3rd Arg: end_time: length in minutes of simulated period
// 4th Arg: Number of Periods
// 5th Arg: Service Rate (services per HOUR per server) HOUR!!! not minutes!!!
// Vector containing number of servers in real life (n entries)
// Vector containing piecewise arrival rate function (n+1 entries)
//ReadCommandLine(*argv[0]);
maxWait = atof(pArgs[1]);
totalReplications = atoi(pArgs[2]);
endTime = atof(pArgs[3]);
totalPeriods = atoi(pArgs[4]);
periodLength=endTime/totalPeriods;
timeToPeriodEnd = periodLength;
serviceRate = atof(pArgs[5])/60.0; //convert it into (services/minutes)/server
numServers.resize(totalPeriods+1);
cout<<"NumServersSize: "<<numServers.size()<<endl;
for(int i=0; i<totalPeriods;i++){
cout<<i+1<<"\t"<<pArgs[i+6]<<endl;
numServers[i]=atoi(pArgs[i+6]);
cout<<i+1<<"\t"<<numServers[i]<<endl;
if(numServers[i]>maxServers) maxServers = numServers[i];
}
//resize arrivalRate[] and initialize values
//assign value to maxArrivalRate
arrivalRate.resize(totalPeriods+1);
for(int i=0; i<totalPeriods+1;i++){
arrivalRate[i]=atof(pArgs[i+totalPeriods+6]);
if(arrivalRate[i]>maxArrivalRate){maxArrivalRate = arrivalRate[i]; }
}
// Resizing statistic counters
successCounter.resize(totalPeriods);
meanSuccess.resize(totalPeriods);
varSuccess.resize(totalPeriods);
meanNumCalls.resize(totalPeriods);
numCalls.resize(totalPeriods);
workRemain.resize(maxServers);
T.resize(maxServers);
delta.resize(maxServers);
for(int i=0; i<maxServers;i++){
delta[i].resize(totalPeriods);
}
vector<double> zeroVect(totalPeriods,0);
meanIPAGrad.resize(totalPeriods);
varIPAGrad.resize(totalPeriods);
for(int i=0;i<totalPeriods;i++){
meanIPAGrad[i]=zeroVect;
varIPAGrad[i]=zeroVect;
}
ShowParameters();
char duh ;
cin>>duh;
//***************END of Parameter initialization***************//
//this is the main simulation loop that loops through all the simulations.
for(; currentReplication<=totalReplications;currentReplication++){
Simulate();
CalculateStats();
}
TraceFile.close();
WriteStatsFile();
Initialize();
meanIPAGrad.clear();
varIPAGrad.clear();
meanNumCalls.clear();
varSuccess.clear();
meanSuccess.clear();
cout<<"DONE"<<endl;
cin>>duh;
return 1;
} //End Main
/******************************************************************************/
void CalculateStats(){
//This function updates the current mean and variance of successful services
//in each period.
//One pass variance/mean calculation algorithms obtained from
//http://mathworld.wolfram.com/SampleVarianceComputation.html
double OldMean;
double OldArrivalMean;
//accumulate service level stats
for(int l=0;l<totalPeriods;l++){
if(currentReplication==1){
meanSuccess[l]=successCounter[l]; //used for 1st moment
meanNumCalls[l]=numCalls[l];
}
else {
OldMean = meanSuccess[l];
OldArrivalMean = meanNumCalls[l];
meanSuccess[l]=((currentReplication-1)*OldMean+successCounter[l])/(currentReplication); //Sample mean after (repnumber+1)th replication
meanNumCalls[l]=((currentReplication-1)*OldArrivalMean+numCalls[l])/(currentReplication); //Sample mean after (repnumber+1)th replication
varSuccess[l] = (1.00-(1.00/(currentReplication-1)))*varSuccess[l] + (currentReplication)*pow((meanSuccess[l]-OldMean),2);
}
}
}//End CalculateStats()
/******************************************************************************/
void Simulate(){
Event nextEvent; // type of the next event
Initialize();
MakeTrace(nextEvent);
while(simTime < endTime || numCallsInSystem != 0){
nextEvent = Transition();
switch(nextEvent){
case ARRIVAL:
Arrival();
MakeTrace(nextEvent);
break;
case COMPLETION:
Completion();
MakeTrace(nextEvent);
break;
case PERIODEND:
PeriodEnd();
MakeTrace(nextEvent);
break;
default:
cout<<"NO SUCH EVENT!!!!"<<endl;
} //end Switch
} //End While
UpdateIPA();
// do the neccessary data collection here....
} // End Simulate
/******************************************************************************/
void Initialize(){
//Initialize values
simTime = 0;
numCallsTotal = 0;
numCallsServiced = 0;
timeToPeriodEnd = periodLength;
currentPeriod = 1;
numCallsInSystem = 0;
timeToNextArrival = NHPP(arrivalStream, simTime);
for(int j=0;j<maxServers;j++){
workRemain[j] = 0;
T[j] = 0;
for(int i=0;i<totalPeriods;i++){
delta[j][i] = 0;
}
}
//Initialize Service Level Variables
for(int i=0;i<totalPeriods;i++){
successCounter[i] = 0;
numCalls[i] = 0;
}
//Clear customer queue
//some new fast way to do this?
// delete customer objects that Customer* is pointing to.
for(list<Customer*>::iterator i =queue.begin() ;i!=queue.end();i++){
delete *i;
*i=0;
}
delKs.clear(); // clear vector of derivatives
xiKs.clear();
Ak.clear();
queue.clear(); // clear the list of all pointers
if(!queue.empty()){
cout<<"Queue is not EMPTY!!!"<<endl;
}
} // End Initialize
/******************************************************************************/
Event Transition(){
double timeToNextEvent; //time in minutes from simTime to next event
double timeToNextCompletion; //time to next service completion
Event nextEvent; //type of event to return
serverLeastWork = 0; //holds index of busy server with the least work
char dummy;
//get the index serverLeastWork of the server with the least amount of work
//among the set of busy servers B, serverLeastWork = 0 if all servers free,
//we assume that workRemain[j] = 0 indicates that server j is free, this holds
//with probability 1
for(int i=0;i<numServers[currentPeriod-1];i++){
//if server i is in B, i is a candidate for serverLeastWork
if(workRemain[i] > 0){
if(workRemain[serverLeastWork] == 0) {serverLeastWork = i;}
else if(workRemain[i] < workRemain[serverLeastWork]) {serverLeastWork = i;}
}
}
//cout<<"Server with least work: "<<serverLeastWork<<", Remaining Work: "<<workRemain[serverLeastWork]<< endl;
//cin>>dummy;
//if all servers are free, set timeToNextCompletion to infinity
if(workRemain[serverLeastWork] == 0) {timeToNextCompletion = 999999;}
//otherwise, set timeToNextCompletion
else {timeToNextCompletion = workRemain[serverLeastWork];}
//cout<<"timeToNextCompletion: "<<timeToNextCompletion<<endl;
//cout<<"timeToPeriodEnd: "<<timeToPeriodEnd<<endl;
//cout<<"timeToNextArrival: "<<timeToNextArrival<<endl;
if(timeToNextCompletion<=timeToPeriodEnd && timeToNextCompletion<=timeToNextArrival){
nextEvent = Event(1);
//cout<<"Next Event is a Service Completion"<<endl;
timeToNextEvent = timeToNextCompletion;
}
else if(timeToNextArrival<timeToPeriodEnd && timeToNextArrival<timeToNextCompletion){
nextEvent = Event(0);
//cout<<"Next Event is an Arrival"<<endl;
timeToNextEvent = timeToNextArrival;
}
else{
nextEvent = Event(2);
//cout<<"Next Event is a end of Period"<<endl;
timeToNextEvent = timeToPeriodEnd;
}
//Update simulation clock
simTime += timeToNextEvent;
timeToNextArrival -= timeToNextEvent;
timeToPeriodEnd -= timeToNextEvent;
// Update remaining work in busy servers Step 1.vi
// Update delay at busy servers Step 1.vii
for(int i=0;i<numServers[currentPeriod-1];i++){
//if server i is in B, update remaining work, update the delay
if(workRemain[i] > 0){
workRemain[i] -= timeToNextEvent;
delta[i][currentPeriod-1] -= (timeToNextEvent/serviceRate);
}
}
return nextEvent;
} //End Transition
/******************************************************************************/
void Arrival(){
numCalls[currentPeriod-1]++; //increment no. of calls in this period
numCallsTotal++; //increment no. of incoming calls
Ak.push_back(simTime);
timeToNextArrival = NHPP(arrivalStream, simTime)-simTime; //update time until next call
double newWork = Exponential(workStream, serviceRate);
bool serverIsAssigned = false;
//if queue is empty and there are free servers, then there is a server
//available to answer the call immediately
if(numCallsInSystem < numServers[currentPeriod-1]){
successCounter[currentPeriod-1]++; //this is a successful service
int assignedServer = 0; //index of server assigned to this call
//make default assignedServer to be the first free server
for(int j=0; j<numServers[currentPeriod-1];j++){
if(workRemain[j]==0){
serverIsAssigned = true;
assignedServer=j;
break;
}
}
//assign call to server who has been available for service the longest
//among all free servers
for(int i=assignedServer;i<numServers[currentPeriod-1];i++){
if(workRemain[i] == 0 && T[i] <T[assignedServer]) {
assignedServer = i;
}
}
if(serverIsAssigned!=true){
cout<<"SERVER NOT ASSIGNED. Replication: "<<currentReplication<<", simTime: "<<simTime<<endl;
ErrorTrace(TraceFile);
}
numCallsServiced++; //increment no. of calls that have started service
workRemain[assignedServer] = newWork; //assign service time and mark server busy
//2.b.iv.f: assign value to xiKs = T[assignedServer]
xiKs.push_back(T[assignedServer]);
//2.b.iv.g: assign values to delKs =delta[assignedServer][]
delKs.push_back(delta[assignedServer]);
//2.b.iv.h: reset delta[assignedServer][] to 0
//reset the accumulated delay of assignedServer to zero
for(int i=0;i<totalPeriods; i++){
delta[assignedServer][i] = 0;
}
}
//if there are no free servers to answer call immediately, put
//customer in the queue
else {
Customer *pNewcust= new Customer(simTime, numCallsTotal, newWork, currentPeriod); //create new customer
queue.push_back(pNewcust);//add customer to the queue w/id
}
CheckIPACounters();
numCallsInSystem++; //increment no. calls in sys
} //End Arrival
/******************************************************************************/
void Completion(){
T[serverLeastWork] = simTime; //record the time the server became free
if(numCallsInSystem > numServers[currentPeriod-1]){
//there are calls waiting to be answered, answer the next call in queue
numCallsServiced++; //increment no. of calls that have started service
Customer * pNextCust = *(queue.begin());
if(simTime - pNextCust->attribute <= maxWait) {
successCounter[pNextCust->arrivalPeriod-1]++;
}
workRemain[serverLeastWork] = pNextCust->workload;
queue.pop_front(); //remove customer from queue
delete pNextCust; // memory management. delete object pointer is
pNextCust = 0; // referring to, then reassign pointer to NULL
//assign T[serverLeastWork] to xiKs
xiKs.push_back(T[serverLeastWork]);
//assign delKs
delKs.push_back(delta[serverLeastWork]);
}
/* deleted in version2 of the IPA algorithm
else{
//reset the accumulated delay
for(int i=0;i<totalPeriods; i++){
delta[serverLeastWork][i] = 0;
}
}
*/
CheckIPACounters();
numCallsInSystem--;
} //End Completion
/******************************************************************************/
void PeriodEnd(){
// Case i(more periods remaining)
if(currentPeriod < totalPeriods) { // if 1
if(traceOn) TraceFile<<"\t PeriodEnd-CASEI"<<endl;
//There are more periods
currentPeriod++;
timeToPeriodEnd = periodLength;
if(numServers[currentPeriod-1] < numServers[currentPeriod-2]){
if(traceOn)TraceFile<<"\t PeriodEnd-CASE I- Need to reduce servers"<<endl;
PeriodEndTrace();
//we need to reduce the number of servers
//let l be a vector such that l(i) is the index of the server with
//the ith most remaining work
//reorder T, delta, and workRemain according to l
int index;
double temp;
for(int i=0;i<maxServers-1;i++ ){
index = i;
for(int j=i+1;j<maxServers;j++){
if(workRemain[j]>workRemain[index]) {index = j;}
}
if(index > i){
//swap workRemain values
temp = workRemain[index];
workRemain[index] = workRemain[i];
workRemain[i] = temp;
//swap T values
temp = T[index];
T[index] = T[i];
T[i] = temp;
//swap delta values for each period
for(int k=0;k<totalPeriods;k++){
temp = delta[index][k];
delta[index][k] = delta[i][k];
delta[i][k] = temp;
} // end for
} //end if
} // end for
//clean up remaining service
for(int r=numServers[currentPeriod-1]+1; r<=numServers[currentPeriod-2];r++){
//call extends to next period, effectively one more server
if(workRemain[r-1]>timeToPeriodEnd) {numServers[currentPeriod-1]++;}
else if(workRemain[r-1]>0){
numCallsInSystem--;
workRemain[r-1] = 0;
}
} //end for
PeriodEndTrace();
} //end if
if(numServers[currentPeriod-1] > numServers[currentPeriod-2]){
if(traceOn)TraceFile<<"\t PeriodEnd-CASE I- need to increase Servers"<<endl;
//we need to increase the number of servers
for(int r=numServers[currentPeriod-2]+1; r<=numServers[currentPeriod-1];r++){
T[r-1] = simTime;
for(int i=0;i<totalPeriods;i++){
delta[r-1][i] = 0;
}
//if a call is waiting to be served
if(traceOn)TraceFile<<"\t r: "<<r<<", NumCallsInSystem: "<<numCallsInSystem<<endl;
if(numCallsInSystem >= r){
if(traceOn)TraceFile<<"Allocating waiting customer to incoming Server"<<endl;
numCallsServiced++; //increment no. calls to start service
Customer * pNextCust = *queue.begin();
if(simTime - pNextCust->attribute <= maxWait) {
successCounter[pNextCust->arrivalPeriod-1]++;
}
workRemain[r-1] = pNextCust->workload;
queue.pop_front(); //remove customer from queue
delete pNextCust; // deleting Customer object
pNextCust = 0; // assigning pointer to NULL
//assign values to xiKs
xiKs.push_back(T[r-1]);
//assign values to delKs
delKs.push_back(delta[r-1]);
} //end if
} // end for
} //end of if 3
} //End 1st If
//Case ii (end of last period)
else if(currentPeriod ==totalPeriods) {
if(traceOn)TraceFile<<"\t PeriodEnd-CASE II"<<endl;
//we have reached the end of the last period
currentPeriod++;
timeToPeriodEnd = maxWait;
numServers[currentPeriod-1] = numServers[currentPeriod-2];
}
//Case iii (time to terminate)
else if(currentPeriod>totalPeriods){
if(traceOn)TraceFile<<"\t PeriodEnd-CASE III"<<endl;
//then we terminate the simulation
//currentPeriod > totalPeriods implies simTime > endTime
Termination();
numCallsInSystem = 0; //this will terminate the simulation
}
}//End PeriodEnd
/******************************************************************************/
void Termination(){
vector<double> temp(totalPeriods,0); //create a 1 by p vector of zeros
if(!queue.empty()){
//queue is not empty which implies that
//no servers available. Set xiKs=infinity and DelKs=0 for all remaining
//customers in the queue
//Update IPA counters for each customer in the queue, and for an imaginary 'next customer'
for(int i = 0;i<queue.size()+1;i++){
xiKs.push_back(99999);
delKs.push_back(temp);// add a 1xp zero vector
numCallsServiced++;
}
}
else {
//assign the IPA counters
//for the imaginary 'next customer', the (NumCallsTotal+1)th customer
bool serverIsAssigned = false;
int assignedServer = 0; //index of server assigned to this call
//if queue is empty and there are free servers, then there is a server
//available to answer the next call immediately
//make default assignedServer to be the first free server
for(int j=0; j<numServers[totalPeriods-1];j++){
if(workRemain[j]==0){
serverIsAssigned = true;
assignedServer=j;
break;
}
}
//assign call to server who has been available for service the longest
//among all free servers
for(int i=assignedServer;i<numServers[currentPeriod-1];i++){
if(workRemain[i] == 0 && T[i] <T[assignedServer]) {
assignedServer = i;
}
}
if(serverIsAssigned!=true){
//the case when the queue is empty, and all servers are busy
xiKs.push_back(99999);
delKs.push_back(temp);// add a 1xp zero vector
}
else{
//there is a server available to service the next 'imaginary' incoming call
xiKs.push_back(T[assignedServer]);
delKs.push_back(delta[assignedServer]);
}
numCallsServiced++;
}
CheckIPACounters();
if(numCallsServiced!=(numCallsTotal+1)){
cout<<"ACCOUNTING ERROR in Termination of replication "<<currentReplication<<endl;
}
}// End Termination()
/******************************************************************************/
void UpdateIPA(){
//update the mean IPA gradients
vector<double> zeroVect(totalPeriods,0);
//dJik_dMu is the IPA gradient obtained from this particular replication
vector< vector<double> > dJik_dMu(totalPeriods,zeroVect); // this creates a totalPeriods by totalPeriods matrix of zeros
double dFdU = 0;
double A;
//Sum up Jik for this replication from k=1 to k=numCallsTotal+1=NumCallsServiced
CheckIPACounters();
for(int k=1;k<=numCallsServiced;k++){
for(int i=1; i<=totalPeriods;i++){
for(int j = 1; j<=totalPeriods;j++){
if( (periodLength*(i-1)+maxWait)< xiKs[k-1] && xiKs[k-1]<(periodLength*(i)+maxWait)){
if(k==1){A = 0;}
else{ A = Ak[k-2]; }
dJik_dMu[i-1][j-1] += -exp(-Integral(A,xiKs[k-1]-maxWait))* delKs[k-1][j-1];
}
else{
//do nothing. this is equivalent to adding zero.
}
}//end looping through j periods
} //end looping through i periods
}//end looping through each call
//Update the mean and variance of the IPA gradient with the results
//of this replication
//One pass variance/mean calculation algorithms obtained from
//http://mathworld.wolfram.com/SampleVarianceComputation.html
double OldMean;
for(int l=0;l<totalPeriods;l++){
for(int j=0; j<totalPeriods;j++){
if(currentReplication==1){
meanIPAGrad[l][j]=dJik_dMu[l][j]; //used for 1st moment
}
else {
OldMean = meanIPAGrad[l][j];
meanIPAGrad[l][j]=((currentReplication-1)*OldMean+dJik_dMu[l][j])/(currentReplication); //Sample mean after (repnumber+1)th replication
varIPAGrad[l][j] = (1.00-(1.00/(currentReplication-1)))*varIPAGrad[l][j] + (currentReplication)*pow((meanIPAGrad[l][j]-OldMean),2);
}
}//End loop through j periods
}//end loop through k periods
}// End UpdateIPA()
/******************************************************************************/
double Integral(double a, double b){
if (a>b){
double temp = b;
b = a;
a= temp;
}
int periodA=int(ceil(a / periodLength));
int periodB=int(ceil(b / periodLength));
double fA = ((a-(periodA-1)*periodLength)*(arrivalRate[periodA]-arrivalRate[periodA-1])/periodLength)+arrivalRate[periodA-1];
double fB = ((b-(periodB-1)*periodLength)*(arrivalRate[periodB]-arrivalRate[periodB-1])/periodLength)+arrivalRate[periodB-1];
double area = 0;
if(periodA == periodB){
area = 0.5*(b-a)*(fA+fB);
}
else{
area = 0.5 * ( (periodA*periodLength-a)*(arrivalRate[periodA]+fA)+ (b-(periodB-1)*periodLength)*(arrivalRate[periodB-1]+fB)) ;
if(periodB-periodA>1){
for(int i = periodA+1;i<periodB;i++){
area += periodLength*0.5*(arrivalRate[i-1]+arrivalRate[i]);
} //end for
} //end if
} //end else
return area;
}//end of Integral
/******************************************************************************/
double NHPP(RngStream &generator, double time){
//this function returns the time (from start of simulation) of the next arrival
//"time" is the current simulation time
bool realArrival = false; // has a real arrival been generated from thinning?
double nextTime = time; // the next possible arrival time, used in thinning
double currentArrivalRate; // actual arrival rate at nextTime
double thinning; // U(0,maxArrivalRate), for thinning
//thinning algorithm
while(!realArrival){
//generate an arrival time, based on maxArrival rate
nextTime+= -(1.0/maxArrivalRate)*log(1-generator.RandU01());
//TraceFile<<"NHPP-nextTime: "<<nextTime<<endl;
// generate thinning variable
thinning = maxArrivalRate*generator.RandU01();
// Determine if the thinning variable is less than current arrival rate
// if so, accept the current nextTime
currentArrivalRate = RateFunct(nextTime);
if(thinning<=currentArrivalRate){
realArrival = true;
}
}
// return next arrival time
if(nextTime>endTime){
nextTime = 999999;
}
return nextTime;
} //End NHPP
/******************************************************************************/
// returns the rate function from a piecewise linear rate function at any given time
double RateFunct(double time){
int periodA=int(ceil(time / periodLength)); //period that 'time' falls into
if(periodA>totalPeriods){
return 999999;
}
//returns the arrival rate at this particular time
return ((time-(periodA-1)*periodLength)*(arrivalRate[periodA]-arrivalRate[periodA-1])/periodLength)+arrivalRate[periodA-1];
}//End RateFunct
/******************************************************************************/
//returns an exp(lambda) r.v. using a specified stream
double Exponential(RngStream &generator, double lambda){
double variable;
variable = -(1.0/lambda)*log(1-generator.RandU01());
return variable;
}//End Exponential
/******************************************************************************/
void ShowParameters(){
cout<<"Number of Periods: "<<totalPeriods<<endl;
cout<<"Number of Replications: "<<totalReplications<<endl;
cout<<"ServiceRate(services/min): "<<serviceRate<<endl;
cout<<"Length of Period: "<<periodLength<<endl;
cout<<"Maximum wait time(minutes): "<<maxWait<<endl;
cout<<"Maximum number of Servers: "<<maxServers<<endl;
cout<<"Number of Servers in each period:"<<endl;
for(int i=0; i<totalPeriods;i++){
cout<<i+1<<" " <<numServers[i]<<endl;
}
cout<<"Max. Arrival Rate: "<<maxArrivalRate<<endl;
cout<<"Arrival Rate at the end of each period:"<<endl;
for(int i=0; i<totalPeriods+1;i++){
cout<<i<<" " <<arrivalRate[i]<<endl;
}
} //End ShowParameters()
/******************************************************************************/
void WriteStatsFile(){
fstream StatsFile("Stats.txt",ios::in | ios::out |ios::trunc);
for(int i = 0;i<totalPeriods;i++){
StatsFile<<meanSuccess[i]<<endl;
}
for(int i = 0;i<totalPeriods;i++){
StatsFile<<varSuccess[i]<<endl;
}
for(int i = 0;i<totalPeriods;i++){
StatsFile<<meanNumCalls[i]<<endl;
}
StatsFile.close();
fstream GradFile("Gradient.txt",ios::in | ios::out |ios::trunc);
for(int i=0;i<totalPeriods;i++){
for(int j=0;j<totalPeriods;j++){
GradFile<<meanIPAGrad[i][j]<<"\t";
}
GradFile<<endl;
}
GradFile.close();
}//End WriteStatsFile()
/******************************************************************************/
void MakeTrace(Event event){
if(traceOn==false||currentReplication!=2) return;
if (!TraceFile.is_open()) { // check for successful opening
cout << "Trace.txt file could not be opened!" << endl;
}
TraceFile<<"********************************************************"<<endl;
TraceFile<<"\t Next EventType: "<<event<<endl;
TraceFile<<"\t Current Replication: "<<currentReplication<<endl;
TraceFile<<"\t Current Period: "<<currentPeriod<<endl;
TraceFile<<"\t Total Periods: "<<totalPeriods<<endl;
TraceFile<<"\t Current Simulation Time: "<<simTime<<endl;
TraceFile<<"\t delKs.size()="<<delKs.size()<<endl;
TraceFile<<"\t xiKs.size()="<<xiKs.size()<<endl;
TraceFile<<"\t Ak.size()="<<Ak.size()<<endl;
TraceFile<<"\t Number of Calls in Total: "<<numCallsTotal<<endl;
TraceFile<<"\t Number of Calls that have started Service: "<<numCallsServiced<<endl;
TraceFile<<"\t Number of Calls in System : "<<numCallsInSystem<<endl;
TraceFile<<"\t Number of Servers in this Period: "<<numServers[currentPeriod-1]<<endl;
TraceFile<<"\t Number of Customers in the queue: "<<queue.size()<<endl;
TraceFile<<"\t Time to end of Period: "<<timeToPeriodEnd<<endl;
TraceFile<<"\t Time to next arrival: "<<timeToNextArrival<<endl;
TraceFile<<endl;
TraceFile<<"\t Work Remaining at each Server: "<<endl;
for(int i=0; i<maxServers;i++){
TraceFile<<"\t "<< i+1 <<"\t"<< workRemain[i]<<endl;
}
TraceFile<<endl;
TraceFile<<"\t Time server became free:"<<endl;
for(int i =0; i<maxServers;i++){
TraceFile<<"\t " <<i+1<<"\t"<<T[i]<<endl;
}
}//End MakeTrace()
/******************************************************************************/
void PeriodEndTrace(){
if(traceOn==false||currentReplication!=2) return;
if (!TraceFile.is_open()) { // check for successful opening
cout << "Trace.txt file could not be opened!" << endl;
}
TraceFile<<"*******Tracing Vector ReOrdering***********"<<endl;
for(int i=0; i<maxServers;i++){
TraceFile<<"\t "<< i+1 <<"\t"<< workRemain[i]<<"\t"<<T[i]<<endl;
}
for(int i=0;i<maxServers;i++){
TraceFile<<"\t "<< i+1 ;
for(int j=0;j<totalPeriods;j++){
TraceFile<<"\t "<<delta[i][j];
}
TraceFile<<endl;
}
} //End PeriodEndTrace()
/*******************************************************************************/
void ErrorTrace(fstream &DataFile){
if (!TraceFile.is_open()) { // check for successful opening
cout << "ErrorTrace file could not be opened!" << endl;
return;
}
DataFile<<"****************ERROR TRACE START************************"<<endl;
DataFile<<"\t Current Replication: "<<currentReplication<<endl;
DataFile<<"\t Current Period: "<<currentPeriod<<endl;
DataFile<<"\t Total Periods: "<<totalPeriods<<endl;
DataFile<<"\t Current Simulation Time: "<<simTime<<endl;
DataFile<<"\t Number of Calls in Total: "<<numCallsTotal<<endl;
DataFile<<"\t Number of Calls that have started Service: "<<numCallsServiced<<endl;
DataFile<<"\t Number of Calls in System : "<<numCallsInSystem<<endl;
DataFile<<"\t Number of Servers in this Period: "<<numServers[currentPeriod-1]<<endl;
DataFile<<"\t Number of Customers in the queue: "<<queue.size()<<endl;
DataFile<<"\t Time to end of Period: "<<timeToPeriodEnd<<endl;
DataFile<<"\t Time to next arrival: "<<timeToNextArrival<<endl;
DataFile<<endl;
DataFile<<"\t Work Remaining at each Server: "<<endl;
for(int i=0; i<maxServers;i++){
DataFile<<"\t "<< i+1 <<"\t"<< workRemain[i]<<endl;
}
DataFile<<endl;
DataFile<<"\t Time server became free:"<<endl;
for(int i =0; i<maxServers;i++){
DataFile<<"\t " <<i+1<<"\t"<<T[i]<<endl;
}
DataFile<<"****************ERROR TRACE END**************************"<<endl;
}//End ErrorTrace()
/******************************************************************************/