**Meet the editor**

Dr Eldin Wee Chuan Lim is a Lecturer in the Department of Chemical & Biomolecular Engineering at the National University of Singapore (NUS). He is a Member of the Editorial Boards for the Elsevier journal, Advanced Powder Technology, the American Journal of Fluid Dynamics published by Scientific & Academic Publishing, and the Journal of Chemical Engineering

& Process Technology published by OMICS Publishing Group. Some of the academic awards that he has received include the HPC Quest Silver Award presented by IBM and the Institute of High Performance Computing Singapore, ERS-IASC Young Researcher Award presented by the European Regional Section of the International Association for Statistical Computing, two Best Paper Awards presented at the 2010 SECAU Security Congress and the Frontier Award presented by the Fluid & Particle Processing Division of The Society of Chemical Engineers, Japan. He is also a recipient of the Teaching Commendation award presented by the Faculty of Engineering of NUS.

Contents

**Preface IX** 

Chapter 2 **The Speedup of Discrete Event** 

Chapter 3 **Sensitivity Analysis in Discrete** 

Wennai Wang and Yi Yang

and Jonathan Daniel Friend

**Section 2 Novel Integration of Discrete Event** 

**Section 3 Applications of Discrete Event** 

Chapter 6 **Discrete-Event Simulation of** 

and Andrey Shorov

Chapter 5 **Human Evacuation Modeling 135** 

**Section 1 Fundamental Development and Analyses** 

**of the Discrete Event Simulation Method 1** 

Giulia Pedrielli, Tullio Tolio, Walter Terkaj and Marco Sacco

**Event Simulation Using Design of Experiments 63**  José Arnaldo Barra Montevechi, Rafael de Carvalho Miranda

**Simulation with Other Modeling Techniques 103** 

**Criteria Decision Analysis as a Decision Support Methodology in Complex Logistics Systems 105**  Thiago Barros Brito, Rodolfo Celestino dos Santos Silva, Edson Felipe Capovilla Trevisan and Rui Carlos Botter

**Simulation Towards Various Systems 133** 

Stephen Wee Hun Lim and Eldin Wee Chuan Lim

**Botnet Protection Mechanisms 143**  Igor Kotenko, Alexey Konovalov

Chapter 4 **Discrete Event Simulation Combined with Multiple** 

Chapter 1 **Distributed Modeling of Discrete Event Systems 3** 

**Simulations by Utilizing CPU Caching 47** 

## Contents

#### **Preface XI**


X Contents

Chapter 7 **Using Discrete Event Simulation for Evaluating Engineering Change Management Decisions 169**  Weilin Li

## Preface

With rapid advancements in computing power, computer modeling and simulations have become an important complement to experimentations in many areas of research as well as industrial applications. The Discrete Event Simulation (DES) method has received widespread attention and acceptance by both researchers and practitioners in recent years. The range of application of DES spans across many different disciplines and research fields. In research, further development and advancements of the basic DES algorithm continue to be sought while various hybrid methods derived by combining DES with other simulation techniques continue to be developed. This book presents state-of-the-art contributions on fundamental development of the DES method, novel integration of the method with other modeling techniques as well as applications towards simulating and analyzing the performances of various types of systems. This book will be of interest to undergraduate and graduate students, researchers as well as professionals who are actively engaged in DES related work.

There are nine chapters in this book that are organized into three sections. The first section comprises three chapters that report recent studies on fundamental development and analyses of the DES method. In Chapter 1, Pedrielli and co-authors introduce a distributed modeling approach that allows complex discrete event systems that would otherwise not be practicable to model using conventional simulation techniques to be modeled efficiently. Wang and Yang discuss in Chapter 2 various approaches for fast event scheduling for simulations of large-scale networks. They report the results of computational experiments that demonstrate the performance of a cache aware algorithm to be better than that of a conventional Calendar Queue. In Chapter 3, Montevechi and co-authors present the application of factorial design statistical techniques for identifying significant variables in discrete event simulation models with a view towards speeding up simulation optimization processes.

The approach of integrating DES with various modeling techniques has also attracted the interests of several researchers throughout the world in recent years. In the second section of this book, two chapters on work conducted in this area are presented. Ortega discusses in Chapter 4 a simulation platform that combines DES with stochastic simulation and multi-agent systems for modeling holonic manufacturing systems. Brito describes in Chapter 5 a decision support system that was developed by combining DES with Multiple Criteria Decision Analysis. The application of such a

#### XII Preface

hybrid decision support system towards analysis of a steel manufacturing plant is illustrated.

The final section of this book is devoted to contributions reporting applications of DES towards various systems. Lim and Lim describe simulations of human evacuation processes using a discrete approach for modeling individual human subjects in Chapter 6. Kotenko and co-authors report the development of a DES based environment for analyses of botnets and evaluation of defense mechanisms against botnet attacks in Chapter 7. In Chapter 8, Li proposes a comprehensive DES model that is able to capture the various complexities associated with new product development projects as well as take into account engineering changes that arise stochastically during the course of such projects. In the final chapter of this book, Klug focuses on project management issues that are of relevance to simulation projects in general.

This book represents the concerted efforts of many individuals. First and foremost, I would like to take this opportunity to acknowledge the efforts of all authors who have contributed to the success of this book project. I would also like to thank the support provided by Ms. Mirna Cvijic of InTech Open Access Publisher, without which the publication of this book would not have been possible. Last but certainly not least, and on behalf of all contributing authors, I wish to express my sincere appreciation and gratitude towards InTech Open Access Publisher for transforming this book project from inception to reality.

> **Eldin Wee Chuan Lim**  Department of Chemical & Biomolecular Engineering, National University of Singapore Singapore

**Fundamental Development and Analyses of the Discrete Event Simulation Method** 

**Chapter 1** 

© 2012 Pedrielli et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution,

**Distributed Modeling of Discrete Event Systems** 

Computer simulation is widely used to support the design of any kind of complex system and to create computer-generated "virtual worlds" where humans and/or physical devices are embedded (e.g. aircraft flight simulators [20]). However, both the generation of simulation models and the execution of simulations can be time and cost expensive. While there are already several ways to increase the speed of a simulation run, the scientific challenge for the simulation of complex systems still resides in the ability to model

A computer simulation is a computation that emulates the behavior of some real or

 *Continuous simulation*. Given the discrete nature of the key parameters of a digital computer, including the number of memory locations, the data structures, and the data representation, continuous simulation may be best approximated on a digital computer through time-based discrete simulation where the time steps are sufficiently small

 *Time-based discrete simulation*. In this case the universal time is organized into a discrete set of monotonically increasing timesteps where the choice of the duration of the timestep interval changes as a result of the external stimuli, any change between two subsequent timesteps must occur atomically within the corresponding timestep interval. Regardless of whether its state incurs and changes, a process and all its

 *Discrete event simulation* [5]. The difference between discrete event simulation and timebased simulation is twofold. Firstly, the process being modeled is understood to advance through events under discrete event conditions. Second, an event (i.e. an activity of the process as determined by the model developer) carries with it the potential for affecting the state of the model and is not necessarily related to the

and reproduction in any medium, provided the original work is properly cited.

conceptual systems over time. There are three main simulation techniques [23]:

Giulia Pedrielli, Tullio Tolio, Walter Terkaj and Marco Sacco

Additional information is available at the end of the chapter

(simulate) those systems in a parallel/distributed way [35].

relative to the process being modeled.

parameters may be examined at every time step.

http://dx.doi.org/10.5772/50350

**1. Introduction** 

## **Distributed Modeling of Discrete Event Systems**

Giulia Pedrielli, Tullio Tolio, Walter Terkaj and Marco Sacco

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50350

## **1. Introduction**

Computer simulation is widely used to support the design of any kind of complex system and to create computer-generated "virtual worlds" where humans and/or physical devices are embedded (e.g. aircraft flight simulators [20]). However, both the generation of simulation models and the execution of simulations can be time and cost expensive. While there are already several ways to increase the speed of a simulation run, the scientific challenge for the simulation of complex systems still resides in the ability to model (simulate) those systems in a parallel/distributed way [35].

A computer simulation is a computation that emulates the behavior of some real or conceptual systems over time. There are three main simulation techniques [23]:


progress of time. In this case, the executable model must necessarily be run corresponding to every event to accurately reflect the reality of the process.

Distributed Modeling of Discrete Event Systems 5

3. Asynchronous interactions between the entities 4. Entities which concur for the use of shared resources

the owners of each subsystem do not want to share the information.

simulation (DS) approach which will be the focus of this chapter.

The simulation of complex systems through the use of traditional simulation tools presents several drawbacks, e.g. the long time required to develop the unique monolithic simulation model, the computational effort required for running the simulation, the impossibility to run the simulation model on a set of geographically distributed computers, the absence of fault tolerance (i.e. the work done is lost if one processor goes down), the impossibility to realize a realistic model of the entire system in the case several subsystems are included and

Most of the aforementioned problems can be effectively addressed by the distributed

The chapter will be organized as follows: Section 2 presents the main concepts and definitions together with a literature review on applications and open issues related to distributed simulation. Section 3 delves into the High Level Architecture [1], i.e. the reference standard supporting the distributed simulation. Section 4 shows an application of distributed simulation on a real industrial case in the manufacturing domain [77]. Finally, Section 5 presents the conclusions and the main topics for future research in the field of

Traditional stand alone simulation is based on a simulation clock and an event list. The interaction of the event list and the simulation clock generates the sequence of the events

The execution of any event might cause an update of the value of the state variables, a modification to the event list and (or) the collection of the statistics. Each event is executed

The idea underlying the distributed simulation is to minimize the sequential aspect of traditional simulation. Distributed simulation can be classified into two major categories: (1)

Parallel and distributed computing refers to technologies that enable a simulation program to be executed on a computing system containing multiple processors, such as personal

The main benefits resulting from the adoption of distributed computing technologies are

 *Reduced execution time*. By decomposing a large simulation computation into many subcomputations and executing the sub-computations concurrently across different

based on the simulation time assigned to it, i.e. the simulation is sequential.

parallel and distributed computing, and (2) distributed modeling.

computers, interconnected by a communication network [20].

processors, one can reduce the global execution time.

5. Connectivity between the entities

distributed simulation.

that have to be simulated.

[20]:

**2. Distributed simulation** 

Since continuous simulation is simply academic and cannot be reproduced on real computers, it is important to comment the difference between time-based simulation and discrete event simulation.

Under time-based simulation, the duration of the timestep interval is determined based on the nature of the specific activity or activities of the process that the model developer considers important and worth modeling and simulating. Similarly, under discrete event simulation, events for a given process are also identified on the basis of the activity or activities the model developer views as important. Whereas time-based simulation constitutes the logical choice for processes in which the activity is distributed over every timestep, discrete event simulation is more efficient when the activity of a process being modeled is sparsely distributed over time. The overhead in discrete event simulation, arising from the additional need to detect and record the events, is higher than in the simpler time-based technique and must be more than compensated by the savings not to have to execute the model at every time step.

A fundamental difference between time-based and discrete event simulations lies in their relationship to the principle of causality. In the time-based approach, while a cause may refer to a process state at a specific timestep, the fact that the state of the process is observed at every subsequent time step reflects the assumption that the effect of the cause is expected. Thus both the cause and the effect refer to the observed state of the process in time-based simulation. In discrete event simulation, both cause and effects refer to events. However, upon execution due to an event, a model may not generate an output event thus appearing to imply that a cause will not necessary be accompanied by corresponding observed facts.

Discrete Event Simulation (DES) has been widely adopted to support system analysis, education and training, organizational change [43] in a range of diverse areas such as commerce [13], manufacturing ([14],[38], [79]), supply chains [24], health services and biomedicine ([3], [18]), simulation in environmental and ecological systems [6], city planning and engineering [45], aerospace vehicle and air traffic simulation [40], business administration and management [16], military applications [17].

All the aforementioned areas are usually characterized by the presence of complex systems. Indeed, a system represented by a simulation model is defined as complex when it is extremely large, i.e. a large number of components characterize it, or a large number of interactions describes the relationships between objects within the system, or it is geographically dispersed. In all cases the dynamics can be hard to describe. The complexity is reflected in the system simulation model that can be characterized according to the following concepts [23]:


have to execute the model at every time step.

discrete event simulation.

following concepts [23]:

simulation [62]

2. Asynchronous behavior of the entities

progress of time. In this case, the executable model must necessarily be run

Since continuous simulation is simply academic and cannot be reproduced on real computers, it is important to comment the difference between time-based simulation and

Under time-based simulation, the duration of the timestep interval is determined based on the nature of the specific activity or activities of the process that the model developer considers important and worth modeling and simulating. Similarly, under discrete event simulation, events for a given process are also identified on the basis of the activity or activities the model developer views as important. Whereas time-based simulation constitutes the logical choice for processes in which the activity is distributed over every timestep, discrete event simulation is more efficient when the activity of a process being modeled is sparsely distributed over time. The overhead in discrete event simulation, arising from the additional need to detect and record the events, is higher than in the simpler time-based technique and must be more than compensated by the savings not to

A fundamental difference between time-based and discrete event simulations lies in their relationship to the principle of causality. In the time-based approach, while a cause may refer to a process state at a specific timestep, the fact that the state of the process is observed at every subsequent time step reflects the assumption that the effect of the cause is expected. Thus both the cause and the effect refer to the observed state of the process in time-based simulation. In discrete event simulation, both cause and effects refer to events. However, upon execution due to an event, a model may not generate an output event thus appearing to imply that a cause will not necessary be accompanied by corresponding observed facts.

Discrete Event Simulation (DES) has been widely adopted to support system analysis, education and training, organizational change [43] in a range of diverse areas such as commerce [13], manufacturing ([14],[38], [79]), supply chains [24], health services and biomedicine ([3], [18]), simulation in environmental and ecological systems [6], city planning and engineering [45], aerospace vehicle and air traffic simulation [40], business

All the aforementioned areas are usually characterized by the presence of complex systems. Indeed, a system represented by a simulation model is defined as complex when it is extremely large, i.e. a large number of components characterize it, or a large number of interactions describes the relationships between objects within the system, or it is geographically dispersed. In all cases the dynamics can be hard to describe. The complexity is reflected in the system simulation model that can be characterized according to the

1. Presence of entity elements that are dynamically created and moved during a

administration and management [16], military applications [17].

corresponding to every event to accurately reflect the reality of the process.

The simulation of complex systems through the use of traditional simulation tools presents several drawbacks, e.g. the long time required to develop the unique monolithic simulation model, the computational effort required for running the simulation, the impossibility to run the simulation model on a set of geographically distributed computers, the absence of fault tolerance (i.e. the work done is lost if one processor goes down), the impossibility to realize a realistic model of the entire system in the case several subsystems are included and the owners of each subsystem do not want to share the information.

Most of the aforementioned problems can be effectively addressed by the distributed simulation (DS) approach which will be the focus of this chapter.

The chapter will be organized as follows: Section 2 presents the main concepts and definitions together with a literature review on applications and open issues related to distributed simulation. Section 3 delves into the High Level Architecture [1], i.e. the reference standard supporting the distributed simulation. Section 4 shows an application of distributed simulation on a real industrial case in the manufacturing domain [77]. Finally, Section 5 presents the conclusions and the main topics for future research in the field of distributed simulation.

#### **2. Distributed simulation**

Traditional stand alone simulation is based on a simulation clock and an event list. The interaction of the event list and the simulation clock generates the sequence of the events that have to be simulated.

The execution of any event might cause an update of the value of the state variables, a modification to the event list and (or) the collection of the statistics. Each event is executed based on the simulation time assigned to it, i.e. the simulation is sequential.

The idea underlying the distributed simulation is to minimize the sequential aspect of traditional simulation. Distributed simulation can be classified into two major categories: (1) parallel and distributed computing, and (2) distributed modeling.

Parallel and distributed computing refers to technologies that enable a simulation program to be executed on a computing system containing multiple processors, such as personal computers, interconnected by a communication network [20].

The main benefits resulting from the adoption of distributed computing technologies are [20]:

 *Reduced execution time*. By decomposing a large simulation computation into many subcomputations and executing the sub-computations concurrently across different processors, one can reduce the global execution time.

	- *Geographical distribution*. Executing the simulation program on a set of geographically distributed computers enables one to create virtual worlds with multiple participants that are physically located at different sites.

Distributed Modeling of Discrete Event Systems 7

**2.1. HLA-standard: An overview** 

the DoD):

simulation)

**Figure 1.** HLA Reference Architecture

**Figure 2.** RTIAmbassador and FederateAmbassador

HLA (IEEE standard 1516) is a software architecture designed to promote the use and interoperation of simulators. HLA was based on the premise that no single simulator could satisfy all uses and applications in the defense industry and it aimed at reducing the time

The HLA architecture (Figure1) defines a *Federation* as a collection of interacting simulators (*federates*), whose communication is orchestrated by a *Runtime Infrastructure* (RTI) and an interface. Federates can be either simulations, surrogates for live players, or tools for distributed simulation. They are defined as having a single point of attachment to the RTI

HLA can combine the following types of simulators (following the taxonomy developed by

*Constructive* - simulated people operating simulated systems (e.g. a discrete event

and cost required to create a synthetic environment for a new purpose.

*Live* - real people operating real systems (e.g. a field test)

and might consist of several processes, perhaps running on several computers.

*Virtual* - real people operating simulated systems (e.g. flight simulations)


The definition of distributed modeling can be given by highlighting the differences compared to the concept of parallel and distributed computing as presented by Fujimoto [20]. If a single simulator is developed and the simulation is executed on multiple processors we talk about *parallel* and *distributed computing*. Whereas if several simulators are combined into a distributed architecture we talk about *distributed modeling*; in this case, the simulation execution requires the synchronization between the different simulators.

The *distributed computing* can be still applied to each simulator in a distributed simulation model [60], but the complexity related to the synchronization of the different models can be such that the performance of the simulation (in terms of speed) can be worse than when a single simulation model is developed. This drawback related to the decrease in the efficiency in terms of speed of simulation leads to the following question: "Why is it useful to develop a distributed simulation model?". The following benefits represent an answer to this question ([57], [77]):


The feasibility of the distributed simulation concept was demonstrated by the SIMNET project (SIMulator NETworking [73]), which ran from 1983 to 1990. As consequence of this project, a set of protocols were developed for interconnecting simulations and the Distributed Interactive Simulation (DIS) standard was the first one. Afterwards, the High Level Architecture (HLA) standard ([1], [15], [27]) was developed by the U.S. Department of Defense (DoD) under the leadership of the Defense Modeling and Simulation Office (DMSO). The next sub-section presents a general overview of the HLA standard for distributed simulation, whereas Section 2.2 gives an overview of distributed simulation in civilian applications.

#### **2.1. HLA-standard: An overview**

6 Discrete Event Simulations – Development and Applications

failure.

this question ([57], [77]):

easier to cope with.

civilian applications.

that are physically located at different sites.

 *Geographical distribution*. Executing the simulation program on a set of geographically distributed computers enables one to create virtual worlds with multiple participants

 *Fault tolerance*. If one processor goes down, it may be possible for other processors to pick up the work of the failed machine allowing the simulation to proceed despite the

The definition of distributed modeling can be given by highlighting the differences compared to the concept of parallel and distributed computing as presented by Fujimoto [20]. If a single simulator is developed and the simulation is executed on multiple processors we talk about *parallel* and *distributed computing*. Whereas if several simulators are combined into a distributed architecture we talk about *distributed modeling*; in this case, the simulation

The *distributed computing* can be still applied to each simulator in a distributed simulation model [60], but the complexity related to the synchronization of the different models can be such that the performance of the simulation (in terms of speed) can be worse than when a single simulation model is developed. This drawback related to the decrease in the efficiency in terms of speed of simulation leads to the following question: "Why is it useful to develop a distributed simulation model?". The following benefits represent an answer to

 *Complexity management*. If the complexity of the system to be simulated grows and the modeling of each sub-system requires various and specific expertise, then the realization of a single monolithic simulation model is not feasible [65]. Under the distributed modeling approach the problem is decomposed in several sub-problems

 *Overcoming the lack of shared information*. The developer of a simulation model can hardly access all the information characterizing the whole system to model, again hindering

 *Reusability*. The development of a simulation model always represents a costly activity, thus the distributed modeling can be seen as a possibility to integrate pre-existing

The feasibility of the distributed simulation concept was demonstrated by the SIMNET project (SIMulator NETworking [73]), which ran from 1983 to 1990. As consequence of this project, a set of protocols were developed for interconnecting simulations and the Distributed Interactive Simulation (DIS) standard was the first one. Afterwards, the High Level Architecture (HLA) standard ([1], [15], [27]) was developed by the U.S. Department of Defense (DoD) under the leadership of the Defense Modeling and Simulation Office (DMSO). The next sub-section presents a general overview of the HLA standard for distributed simulation, whereas Section 2.2 gives an overview of distributed simulation in

the feasibility of developing a unique and monolithic simulation model.

simulators and to avoid the realization of new models.

*Integration* of simulators that execute on machines from different manufacturers.

execution requires the synchronization between the different simulators.

HLA (IEEE standard 1516) is a software architecture designed to promote the use and interoperation of simulators. HLA was based on the premise that no single simulator could satisfy all uses and applications in the defense industry and it aimed at reducing the time and cost required to create a synthetic environment for a new purpose.

The HLA architecture (Figure1) defines a *Federation* as a collection of interacting simulators (*federates*), whose communication is orchestrated by a *Runtime Infrastructure* (RTI) and an interface. Federates can be either simulations, surrogates for live players, or tools for distributed simulation. They are defined as having a single point of attachment to the RTI and might consist of several processes, perhaps running on several computers.

HLA can combine the following types of simulators (following the taxonomy developed by the DoD):


**Figure 1.** HLA Reference Architecture

**Figure 2.** RTIAmbassador and FederateAmbassador

The HLA standard provides four main components for the realization and management of a federation:

Distributed Modeling of Discrete Event Systems 9

Data Distribution Management. These services can reduce unnecessary information

HLA overcame the shortcomings of the DIS standard by being simulation-domain neutral (it was not developed referred to any specific language, therefore HLA provides means to describe any data exchange format as required and specifying functionalities for time

HLA provides Application Programming Interfaces (APIs) for all the classes of services just mentioned, but the RTI software and algorithms are not defined by HLA. Also the operations in the *FederateAmbassador* need to be implemented at the federate level, as part of

These facts have caused the growth of multiple HLA-RTI implementations (e.g. [80], [81]) and the development of *ad-hoc* solutions for the adapters on the federate side [25]. In particular the last aspect represents one of the most relevant criticalities in applying HLA for distributed simulation: the lack of a standardized approach to adapt a simulator within an HLA-based distributed architecture, makes a distributed simulation project time expensive since a lot of implementation is required in addition to the effort to build the simulation

This consideration represents one of the leading arguments for the research community in the direction of the development of additional complementary standards (Section 3) to ease

It is the objective of the next section to analyze the state of the art on the adoption and

Herein the attention is focused on distributed modeling of complex systems in civilian

HLA constitutes an enabler for implementing the distributed simulation. The standard, though, was conceived for military applications and several problems arise when trying to interoperate heterogeneous simulators in civilian applications (the terminology Commercial off-the-shelf discrete-event simulation packages CSPs [62] will be used to describe

Boer [12] investigated the main benefits and criticalities related to the industrial application of HLA by interviewing the actors involved in the problem (e.g. simulation model developers, software houses, HLA experts, [9]-[11]). The results of the survey showed that CSPs vendors do not see direct benefits in using distributed simulation, whereas in industry HLA is considered troublesome because of the lack of experienced users and the complexity of the standard. In addition, as suggested in [49], although the approaches and general methods used in military and civilian simulation communities have similarities, the

the creation and management of an HLA-based distributed simulation.

advancements in the use of HLA-based distributed simulation technique.

commercially available simulators for the analysis of Discrete Event Systems).

**2.2. Distributed simulation in civilian applications** 

transfer between federates by filtering out irrelevant data.

management and bandwidth control (see the FIS module).

the federate code or some interface service (*adapter*).

model.

domain.


The federates cannot directly exchange information throughout the federation, instead the RTI plays the role of the operating system of the distributed simulation, providing a set of general-purpose services for federation management and enabling the federates in carrying out federate-to-federate *interactions*. In particular interactions represent an explicit action taken by a federate that may have some effect on another federate within a federation execution, such action can be tied with a specific time defined as *interactionTime*, when the action takes place.

Each federate is endowed with an *RTIAmbassador* and a *FederateAmbassador* (Figure 2) to access the services offered by the RTI. Operations on the *RTIAmbassador* are called by the federate whenever it needs an RTI service (e.g. a request to advance simulation time). In the reverse direction, the RTI invokes an operation on the *FederateAmbassador* whenever it needs to pass data to the federate (e.g. to inform the federate that the request to advance simulation time has been granted). Six classes of services (Figure 1) have to be provided by the RTI to be defined HLA-compliant. These classes are specified within the FIS and they can be summarized as follows:


 Data Distribution Management. These services can reduce unnecessary information transfer between federates by filtering out irrelevant data.

HLA overcame the shortcomings of the DIS standard by being simulation-domain neutral (it was not developed referred to any specific language, therefore HLA provides means to describe any data exchange format as required and specifying functionalities for time management and bandwidth control (see the FIS module).

HLA provides Application Programming Interfaces (APIs) for all the classes of services just mentioned, but the RTI software and algorithms are not defined by HLA. Also the operations in the *FederateAmbassador* need to be implemented at the federate level, as part of the federate code or some interface service (*adapter*).

These facts have caused the growth of multiple HLA-RTI implementations (e.g. [80], [81]) and the development of *ad-hoc* solutions for the adapters on the federate side [25]. In particular the last aspect represents one of the most relevant criticalities in applying HLA for distributed simulation: the lack of a standardized approach to adapt a simulator within an HLA-based distributed architecture, makes a distributed simulation project time expensive since a lot of implementation is required in addition to the effort to build the simulation model.

This consideration represents one of the leading arguments for the research community in the direction of the development of additional complementary standards (Section 3) to ease the creation and management of an HLA-based distributed simulation.

It is the objective of the next section to analyze the state of the art on the adoption and advancements in the use of HLA-based distributed simulation technique.

### **2.2. Distributed simulation in civilian applications**

8 Discrete Event Simulations – Development and Applications

supposed to interact with the RTI.

develop and execute their federations.

federation:

action takes place.

can be summarized as follows:

the federates.

and produce and receive data.

object data during the federation execution.

The HLA standard provides four main components for the realization and management of a

HLA rules (IEEE 1516.0, 2000) representing a set of 10 rules that the simulators

Federate Interface Specification (FIS) (IEEE 1516.2, 2000) defining how simulators are

 Object Model Template (OMT) (IEEE 1516.1, 2000) specifying what kind of information is communicated between simulators and how simulations are documented. Following the OMT each federate defines the data that it is willing to share (publish) with other federates and the data it requires from other federates (subscribe). The resulting object models related to each federate are called simulation object models (SOMs). The federation object model (FOM) combines the federate SOMs into a single object model for

the federation to define the overall data to be exchanged within the federation. Federate Development Process (FEDEP) (IEEE 1516.3, 2004) defining the recommended practice processes and procedures that should be followed by users of the HLA to

The federates cannot directly exchange information throughout the federation, instead the RTI plays the role of the operating system of the distributed simulation, providing a set of general-purpose services for federation management and enabling the federates in carrying out federate-to-federate *interactions*. In particular interactions represent an explicit action taken by a federate that may have some effect on another federate within a federation execution, such action can be tied with a specific time defined as *interactionTime*, when the

Each federate is endowed with an *RTIAmbassador* and a *FederateAmbassador* (Figure 2) to access the services offered by the RTI. Operations on the *RTIAmbassador* are called by the federate whenever it needs an RTI service (e.g. a request to advance simulation time). In the reverse direction, the RTI invokes an operation on the *FederateAmbassador* whenever it needs to pass data to the federate (e.g. to inform the federate that the request to advance simulation time has been granted). Six classes of services (Figure 1) have to be provided by the RTI to be defined HLA-compliant. These classes are specified within the FIS and they

Federation Management. These services allow federates to create and destroy

Declaration Management. These services allow federates to publish federate data and

Object Management. These services allow federate to create and delete object instances,

Ownership Management. These services allow federates to transfer the ownership of

Time Management. These services coordinate the advancement of simulation time of

federation execution and join or resign from an existing federation.

subscribe to updated data produced by other federates.

(federates) have to follow in order to be defined HLA-compliant.

Herein the attention is focused on distributed modeling of complex systems in civilian domain.

HLA constitutes an enabler for implementing the distributed simulation. The standard, though, was conceived for military applications and several problems arise when trying to interoperate heterogeneous simulators in civilian applications (the terminology Commercial off-the-shelf discrete-event simulation packages CSPs [62] will be used to describe commercially available simulators for the analysis of Discrete Event Systems).

Boer [12] investigated the main benefits and criticalities related to the industrial application of HLA by interviewing the actors involved in the problem (e.g. simulation model developers, software houses, HLA experts, [9]-[11]). The results of the survey showed that CSPs vendors do not see direct benefits in using distributed simulation, whereas in industry HLA is considered troublesome because of the lack of experienced users and the complexity of the standard. In addition, as suggested in [49], although the approaches and general methods used in military and civilian simulation communities have similarities, the

terminology turns out to be completely different [36]. For instance, terms like live simulation and virtual emulator are rarely used in civilian applications although equivalent techniques are commonly applied.

Distributed Modeling of Discrete Event Systems 11

Systems, Discrete Event Simulation, Manufacturing Applications, Industrial Application and Civilian Applications. These keywords brought to the identification of 26 core papers based on the number of citations ([4], [12], [8], [11], [9], [10], [20], [28], [29], [30], [33], [74], [75] , [48], [50], [47], [49], [58], [53], [59], [56], [68], [70], [68], [73] and [71]). These papers can be considered as introductory to the topic of distributed simulation in civilian domain. Starting from these articles the bibliographic search followed the path of the citations, i.e. works cited by the core papers and papers citing the core ones were considered. This search brought to the selection of 83 further articles. The overall 109 papers were published mainly in the following journals and conference proceedings: Advanced Simulation Technologies Conference, European Simulation Interoperability Workshop, European Simulation Symposium, Information Sciences, International Journal of Production Research, Journal of the Operational Society, Journal of Simulation, Workshop on Principles of Advanced and

More than 60% (Figure 3) of the analyzed papers propose an application in a specific field of the civilian domain (e.g. [72], [42]). As stated in [46], transportation and logistics are typical application areas of simulation and also the first areas where HLA has been tested by the civilian simulation community. Manufacturing and health care are acquiring increasing importance because of the growth of the extended enterprise and the increase in attention

The main fields of application of DS (Figure 4) are Supply Chain Management (33% of the papers stating the domain of interest) (e.g. [64], [22], [42]), Manufacturing (29% of the papers) (e.g. [69], [77]), Health Care (e.g. [34]) and Production Scheduling & Maintenance

A further analysis was carried out by considering only the articles related to the manufacturing domain, aiming at evaluating whether the contributions addressed real industrial case applications or test cases applications. Only 22% of the articles address a real case, thus confirming the outcomes obtained by Boer [8] in the analysis of the adoption of

Distributed Simulation and Winter Simulation Conference.

**Figure 3.** Overall Classification Criteria

for bio-pharmaceutical supply chains respectively.

(e.g. [72]), 17% of the articles are related to Health Care.

*2.2.2. Domain of application* 

The major difference between military and civilian domain resides in the way simulation models are developed and what are the goals to meet when starting a simulation development process. In the military community where time and budget constraints are not the key elements leading the building process of a simulation tool, languages such as C++ and Java are usually adopted because of their flexibility. On the other hand, in the civilian simulation community, the use of commercial simulation tools (e.g. Arena, Automod, Simio, ProModel, Simple++, SLX, etc.) is the common practice. These tools satisfy the need of rapidly and cost-effectively developing the simulation models.

The use of commercial simulation tools hinders the applicability of the HLA standard for the realization of a distributed simulation model, because the direct access to the HLA APIs (Section 2.1.) from the commercial simulation software tools is not usually possible. Therefore, the enhancement of HLA with additional complementary standards [51] and the definition of a standard language for CSPs represent relevant and not yet solved technical and scientific challenges ([25], [49], [50]). Recently, the COTS Simulation Package Interoperability-Product Development Group (CSPI-PDG), within the Simulation Interoperability Standards Organization (SISO), worked on the definition of the CSP interoperability problem (Interoperability Reference Models, IRMs) [74] and on a draft proposal for a standard to support the CSPs interoperability (Entity Transfer Specification, ETS) [61].

#### *2.2.1. Literature review*

The application of distributed simulation in the civilian domain has been studied by reviewing the available literature with the purpose to individuate which civilian domain distributed simulation is generally called, which motivations underlie the adoption of the distributed technique, which technical and scientific challenges have been faced and which solutions have been proposed so far. More than 100 papers have been analyzed and classified according to three criteria:


Most of the articles can be classified according to more than one criterion and Figure 3 shows the percentage of articles falling in each category.

The bibliographic search was carried out by considering the following keywords: Distributed Simulation, Operations Research and Management, Commercial Simulation Packages, Interoperability Reference Models, High Level Architecture, Manufacturing Systems, Discrete Event Simulation, Manufacturing Applications, Industrial Application and Civilian Applications. These keywords brought to the identification of 26 core papers based on the number of citations ([4], [12], [8], [11], [9], [10], [20], [28], [29], [30], [33], [74], [75] , [48], [50], [47], [49], [58], [53], [59], [56], [68], [70], [68], [73] and [71]). These papers can be considered as introductory to the topic of distributed simulation in civilian domain. Starting from these articles the bibliographic search followed the path of the citations, i.e. works cited by the core papers and papers citing the core ones were considered. This search brought to the selection of 83 further articles. The overall 109 papers were published mainly in the following journals and conference proceedings: Advanced Simulation Technologies Conference, European Simulation Interoperability Workshop, European Simulation Symposium, Information Sciences, International Journal of Production Research, Journal of the Operational Society, Journal of Simulation, Workshop on Principles of Advanced and Distributed Simulation and Winter Simulation Conference.

**Figure 3.** Overall Classification Criteria

#### *2.2.2. Domain of application*

10 Discrete Event Simulations – Development and Applications

rapidly and cost-effectively developing the simulation models.

support the CSPs interoperability (Entity Transfer Specification, ETS) [61].

are commonly applied.

*2.2.1. Literature review* 

classified according to three criteria:

terminology turns out to be completely different [36]. For instance, terms like live simulation and virtual emulator are rarely used in civilian applications although equivalent techniques

The major difference between military and civilian domain resides in the way simulation models are developed and what are the goals to meet when starting a simulation development process. In the military community where time and budget constraints are not the key elements leading the building process of a simulation tool, languages such as C++ and Java are usually adopted because of their flexibility. On the other hand, in the civilian simulation community, the use of commercial simulation tools (e.g. Arena, Automod, Simio, ProModel, Simple++, SLX, etc.) is the common practice. These tools satisfy the need of

The use of commercial simulation tools hinders the applicability of the HLA standard for the realization of a distributed simulation model, because the direct access to the HLA APIs (Section 2.1.) from the commercial simulation software tools is not usually possible. Therefore, the enhancement of HLA with additional complementary standards [51] and the definition of a standard language for CSPs represent relevant and not yet solved technical and scientific challenges ([25], [49], [50]). Recently, the COTS Simulation Package Interoperability-Product Development Group (CSPI-PDG), within the Simulation Interoperability Standards Organization (SISO), worked on the definition of the CSP interoperability problem (Interoperability Reference Models, IRMs) [74] and on a draft proposal for a standard to

The application of distributed simulation in the civilian domain has been studied by reviewing the available literature with the purpose to individuate which civilian domain distributed simulation is generally called, which motivations underlie the adoption of the distributed technique, which technical and scientific challenges have been faced and which solutions have been proposed so far. More than 100 papers have been analyzed and

*Domain of application*, i.e. the specific civilian sector where the distributed simulation

*Motivation* underlying the adoption of the distributed simulation, i.e. the main problem

*Technical issue* faced, i.e. the solutions to integration issue or enhancement to services of

Most of the articles can be classified according to more than one criterion and Figure 3

The bibliographic search was carried out by considering the following keywords: Distributed Simulation, Operations Research and Management, Commercial Simulation Packages, Interoperability Reference Models, High Level Architecture, Manufacturing

was applied (e.g. manufacturing domain, health care, emergency, etc.).

leading to the adoption of the distributed simulation architecture.

the HLA architecture proposed within the considered article.

shows the percentage of articles falling in each category.

More than 60% (Figure 3) of the analyzed papers propose an application in a specific field of the civilian domain (e.g. [72], [42]). As stated in [46], transportation and logistics are typical application areas of simulation and also the first areas where HLA has been tested by the civilian simulation community. Manufacturing and health care are acquiring increasing importance because of the growth of the extended enterprise and the increase in attention for bio-pharmaceutical supply chains respectively.

The main fields of application of DS (Figure 4) are Supply Chain Management (33% of the papers stating the domain of interest) (e.g. [64], [22], [42]), Manufacturing (29% of the papers) (e.g. [69], [77]), Health Care (e.g. [34]) and Production Scheduling & Maintenance (e.g. [72]), 17% of the articles are related to Health Care.

A further analysis was carried out by considering only the articles related to the manufacturing domain, aiming at evaluating whether the contributions addressed real industrial case applications or test cases applications. Only 22% of the articles address a real case, thus confirming the outcomes obtained by Boer [8] in the analysis of the adoption of

distributed simulation in the manufacturing environment. Although solutions have been developed for the manufacturing domain, this technique is still far from being adopted as an evaluation tool by industrial companies because the end-users perceive HLA and distributed simulation as an additional trouble rather than a promising approach [10]. As a consequence, a lot of effort is put in the development of decision support systems that hide the complexity of a distributed environment to the end user [41].

Distributed Modeling of Discrete Event Systems 13

Over 70% of the articles deal with technical issues, thus showing that HLA and DS experts are putting a lot of effort in the enhancement and extensions of HLA-based DS to face civilian application problems. Indeed, the application of distributed simulation to civilian domain still presents several technical issues. In particular four main research areas can be

1. *Integration of commercial discrete event simulators (CSP)*. Several CSPs are put together and

2. *Interoperability reference models and entity transfer*. The papers in this category work in the standardization of the communication between federates within an HLA-compliant

3. *Time management enhancement*. The issues related to the time synchronization of

4. *RTI-services extension*. In this case the services listed in Section 2.1. are enhanced for

The outcome of the review was that the integration of CSPs is the most addressed technical issue, (45% of the papers) nonetheless the integration of real CSPs (i.e. not general purpose programming languages) still represents a challenging topic. Figure 5 gives a picture of the main CSP solutions that have been adopted in the literature. In particular, the *y-axis* reports the percentage of articles that use one of the listed CSPs (e.g. [37], [21], [67]) within the papers that deal with the interoperation of simulators. It can be noticed that CSP Emulators (e.g. [68], [38]) are still one of the most used solutions because of the problems related to interoperating real CSPs. These problems are mainly caused by the lack of data and information mapping between simulators and the difficulty in interacting (e.g. send and

receive information, share the internal event list) while the simulation is running.

The enhancement of the Run Time Infrastructure services is another key research topic ([2], [19]). In particular, the scientific articles deal with two open issues: (1) Time management (e.g. [39], [31]), (2) Data Distribution Management (e.g. [66], [73]). Time Management has

synchronized by means of the services offered by the HLA infrastructure.

*2.2.4. Technical issue faced* 

federation (Section 3.).

specific applications [82].

federates are faced.

**Figure 5.** CSP adopted

identified:

**Figure 4.** Distributed Simulation Main fields of application

#### *2.2.3. Motivations underlying the adoption of the distributed simulation*

As Boer stated in [8], if a problem can be solved by a monolithic simulation model created in a single COTS simulation package and a distributed approach is not explicitly required the simulation practitioner should certainly choose the monolithic solution in the selected CSP. Similarly, Strassburger [49] suggests that if a maintainable and reusable monolithic application can be built, then there is no point in building it in a distributed platform.

However, there are simulation projects where the distributed solution seems more advantageous and straightforward [38] because it enables to cope with:


All the papers stating a motivation for using DS mention the system complexity (e.g. [22], [72], [30], [32]), whereas 44% of the papers the demand for reuse [78]. The low percentage (around 5%) of papers using DS to cope with lack of shared information can be partially traced back to the lack of real industrial applications that still characterizes DS in civilian environment [76].

#### *2.2.4. Technical issue faced*

12 Discrete Event Simulations – Development and Applications

the complexity of a distributed environment to the end user [41].

**Figure 4.** Distributed Simulation Main fields of application

related to their subsystems.

whole simulation model.

environment [76].

*2.2.3. Motivations underlying the adoption of the distributed simulation* 

advantageous and straightforward [38] because it enables to cope with:

As Boer stated in [8], if a problem can be solved by a monolithic simulation model created in a single COTS simulation package and a distributed approach is not explicitly required the simulation practitioner should certainly choose the monolithic solution in the selected CSP. Similarly, Strassburger [49] suggests that if a maintainable and reusable monolithic

However, there are simulation projects where the distributed solution seems more

1. Demand for *reusability* of the simulator output of the simulation project. Here the word reusability is adopted both in terms of the possibility to reuse simulators already

3. *System complexity*. In this case no one modeler has enough knowledge to realize the

All the papers stating a motivation for using DS mention the system complexity (e.g. [22], [72], [30], [32]), whereas 44% of the papers the demand for reuse [78]. The low percentage (around 5%) of papers using DS to cope with lack of shared information can be partially traced back to the lack of real industrial applications that still characterizes DS in civilian

developed and of building new simulators that can be readopted in the future. 2. *Lack of Shared information*. This is the case when no one modeler has enough information to develop the simulator. This condition holds when the whole system to be modeled is divided into subsystems owned by different actors that do not want to share data

application can be built, then there is no point in building it in a distributed platform.

distributed simulation in the manufacturing environment. Although solutions have been developed for the manufacturing domain, this technique is still far from being adopted as an evaluation tool by industrial companies because the end-users perceive HLA and distributed simulation as an additional trouble rather than a promising approach [10]. As a consequence, a lot of effort is put in the development of decision support systems that hide

Over 70% of the articles deal with technical issues, thus showing that HLA and DS experts are putting a lot of effort in the enhancement and extensions of HLA-based DS to face civilian application problems. Indeed, the application of distributed simulation to civilian domain still presents several technical issues. In particular four main research areas can be identified:


The outcome of the review was that the integration of CSPs is the most addressed technical issue, (45% of the papers) nonetheless the integration of real CSPs (i.e. not general purpose programming languages) still represents a challenging topic. Figure 5 gives a picture of the main CSP solutions that have been adopted in the literature. In particular, the *y-axis* reports the percentage of articles that use one of the listed CSPs (e.g. [37], [21], [67]) within the papers that deal with the interoperation of simulators. It can be noticed that CSP Emulators (e.g. [68], [38]) are still one of the most used solutions because of the problems related to interoperating real CSPs. These problems are mainly caused by the lack of data and information mapping between simulators and the difficulty in interacting (e.g. send and receive information, share the internal event list) while the simulation is running.

The enhancement of the Run Time Infrastructure services is another key research topic ([2], [19]). In particular, the scientific articles deal with two open issues: (1) Time management (e.g. [39], [31]), (2) Data Distribution Management (e.g. [66], [73]). Time Management has

received more attention (91% of the papers dealing with enhancement of RTI services) because it strongly influences the computational performance of the distributed simulation.

Distributed Modeling of Discrete Event Systems 15

An initial set of interoperability problems identified by the CSPI-PDG have been divided into a series of problem types that are represented by IRMs. The purpose of an IRM can be to [54]:

Clearly identify the CSP/model interoperability capabilities of an existing distributed

Clearly specify the CSP/model interoperability requirements of a proposed distributed

The literature review showed that around 21% of the articles dealing with technical issues (Section 2.2.1.) taken into consideration deal with IRMs (e.g., [39], [56], [51], [63], [55] and

IRM Type A Entity Transfer represents interoperability problems that can occur when transferring an entity from one simulation model to another. This IRM type is the most formalized at the present moment, since the need to transfer entities between simulators has

Figure 6 shows an illustrative example of the problem of Entity Transfer where an entity *e1* leaves activity A1 in model M1 at time T1 and arrives at queue Q2 in model M2 at time T2. For example, if M1 is a car production line and M2 is a paint shop, then the entity transfer happens when a car leaves M1 at T1 and arrives in a buffer in M2 at T2 to wait for painting.

 IRM Type A.1 General Entity Transfer is defined as the transfer of entities from one model to another such that an entity *e1* leaves model M1 at T1 from a given place and arrives at model M2 at T2 at a given place and T1 ≤ T2 or T1 < T2. The place of

 IRM Type A.2 Bounded Receiving Element. The IRM Type A.2 is defined as the relationship between an activity A in a model M1 and a bounded queue Q2 in a model M2 such that if an entity e is ready to leave activity A at T1 and attempts to arrive at the

If the bounded queue Q2 is empty, the entity e can leave activity A at T1 and arrive

If the bounded queue Q2 is full, the entity e cannot leave activity A at T1; activity A

 When the bounded queue Q2 becomes not full at T3, entity e must leave A at T3 and arrive at Q2 at T4, activity A becomes unblocked and may receive new entities at T3.

may then block if appropriate and must not accept any more entities.

been the most popular feature requested from the distributed simulation users so far.

simulation.

simulation.

[34]).

There are four types of IRM:

*3.1.1. IRm type A: Entity transfer* 

There are three subtypes of IRM Type A:

bounded queue Q2 at T2 then:

at Q2 at T2

departure and arrival will be a queue, workstation, etc.

 Type A - Entity Transfer (Section 3.1.1.) Type B - Shared Resource (Section 3.1.2.) Type C - Shared Events (Section 3.1.3.)

Type D - Shared Data Structure (Section 3.1.4)

The following conclusions can be drawn from the literature analysis:


The issues faced to model complex systems give raise to problems in the distributed simulation realization that are strongly dependent on the specific application case [57] and the solutions to those needs can be implemented through an RTI in many different (incompatible) ways. Each way can be promising in its own context, but the lack of a standardized approach means it is difficult for end users and CSP vendors to choose a solution thus slowing down the spreading of the distributed simulation technique.

## **3. A standard based approach for distributed simulation**

The main contribution in standardization of distributed modeling has to be credited to Simulation Interoperability Standard Organization (SISO) and in particular to the High Level Architecture Simulation Package Interoperability Forum (HLA-CSPIF). HLA-CSPIF and, then, COTS Simulation Package Interoperability Product Development Group (CSPI-PDG) were created in an attempt to produce a generalizable solution to the problem of integrating distributed heterogeneous CSPs.

As highlighted at the end of Section 2.2., a standardized approach is fundamental to increase the use of distributed simulation in civilian applications. This led to formalize the problem of the interaction between simulators in civilian applications and to standardize the way data are exchanged between federates within the federation.

The main results of the standardization effort are the Interoperability Reference Models (IRMs) and the Entity Transfer Specification (ETS), that will be presented in Section 3.1 and 3.2, respectively. Section 3.3 shows the distributed simulation communication protocol presented in [77], based on IRMs and the extended ETS proposed in [38].

#### **3.1. Interoperability reference model**

An interoperability problem type is meant to capture a general class of interoperability problems, whereas an IRM is meant to capture a specific problem within that class.

The creation of the IRMs has proved to be a powerful tool in the development of standards in the distributed simulation research community, as it is now possible to create solutions for specific integration problems.

An initial set of interoperability problems identified by the CSPI-PDG have been divided into a series of problem types that are represented by IRMs. The purpose of an IRM can be to [54]:


There are four types of IRM:

14 Discrete Event Simulations – Development and Applications

adapted to civilian applications.

integrating distributed heterogeneous CSPs.

**3.1. Interoperability reference model** 

for specific integration problems.

data are exchanged between federates within the federation.

presented in [77], based on IRMs and the extended ETS proposed in [38].

environments.

received more attention (91% of the papers dealing with enhancement of RTI services) because it strongly influences the computational performance of the distributed simulation.

There is a lack of distributed simulation applications in real manufacturing

The HLA architecture components (in particular RTI services) must be extended and

The issues faced to model complex systems give raise to problems in the distributed simulation realization that are strongly dependent on the specific application case [57] and the solutions to those needs can be implemented through an RTI in many different (incompatible) ways. Each way can be promising in its own context, but the lack of a standardized approach means it is difficult for end users and CSP vendors to choose a

The main contribution in standardization of distributed modeling has to be credited to Simulation Interoperability Standard Organization (SISO) and in particular to the High Level Architecture Simulation Package Interoperability Forum (HLA-CSPIF). HLA-CSPIF and, then, COTS Simulation Package Interoperability Product Development Group (CSPI-PDG) were created in an attempt to produce a generalizable solution to the problem of

As highlighted at the end of Section 2.2., a standardized approach is fundamental to increase the use of distributed simulation in civilian applications. This led to formalize the problem of the interaction between simulators in civilian applications and to standardize the way

The main results of the standardization effort are the Interoperability Reference Models (IRMs) and the Entity Transfer Specification (ETS), that will be presented in Section 3.1 and 3.2, respectively. Section 3.3 shows the distributed simulation communication protocol

An interoperability problem type is meant to capture a general class of interoperability

The creation of the IRMs has proved to be a powerful tool in the development of standards in the distributed simulation research community, as it is now possible to create solutions

problems, whereas an IRM is meant to capture a specific problem within that class.

The interoperability of CSPs still represents a technical challenging problem.

solution thus slowing down the spreading of the distributed simulation technique.

**3. A standard based approach for distributed simulation** 

The following conclusions can be drawn from the literature analysis:


The literature review showed that around 21% of the articles dealing with technical issues (Section 2.2.1.) taken into consideration deal with IRMs (e.g., [39], [56], [51], [63], [55] and [34]).

#### *3.1.1. IRm type A: Entity transfer*

IRM Type A Entity Transfer represents interoperability problems that can occur when transferring an entity from one simulation model to another. This IRM type is the most formalized at the present moment, since the need to transfer entities between simulators has been the most popular feature requested from the distributed simulation users so far.

Figure 6 shows an illustrative example of the problem of Entity Transfer where an entity *e1* leaves activity A1 in model M1 at time T1 and arrives at queue Q2 in model M2 at time T2. For example, if M1 is a car production line and M2 is a paint shop, then the entity transfer happens when a car leaves M1 at T1 and arrives in a buffer in M2 at T2 to wait for painting.

There are three subtypes of IRM Type A:

	- If the bounded queue Q2 is empty, the entity e can leave activity A at T1 and arrive at Q2 at T2
	- If the bounded queue Q2 is full, the entity e cannot leave activity A at T1; activity A may then block if appropriate and must not accept any more entities.
	- When the bounded queue Q2 becomes not full at T3, entity e must leave A at T3 and arrive at Q2 at T4, activity A becomes unblocked and may receive new entities at T3.


Distributed Modeling of Discrete Event Systems 17

*3.1.2. IRM type B: Shared resource* 

status of the resource..

same at T1.

time

*3.1.3. Type C: Shared event* 

*3.1.4. Type D: Shared data structures* 

**3.2. Entity transfer specification** 

refers to the architecture shown in Figure 8.

such that:

IRM Type B deals with the problem of sharing resources across two or more models in a distributed simulation. A modeler can specify if an activity requires a resource (such as machine operators, conveyor, pallets, etc.) of a particular type to begin. If an activity does require a resource, when an entity is ready to start that activity, it must therefore be determined if there is a resource available. If it is available then the resource is secured by the activity and held until the activity ends. A resource shared by two or more models leads to a problem of maintaining the consistency change the deleted part with: related to the

Currently there is only one IRM Type B subtype. The IRM Type B.1 General Shared Resources is defined as the maintenance of consistency of all copies of a shared resource

 if a model M1 wishes to change its copy of the resource at time T1 then the value of all other copies of the same resource present in other models will be guaranteed to be the

if two or more models wish to change their copies of the resource at the same time T1,

IRM Type C deals with the problem of sharing events (such as an emergency signal,

There is currently one IRM Type C sub-type. The IRM Type C.1 General Shared Event is

 if a model M1 wishes to schedule a shared event E at T1, then its local copies (i.e. the events scheduled within each simulator) will be guaranteed to be executed at the same

 if two or more models wish to schedule shared events E1, E2, etc. at T1, then all local copies of all shared events will be guaranteed to be executed at the same time T1.

IRM Type D Shared Data Structure deals with the sharing of variables and data structures across simulation models that are semantically different to resources (e.g. while referring to the supply chain environment, this is the case with bill of material information management

CSPI-PDG proposed the Entity Transfer Specification (ETS) Protocol ([51], [52], [62]) which

Each federate consists of a COTS simulation package (CSP), a model that is executed by the CSP, and the middleware that is a sort of adaptor interfacing the CSP with the Run Time

defined as the guaranteed execution of all local copies of a shared event E such that:

then all copies will be guaranteed to be the same at T1.

explosion, etc.) across two or more models in a distributed simulation.

or shared inventory). There is currently one IRM Type D sub-type.

Note that the IRM sub-types are intended to be cumulative, i.e. a distributed simulation that correctly transfers entities from one model to a bounded buffer in another model should be compliant with both IRM Type A.1 General Entity Transfer and IRM Type A.2 Bounded Receiving Element.

The largest part of the papers analyzed in the literature review deal with the basic IRM that is the general entity transfer (85% of the papers dealing with IRMs), since most of the applications are related to Supply Chain Management and queues can often be modeled as infinite capacity as they represent inventory, production or distribution centers.

The situation is slightly different if the manufacturing domain is considered. Indeed, IRM Type A.2 (70% of the articles dealing with IRMs), is largely adopted when a production system is modeled, because the decoupling buffers between workstations must be usually represented as finite capacity queues.

**Figure 6.** Typr A.1 IRM

**Figure 7.** Type A.3 IRM

#### *3.1.2. IRM type B: Shared resource*

16 Discrete Event Simulations – Development and Applications

priority ordering is dynamic or specialized.

If activity A is blocked then the simulation of model M1 must continue.

 IRM Type A.3 Multiple Input Prioritization. As shown in Figure 7, the IRM Type A.3 Multiple Input Prioritization represents the case where a model element such as queue Q1 (or workstation) can receive entities from multiple places. Let us assume that there are two models M2 and M3 which are capable of sending entities to Q1 and that Q1 has a First-In-First-Out (FIFO) queuing discipline. If an entity *e1* is sent from M2 at T1 and arrives at Q1 at T2 and an entity *e2* is sent from M3 at T3 and arrives at Q1 at T4, then if T2 < T4 we would expect the order of entities in Q1 would be e1, e2. A problem arises when both entities arrive at the same time, i.e. when T2 = T4. Depending on implementation, the order of entities would either be *e1*, *e2* or *e2*, *e1*. In some modeling situations it is possible to specify the priority order if such a conflict arises, e.g. it can be specified that model M1 entities will always have a higher priority than model M2 (and therefore require the entity order *e1*, *e2* if T2 = T4). Furthermore, it is possible that this

Note that the IRM sub-types are intended to be cumulative, i.e. a distributed simulation that correctly transfers entities from one model to a bounded buffer in another model should be compliant with both IRM Type A.1 General Entity Transfer and IRM Type A.2 Bounded

The largest part of the papers analyzed in the literature review deal with the basic IRM that is the general entity transfer (85% of the papers dealing with IRMs), since most of the applications are related to Supply Chain Management and queues can often be modeled as

The situation is slightly different if the manufacturing domain is considered. Indeed, IRM Type A.2 (70% of the articles dealing with IRMs), is largely adopted when a production system is modeled, because the decoupling buffers between workstations must be usually

infinite capacity as they represent inventory, production or distribution centers.

T1 ≤ T2 and T3 ≤ T4.

Receiving Element.

**Figure 6.** Typr A.1 IRM

**Figure 7.** Type A.3 IRM

represented as finite capacity queues.

IRM Type B deals with the problem of sharing resources across two or more models in a distributed simulation. A modeler can specify if an activity requires a resource (such as machine operators, conveyor, pallets, etc.) of a particular type to begin. If an activity does require a resource, when an entity is ready to start that activity, it must therefore be determined if there is a resource available. If it is available then the resource is secured by the activity and held until the activity ends. A resource shared by two or more models leads to a problem of maintaining the consistency change the deleted part with: related to the status of the resource..

Currently there is only one IRM Type B subtype. The IRM Type B.1 General Shared Resources is defined as the maintenance of consistency of all copies of a shared resource such that:


#### *3.1.3. Type C: Shared event*

IRM Type C deals with the problem of sharing events (such as an emergency signal, explosion, etc.) across two or more models in a distributed simulation.

There is currently one IRM Type C sub-type. The IRM Type C.1 General Shared Event is defined as the guaranteed execution of all local copies of a shared event E such that:


#### *3.1.4. Type D: Shared data structures*

IRM Type D Shared Data Structure deals with the sharing of variables and data structures across simulation models that are semantically different to resources (e.g. while referring to the supply chain environment, this is the case with bill of material information management or shared inventory). There is currently one IRM Type D sub-type.

#### **3.2. Entity transfer specification**

CSPI-PDG proposed the Entity Transfer Specification (ETS) Protocol ([51], [52], [62]) which refers to the architecture shown in Figure 8.

Each federate consists of a COTS simulation package (CSP), a model that is executed by the CSP, and the middleware that is a sort of adaptor interfacing the CSP with the Run Time

Infrastructure (RTI) (Figure 8). The relationship between CSP, the middleware and the RTI consists of two communication flows: (1) middleware-RTI, (2) CSP-middleware. The middleware translates the simulation information into a common format so that the RTI can share it with the federation. In addition, the middleware receives and sends information from/to the CSP. The CSP communicates with its middleware by means of Simulation Messages (Section 3.3.1.) [77]. The presence of Simulation Messages is the main difference between the reference architecture in Figure 8 and the architecture proposed by Taylor et al [51].

Distributed Modeling of Discrete Event Systems 19

 The representation of an entity depends on how the simulation model is designed and implemented in a CSP. Indeed, the names that the modelers use to represent the same entity might be different. A similar problem can arise for the definition of simple datatypes. For example, some CSPs use 32-bit real numbers while others use 64-bit [62].

 ETS suggested interaction hierarchy does not work: a federate subscribing to the superclass will never receive the values transmitted in the interaction parameters. The specification of user defined attributes is placed into a complex datatype, this introduces new room for interoperability challenges as all participating federates have

 There are some possibilities for misinterpretation in the definition of Entity and EntityType introducing changes in FOMs whenever a new entity type is talked about.

Furthermore, the ETS was not designed to manage the Type A.2. IRM and the interaction class hierarchy refers to the entity transfer without taking into account any information on

One of the most recent contributions in ETS was presented by Pedrielli et al. [77] and consists in the proposal of a new class hierarchy. In particular, different subclasses of the transferEntity class were defined to enable the differentiation of multiple connections between models and the Type A.2. IRM management. After developing the interaction class hierarchy, following the HLA standard, the Simulation Object Model (SOM) and Federation Object Model (FOM) were developed to include the novel interactions and their parameters. In particular, extensions were proposed to the Interaction Class Table (part of the OMT, Section 2.1) to include the novel interaction classes and define them as publish and/or subscribe. The Parameter Table (part of the OMT, Section 2.1) was modified to include the

proposed parameters for the interactions and the Datatype table was also modified.

 *transferEntity*, as already defined in the ETS protocol. This superclass allows the federate subscribing to all the instances of entity transfer. The instantiation of this class

 *transferEntityFromFedSourceEx* is a novel subclass defined for every exit point, where FedSourceEx stands for the name or abbreviation of a specific exit point in the sending model. This class is useful to group the instances of the transferEntity that are related to the source federate, so that the FedSourceEx can subscribe to all these instances without

 *transferEntityFromFedSourceExToFedDestEn* is a novel subclass defined for each pair of exit point (Ex) of the source federate (FedSource) and entry point (En) of the receiving federate (FedDest). This class is instantiated both when a sending model needs to transfer a part to a specific entry point in the receiving model, and when a receiving model needs to share information about a buffer state or about the receipt of a part from

The resulting class hierarchy consists of the following classes [38]:

is related to visualization and monitoring tasks.

explicitly naming them.

Straburger [50] highlighted some relevant drawbacks in the ETS standard proposal: It is not possible to differentiate multiple connections between any two models.

to be able to interpret all of the attributes.

the state of the receiving buffer (e.g. Q2 in Figure 6).

ETS defines the communication between the sending model and the receiving model (*ModelA* and *ModelB* in Figure 8, respectively) at RTI level. In particular the way the middleware of each federate and the RTI exchange information is formalized by means of a special hierarchy of interaction classes. An interaction class is defined as a template for a set of characteristics (parameters) that are shared by a group of interactions (refer to IEEE HLA standard, 2000). The middleware of the sending model instantiates a specific interaction class and sends it to the RTI whenever an entity has to be transferred.

**Figure 8.** ETS refrence Architecture

Two main issues arise when the simulation information is translated for the RTI:

 A common time definition and resolution is necessary. For example, the time should be defined as being the time when an entity exits a source model and then instantaneously arrives at the destination model (i.e. the definition of time implies zero transit time) [62]. Alternatively, it should be defined including the notion of travel time and the entity would arrive at destination with a delay equal to the transfer time.

 The representation of an entity depends on how the simulation model is designed and implemented in a CSP. Indeed, the names that the modelers use to represent the same entity might be different. A similar problem can arise for the definition of simple datatypes. For example, some CSPs use 32-bit real numbers while others use 64-bit [62].

Straburger [50] highlighted some relevant drawbacks in the ETS standard proposal:

18 Discrete Event Simulations – Development and Applications

**Figure 8.** ETS refrence Architecture

[51].

Infrastructure (RTI) (Figure 8). The relationship between CSP, the middleware and the RTI consists of two communication flows: (1) middleware-RTI, (2) CSP-middleware. The middleware translates the simulation information into a common format so that the RTI can share it with the federation. In addition, the middleware receives and sends information from/to the CSP. The CSP communicates with its middleware by means of Simulation Messages (Section 3.3.1.) [77]. The presence of Simulation Messages is the main difference between the reference architecture in Figure 8 and the architecture proposed by Taylor et al

ETS defines the communication between the sending model and the receiving model (*ModelA* and *ModelB* in Figure 8, respectively) at RTI level. In particular the way the middleware of each federate and the RTI exchange information is formalized by means of a special hierarchy of interaction classes. An interaction class is defined as a template for a set of characteristics (parameters) that are shared by a group of interactions (refer to IEEE HLA standard, 2000). The middleware of the sending model instantiates a specific interaction

class and sends it to the RTI whenever an entity has to be transferred.

Two main issues arise when the simulation information is translated for the RTI:

would arrive at destination with a delay equal to the transfer time.

 A common time definition and resolution is necessary. For example, the time should be defined as being the time when an entity exits a source model and then instantaneously arrives at the destination model (i.e. the definition of time implies zero transit time) [62]. Alternatively, it should be defined including the notion of travel time and the entity


Furthermore, the ETS was not designed to manage the Type A.2. IRM and the interaction class hierarchy refers to the entity transfer without taking into account any information on the state of the receiving buffer (e.g. Q2 in Figure 6).

One of the most recent contributions in ETS was presented by Pedrielli et al. [77] and consists in the proposal of a new class hierarchy. In particular, different subclasses of the transferEntity class were defined to enable the differentiation of multiple connections between models and the Type A.2. IRM management. After developing the interaction class hierarchy, following the HLA standard, the Simulation Object Model (SOM) and Federation Object Model (FOM) were developed to include the novel interactions and their parameters. In particular, extensions were proposed to the Interaction Class Table (part of the OMT, Section 2.1) to include the novel interaction classes and define them as publish and/or subscribe. The Parameter Table (part of the OMT, Section 2.1) was modified to include the proposed parameters for the interactions and the Datatype table was also modified.

The resulting class hierarchy consists of the following classes [38]:


a specific exit point in a sending model. The models both publish and subscribe to this subclass that was designed to create a private communication channel between the sending and the receiving model. Therefore, if an entry point in the receiving model is connected with multiple federates/exit points, then the receiving federate has to inform about the state of the entry point by means of multiple interactions, each dedicated to a specific federate/exit point. This communication strategy is not the most efficient in a generic case, but it offers the possibility to deliver customized information and adopt different priorities for the various federates/exit points. This becomes fundamental in real industrial applications where information sharing among different subsystems is seen as a threat, thus rising the need to design a protocol that creates a one to one communication between each pair of exit/entry point inside the corresponding sending/receiving model.

Distributed Modeling of Discrete Event Systems 21

Interaction Parameter DataType Transportation Order

Record Name Field Encoding

**3.3. Communication within the HLA-based integration infrastructure** 

Pedrielli et al. [77] proposed a communication protocol (see Section 3.3.2) based on messages to manage the communication between a CSP and its middleware (or adapter). The communication protocol was conceived for the distributed simulation of network of Discrete Event Manufacturing Systems characterized by the transfer of parts in the presence of buffers with finite capacity, with the objective to minimize the use of zero-lookahead [62]

Before illustrating the communication protocol, Section 3.3.1. presents the concept and functioning of Simulation Messages created to support the communication between a CSP and the middleware. The communication protocol between federates is then explained in Section 3.3.2., whereas Section 3.3.3. defines the hypotheses needed to minimize the zero

The function of the simulation messages depends on the role played by the federate. The sending federate uses the message for communications concerning the need of sending an entity to another model (outgoing communication) and/or information on the availability of the target receiving federate (incoming communication). The receiving federate uses the message for communications concerning the buffer state and/or the acceptance of an entity (outgoing communication) and/or the receipt of an entity from other models (incoming communication). Simulation Messages are implemented as a class that is characterized by

 time referring to the simulation time when the message is sent to the middleware from the CSP. This attribute is used by the middleware to determine the TimeStamp of the

EntityType EntityName HLAASCIIString EntityNumber HLAInteger32BE

Name Type HLAfixedRecord

Entity EntityType HLAreliable TimeStamp ReceivedEntity EntityType HLAreliable TimeStamp Buffer\_Availability HLAInteger32BE HLAreliable TimeStamp SourcePriority HLAInteger32BE HLAreliable TimeStamp EntityTransferTime HLAFloat32BE HLAreliable TimeStamp

TransferEnti tyFromFedS ourceExAtoF edDestEnC (P/S)

**Table 2.** Parameter Table

**Table 3.** Fixed Record Datatype table

for the synchronization of federates.

*3.3.1. Simulation messages* 

the following attributes:

lookahead when applying the proposed protocol.

interaction that will be sent to the RTI.

The ETS Interaction class table was modified to represent the *transferEntityFromFedSourceEx* and *transferEntityFromFedSourceExToFedDestEn* subclasses. The Parameter Table was modified to include the parameters of the novel interaction class *transferEntity FromFedSourceExToFedDestEn*. The introduced parameters are presented below. The similarities with the parameters included in the ETS Parameter Table are highlighted where present.



The resulting tables are shown in Tables 1, 2 and 3.

**Table 1.** Table 1. Interaction Class Table

Distributed Modeling of Discrete Event Systems 21


**Table 2.** Parameter Table

20 Discrete Event Simulations – Development and Applications

present.

type of the parameter Entity.

entity instantaneously arrives at destination.

TransferEntity

The resulting tables are shown in Tables 1, 2 and 3.

(N/S)

**Table 1.** Table 1. Interaction Class Table

availability.

HLAinteractionRoo

t(N)

IRM (Section 3.1)

a specific exit point in a sending model. The models both publish and subscribe to this subclass that was designed to create a private communication channel between the sending and the receiving model. Therefore, if an entry point in the receiving model is connected with multiple federates/exit points, then the receiving federate has to inform about the state of the entry point by means of multiple interactions, each dedicated to a specific federate/exit point. This communication strategy is not the most efficient in a generic case, but it offers the possibility to deliver customized information and adopt different priorities for the various federates/exit points. This becomes fundamental in real industrial applications where information sharing among different subsystems is seen as a threat, thus rising the need to design a protocol that creates a one to one communication between each pair of exit/entry point inside the corresponding sending/receiving model.

The ETS Interaction class table was modified to represent the *transferEntityFromFedSourceEx* and *transferEntityFromFedSourceExToFedDestEn* subclasses. The Parameter Table was modified to include the parameters of the novel interaction class *transferEntity FromFedSourceExToFedDestEn*. The introduced parameters are presented below. The similarities with the parameters included in the ETS Parameter Table are highlighted where

 *Entity*. It is a parameter of complex datatype containing the *EntityName* that is used to communicate the type of the entity, and the *EntityNumber* that is used to communicate the number of entities to be transferred. The *EntityName* and *EntityNumber* play the role

*ReceivedEntity*. It refers to the entity received by the receiving federate and has the same

*Buffer\_Availability*. It was designed to enable the communication about the buffer

 *SourcePriority*. This parameter was designed to communicate the priority assigned to the entity source, so that the infrastructure can be further extended to manage Type A.3

 *EntityTransferTime*. It defines the simulation time when the entity is transferred to the destination point, i.e. the arrival time. Herein the entity leaves the source node and reaches the destination node at the same time, since it is assumed that the transferred

> TransferEntityFrom FromFedSourceExA

TransferEntity

nC(PS)

nC(PS)

FromFedSourceExAToFedDestE

FromFedSourceExBToFedDestE

(N/S)

TransferEntity

of the *EntityName* and *EntitySize* defined in ETS, respectively [51], [62].


**Table 3.** Fixed Record Datatype table

#### **3.3. Communication within the HLA-based integration infrastructure**

Pedrielli et al. [77] proposed a communication protocol (see Section 3.3.2) based on messages to manage the communication between a CSP and its middleware (or adapter). The communication protocol was conceived for the distributed simulation of network of Discrete Event Manufacturing Systems characterized by the transfer of parts in the presence of buffers with finite capacity, with the objective to minimize the use of zero-lookahead [62] for the synchronization of federates.

Before illustrating the communication protocol, Section 3.3.1. presents the concept and functioning of Simulation Messages created to support the communication between a CSP and the middleware. The communication protocol between federates is then explained in Section 3.3.2., whereas Section 3.3.3. defines the hypotheses needed to minimize the zero lookahead when applying the proposed protocol.

#### *3.3.1. Simulation messages*

The function of the simulation messages depends on the role played by the federate. The sending federate uses the message for communications concerning the need of sending an entity to another model (outgoing communication) and/or information on the availability of the target receiving federate (incoming communication). The receiving federate uses the message for communications concerning the buffer state and/or the acceptance of an entity (outgoing communication) and/or the receipt of an entity from other models (incoming communication). Simulation Messages are implemented as a class that is characterized by the following attributes:

 time referring to the simulation time when the message is sent to the middleware from the CSP. This attribute is used by the middleware to determine the TimeStamp of the interaction that will be sent to the RTI.

	- BoundedBuffer containing the information about the availability of the bounded buffer in the receiving model.

Distributed Modeling of Discrete Event Systems 23

the middleware that another federate needs to transfer an entity, the receiving model actually simulates the arrival of the entity only if the buffer is not full, otherwise the arrival

**Example**. The application of the Simulation Messages can be better appreciated by presenting an example (see Figure 9) that is characterized as follows: (1) the reference production system is represented in Figure 10, (2) the buffer Q2a at time t accommodates a number of parts that is greater than zero and less than the buffer capacity and an entity enters workstation W1a, (3) a departure event from workstation W1a is scheduled for time t' = t + p, where p represents the processing time of the leaving entity at station W1a, (4) during the time interval (t; t'), no event happening in the federate M2 (local event) influences the state of the buffer Q2a. Since W1a is the last machine in model M1, the departure event is also a *TransferEntityEvent*. Therefore, the CSP sends a message to its middleware containing time (t) and the *TransferEntityEvent* attributes. After receiving the

Once the RTI time advances to time t, the middleware of the receiving model receives the information about the need of the sending model to transfer an entity at time t'. Then, the middleware sends to the receiving model a simulation message containing the *ExternalArrivalEvent*. The receiving model simulates the external arrival as soon as the simulation time advances to t' and all local events for that time have been simulated (since

is not simulated and the workstation in the sending model becomes blocked.

message, the middleware of the sending model informs the RTI via interaction.

**Figure 9.** Communication Protocol

**Figure 10.** Reference Production System


### *3.3.2. Communication protocol*

Herein, the behavior of the sending federate will be analyzed at first, then the receiving federate will be taken under consideration. Finally an example will be described to clarify how the protocol works.

**Sending Federate**. The CSP of the sending federate sends a message to its middleware whenever a *TransferEntityEvent* is scheduled, i.e. the departure event of an entity from the last workstation of the sending model is added to the simulation event list. Then, the middleware uses the attributes time and *TransferEntityEvent* to inform the RTI about the need of passing an entity, while the simulation keeps on running (the *TransferEntityEvent* time corresponds to the *EntityTransferTime* presented in Section 3.2.).

The request to advance to *EntityTransferTime* is sent by the middleware to the RTI as soon as all local events scheduled for that time instant have been simulated.

After the time has advanced, the middleware can inform the CSP of the sending model about the state of the receiving buffer in the receiving model. If the receiving buffer is not full, then the workstation can simulate the *TransferEntityEvent*, otherwise it becomes blocked. From the blocking instant until when the middleware informs the sending model that the receiving buffer is not full, the model keeps on sending requests for time advance at the lookahead value.

**Receiving Federate**. The CSP of the receiving federate sends a message to its middleware whenever a change in the buffer availability occurs. This message contains the updated value of the attribute *boundedBuffer* representing the availability of the buffer, i.e. the number of available slots. Then, the middleware communicates this information to the RTI via interactions. In particular the information on the availability of the buffer represents a field of the timestamped interaction *transferEntityFromFedSourceExToFedDestEn* (Section 3.2.).

If the change in the buffer availability is due to the arrival of an entity from another model, then the update of the information does not imply zero lookahead and the communication is characterized by defining the entity that has been accepted (i.e. the *ReceivedEntity* attribute). If the buffer state change is not related to an external arrival, then the update of the buffer information may imply a zero lookahead whenever it is not possible to determine an advisable a-priori lookahead for the federation (Section 3.3.3) [62]. After being informed by the middleware that another federate needs to transfer an entity, the receiving model actually simulates the arrival of the entity only if the buffer is not full, otherwise the arrival is not simulated and the workstation in the sending model becomes blocked.

**Example**. The application of the Simulation Messages can be better appreciated by presenting an example (see Figure 9) that is characterized as follows: (1) the reference production system is represented in Figure 10, (2) the buffer Q2a at time t accommodates a number of parts that is greater than zero and less than the buffer capacity and an entity enters workstation W1a, (3) a departure event from workstation W1a is scheduled for time t' = t + p, where p represents the processing time of the leaving entity at station W1a, (4) during the time interval (t; t'), no event happening in the federate M2 (local event) influences the state of the buffer Q2a. Since W1a is the last machine in model M1, the departure event is also a *TransferEntityEvent*. Therefore, the CSP sends a message to its middleware containing time (t) and the *TransferEntityEvent* attributes. After receiving the message, the middleware of the sending model informs the RTI via interaction.

**Figure 9.** Communication Protocol

22 Discrete Event Simulations – Development and Applications

in the receiving model.

scheduled time for the event.

scheduled time for the event.

*3.3.2. Communication protocol* 

how the protocol works.

the lookahead value.

accepted by the receiving model.

BoundedBuffer containing the information about the availability of the bounded buffer

 TransferEntityEvent representing the entity transfer event scheduled in the sending model event list and contains the information about the entity to be transferred and the

 ExternalArrivalEvent representing the external arrival event that is scheduled in the receiving model. It contains the information about the entity to be received and the

ReceivedEntity representing the information about the entity that was eventually

Herein, the behavior of the sending federate will be analyzed at first, then the receiving federate will be taken under consideration. Finally an example will be described to clarify

**Sending Federate**. The CSP of the sending federate sends a message to its middleware whenever a *TransferEntityEvent* is scheduled, i.e. the departure event of an entity from the last workstation of the sending model is added to the simulation event list. Then, the middleware uses the attributes time and *TransferEntityEvent* to inform the RTI about the need of passing an entity, while the simulation keeps on running (the *TransferEntityEvent*

The request to advance to *EntityTransferTime* is sent by the middleware to the RTI as soon as

After the time has advanced, the middleware can inform the CSP of the sending model about the state of the receiving buffer in the receiving model. If the receiving buffer is not full, then the workstation can simulate the *TransferEntityEvent*, otherwise it becomes blocked. From the blocking instant until when the middleware informs the sending model that the receiving buffer is not full, the model keeps on sending requests for time advance at

**Receiving Federate**. The CSP of the receiving federate sends a message to its middleware whenever a change in the buffer availability occurs. This message contains the updated value of the attribute *boundedBuffer* representing the availability of the buffer, i.e. the number of available slots. Then, the middleware communicates this information to the RTI via interactions. In particular the information on the availability of the buffer represents a field of

If the change in the buffer availability is due to the arrival of an entity from another model, then the update of the information does not imply zero lookahead and the communication is characterized by defining the entity that has been accepted (i.e. the *ReceivedEntity* attribute). If the buffer state change is not related to an external arrival, then the update of the buffer information may imply a zero lookahead whenever it is not possible to determine an advisable a-priori lookahead for the federation (Section 3.3.3) [62]. After being informed by

the timestamped interaction *transferEntityFromFedSourceExToFedDestEn* (Section 3.2.).

time corresponds to the *EntityTransferTime* presented in Section 3.2.).

all local events scheduled for that time instant have been simulated.

**Figure 10.** Reference Production System

Once the RTI time advances to time t, the middleware of the receiving model receives the information about the need of the sending model to transfer an entity at time t'. Then, the middleware sends to the receiving model a simulation message containing the *ExternalArrivalEvent*. The receiving model simulates the external arrival as soon as the simulation time advances to t' and all local events for that time have been simulated (since

the buffer Q2a is not full according to the example settings). A message is sent to the middleware of the receiving model containing the updated level of Q2a (attribute *BoundedBuffer*) together with the information concerning the recently accommodated part (attribute *ReceivedEntity*).

Distributed Modeling of Discrete Event Systems 25

executed. When this happens it is possible to minimize the use of the zero lookahead for the

The federate sending the message can communicate with t < t' under the following

 The Entity transfer event is scheduled when the part enters the machine in the sending model. In this case the event is put into the event list a number of time units before it must be simulated that is at least equal to the processing time of the workstation under analysis. In the case the event is scheduled when the part leaves the workstation, then the condition holds if there exists a transfer time between the sending and the receiving model that is larger than zero and no events affect the arrival of the part once the transfer has started. The conditions aforementioned are not unrealistic when a manufacturing plant is simulated: both in the case the event is scheduled before or after the processing activity, the time between the departure from the exit point and the arrival to the entry point is in general not negligible. Nonetheless, in both the aforementioned cases, it is required that no other external events are scheduled by the same exit point during the interval (t; t'). This can happen when, after a leaving event has been scheduled, a failure affects the machine. In this case the information related to the part to be transferred has already been delivered and cannot be updated. As a consequence an external arrival event will be scheduled in the receiving model although the sending model will not be able to deliver the part because of the machine

 It is possible to communicate in advance the Buffer availability change event if the workstation processing the part schedules the leave event in advance to its realization and no other events are scheduled by the same workstation during the interval (t; t'). However, the zero lookahead cannot be avoided by the sending federate which cannot be aware of the downstream buffer changes and then it will send update request at the

 The zero lookahead can be avoided if the middleware of the receiving model can schedule the External Arrival event in advance and then inform the target federate(s) on the availability of the buffer in advance. This condition can be satisfied based on the

In the case one or more of the conditions aforementioned do not hold than the communication protocol shown in the Section 3.3.2. implies the use of zero lookahead. If the hypothesis that no additional external events must be scheduled by the same exit (entry) point in a federate (sending or receiving) within the time interval (t; t') is relaxed, then the middleware should be able to arrange incoming events in a queue and wait before delivering the information to the simulator until when the most updated information has been received. However it is quite straightforward to show that, in the worst case, the middleware should wait until when the simulation time reaches t', and therefore all the time advance requests would be performed at the zero lookahead. This relaxation is under

failure. A solution to this issue is part of present research.

communication between federates.

lookahead value.

analysis.

entity transfer event characteristics.

conditions:

Afterwards, the middleware sends two interactions to the RTI: one is with a TimeStamp equal to t' and contains the updated state of the buffer Q2a and the receipt of the entity, the other contains the request of time advance to time t'. Once the RTI reaches time t', the middleware of the sending model receives the information regarding the state of Q2a and the received entity by means of the RTI. Since the entity has been delivered to the receiving model, the station W1a is not blocked by the middleware.

#### *3.3.3. Formal characterization of the communication protocol*

This section defines which hypotheses are needed to minimize the occurrence of zero lookahead if the communication protocol afore presented is adopted.

Let represent an external event scheduled in the i-th federate j-th exit (entry) point at simulation time t, where t can be, in general, smaller or equal to t' that represents the simulation time when the event is supposed to be simulated. An event scheduled into the event list of a simulator is defined as external if one of the three following conditions holds:


Herein three types of external events are taken into consideration:


If t < t' it means that the simulation message can be sent by the sending (receiving) model and received by the target federate before the event contained in the message has to be executed. When this happens it is possible to minimize the use of the zero lookahead for the communication between federates.

24 Discrete Event Simulations – Development and Applications

model, the station W1a is not blocked by the middleware.

the queue of the receiving federate.

(W1a can be idle or blocked).

to a receiving federate.

of the External Arrival Event (Section 3.3.1.).

Herein three types of external events are taken into consideration:

*3.3.3. Formal characterization of the communication protocol* 

lookahead if the communication protocol afore presented is adopted.

(attribute *ReceivedEntity*).

the buffer Q2a is not full according to the example settings). A message is sent to the middleware of the receiving model containing the updated level of Q2a (attribute *BoundedBuffer*) together with the information concerning the recently accommodated part

Afterwards, the middleware sends two interactions to the RTI: one is with a TimeStamp equal to t' and contains the updated state of the buffer Q2a and the receipt of the entity, the other contains the request of time advance to time t'. Once the RTI reaches time t', the middleware of the sending model receives the information regarding the state of Q2a and the received entity by means of the RTI. Since the entity has been delivered to the receiving

This section defines which hypotheses are needed to minimize the occurrence of zero

Let represent an external event scheduled in the i-th federate j-th exit (entry) point at simulation time t, where t can be, in general, smaller or equal to t' that represents the simulation time when the event is supposed to be simulated. An event scheduled into the event list of a simulator is defined as external if one of the three following conditions holds: The realization of the event depends on the state of a federate that is, in general, different from the one that scheduled the event. One example of external event is when the sending federate (model M1) wants to transfer a part to the receiving federate (model M2), the possibility for the leaving event to be simulated depends on the state of

 The simulation of the event leads to changes into the state of other federates in the federation. This is the case when the downstream machine to the first buffer in the receiving model takes a part from the buffer thus changing its availability, this information must be delivered to sending models that are willing to transfer an entity, the state of the sending federate(s) will change depending on the information delivered

 The event is not scheduled by the simulator that will simulate it, but is put into the simulation event list by the middleware associated with the simulator. This is the case

*Entity transfer event*, this event happens when a sending federate wants to transfer a part

*Buffer\_availability* change event, this is a departure event from the workstation

*External Arrival event*, this event is scheduled by the middleware inside the simulation

If t < t' it means that the simulation message can be sent by the sending (receiving) model and received by the target federate before the event contained in the message has to be

downstream the buffer representing the entry point of the receiving model.

event list of the receiving federate every time a part has to be transferred.

The federate sending the message can communicate with t < t' under the following conditions:


In the case one or more of the conditions aforementioned do not hold than the communication protocol shown in the Section 3.3.2. implies the use of zero lookahead. If the hypothesis that no additional external events must be scheduled by the same exit (entry) point in a federate (sending or receiving) within the time interval (t; t') is relaxed, then the middleware should be able to arrange incoming events in a queue and wait before delivering the information to the simulator until when the most updated information has been received. However it is quite straightforward to show that, in the worst case, the middleware should wait until when the simulation time reaches t', and therefore all the time advance requests would be performed at the zero lookahead. This relaxation is under analysis.

## **4. Distributed simulation in industry: A case study**

This section presents an application of the architecture and communication protocol proposed in Section 3.3. The aim is to evaluate whether the use of distributed simulation can help to better analyze the dynamics of complex manufacturing systems, whereas the comparison between the HLA-based distributed simulation and a monolithic simulation in terms of computational efficiency is out of scope.

Distributed Modeling of Discrete Event Systems 27

between different roll types, since a roll can be sent to the Roll Shop depending on the

Even if separate simulators are developed, the Roll Shop designer still has to evaluate the performance of the whole system while taking into account the influence of the Roll Milling System related to (1) the arrival rate of worn out rolls from the Roll Milling System that is estimated from the yearly aggregate demand of reconditioned rolls and (2) the acceptance of

In addition to this the Roll Shop designer has to guarantee that the Roll Milling System the Roll Shop is being designed for, never waits for reconditioned rolls interrupting the

Hence, even if simulation models are available, usually the Roll Shop designer overdimensions the number of rolls that have to populate the whole system to avoid any waiting time at the Roll Milling System. For this reason we focused our first analysis on the effect of the number of rolls over the performance of the whole system (Roll Milling System

Sections 4.1 and 4.2 will give the main details characterizing the simulation models, whereas Section 4.3 will present (Section 4.3.1 and Section 4.3.2) and compare (Section 4.3.3) two

In this section the simulator of the Roll Shop developed for the case of interest will be

The Roll Shop simulator has been developed in C++ language using the object oriented paradigm. The C++ based simulator emulates a COTS, following the approach showed in Wang et al. [67], [68]. Figure 11 gives a pictorial representation of the simulation model for

*Buffer areas* where the rolls are kept while waiting to be transferred to the Roll Shop or

the reconditioned rolls sent by the Roll Shop (closed loop model).

behavior of other roll types.

production of sheet metal.

approaches for the system analysis.

**4.1. The Roll Shop simulator** 

the Roll Shop under analysis.

**Figure 11.** The Roll Shop System representation

the Milling system. The buffer areas can be:

The Roll Shop is composed by the following elements (Figure 11).

and Roll Shop).

explained in detail.

Herein the attention is focused on the industrial field represented by sheet metal production. In this industrial field, the production systems are characterized by the presence of at least two subsystems interacting with each other: the Roll Milling System and the Roll Shop. The Roll Milling System produces sheet metal using milling rolls that are subject to wearing out process. The rolls must be replaced when worn out to avoid defects in the laminates. Then the Roll Shop performs the grinding process to recondition the worn out rolls.

The following types of rolls have been considered in the case study:


If the attention is focused on the rolls, then the resulting production system is a closed loop: the Roll Milling System sends batches of worn out rolls to the Roll Shop following a given policy and receives reconditioned rolls back. Both the Roll Milling System and the Roll Shop have finite capacity buffers, therefore it is necessary to check whether the buffer in the system receiving the rolls has enough free slots. The deadlock in the closed loop is avoided because the number of rolls circulating in the system is less than the number of available slots (taking into account also the machines) and it is constant.

The two subsystems forming a closed loop are strongly related and their reciprocal influence should be considered to properly evaluate the performance of the whole factory by means of a comprehensive simulation model. However, both the Roll Shop designer and the Roll Milling System owner usually develop their own detailed simulator to evaluate the performance of their subsystem, because of the lack of shared information between the owner of the Roll Milling System and the Roll Shop designer. Indeed, the owner of the Roll Milling System usually provides the Roll Shop designer only with aggregated data about the yearly average demand of worn out rolls to be reconditioned. Moreover, the Roll Milling System works according to specific roll changing policies that are not shared with the Roll Shop designer even if they play a key role in the dynamics of the whole factory. For instance, when a roll is worn out, also the other rolls are checked and if their remaining duration is under a predefined threshold, then they are sent to the Roll Shop together with the completely worn out rolls. The presence of roll changing policies determines a relation between different roll types, since a roll can be sent to the Roll Shop depending on the behavior of other roll types.

Even if separate simulators are developed, the Roll Shop designer still has to evaluate the performance of the whole system while taking into account the influence of the Roll Milling System related to (1) the arrival rate of worn out rolls from the Roll Milling System that is estimated from the yearly aggregate demand of reconditioned rolls and (2) the acceptance of the reconditioned rolls sent by the Roll Shop (closed loop model).

In addition to this the Roll Shop designer has to guarantee that the Roll Milling System the Roll Shop is being designed for, never waits for reconditioned rolls interrupting the production of sheet metal.

Hence, even if simulation models are available, usually the Roll Shop designer overdimensions the number of rolls that have to populate the whole system to avoid any waiting time at the Roll Milling System. For this reason we focused our first analysis on the effect of the number of rolls over the performance of the whole system (Roll Milling System and Roll Shop).

Sections 4.1 and 4.2 will give the main details characterizing the simulation models, whereas Section 4.3 will present (Section 4.3.1 and Section 4.3.2) and compare (Section 4.3.3) two approaches for the system analysis.

## **4.1. The Roll Shop simulator**

26 Discrete Event Simulations – Development and Applications

terms of computational efficiency is out of scope.

bigger roll type) and IMR2.

rolls.

WR2.

**4. Distributed simulation in industry: A case study** 

The following types of rolls have been considered in the case study:

slots (taking into account also the machines) and it is constant.

This section presents an application of the architecture and communication protocol proposed in Section 3.3. The aim is to evaluate whether the use of distributed simulation can help to better analyze the dynamics of complex manufacturing systems, whereas the comparison between the HLA-based distributed simulation and a monolithic simulation in

Herein the attention is focused on the industrial field represented by sheet metal production. In this industrial field, the production systems are characterized by the presence of at least two subsystems interacting with each other: the Roll Milling System and the Roll Shop. The Roll Milling System produces sheet metal using milling rolls that are subject to wearing out process. The rolls must be replaced when worn out to avoid defects in the laminates. Then the Roll Shop performs the grinding process to recondition the worn out

 Intermediate Rolls (IMR) representing back-up rolls that are not in contact with the laminate. Depending on the size of the roll, the IMR roll will be referred to as IMR1 (the

 Work Roll (WR) representing the rolls directly in contact with the laminate. Also in this case there are two subtypes that will be referred to as WR1 (the bigger roll type) and

If the attention is focused on the rolls, then the resulting production system is a closed loop: the Roll Milling System sends batches of worn out rolls to the Roll Shop following a given policy and receives reconditioned rolls back. Both the Roll Milling System and the Roll Shop have finite capacity buffers, therefore it is necessary to check whether the buffer in the system receiving the rolls has enough free slots. The deadlock in the closed loop is avoided because the number of rolls circulating in the system is less than the number of available

The two subsystems forming a closed loop are strongly related and their reciprocal influence should be considered to properly evaluate the performance of the whole factory by means of a comprehensive simulation model. However, both the Roll Shop designer and the Roll Milling System owner usually develop their own detailed simulator to evaluate the performance of their subsystem, because of the lack of shared information between the owner of the Roll Milling System and the Roll Shop designer. Indeed, the owner of the Roll Milling System usually provides the Roll Shop designer only with aggregated data about the yearly average demand of worn out rolls to be reconditioned. Moreover, the Roll Milling System works according to specific roll changing policies that are not shared with the Roll Shop designer even if they play a key role in the dynamics of the whole factory. For instance, when a roll is worn out, also the other rolls are checked and if their remaining duration is under a predefined threshold, then they are sent to the Roll Shop together with the completely worn out rolls. The presence of roll changing policies determines a relation In this section the simulator of the Roll Shop developed for the case of interest will be explained in detail.

The Roll Shop simulator has been developed in C++ language using the object oriented paradigm. The C++ based simulator emulates a COTS, following the approach showed in Wang et al. [67], [68]. Figure 11 gives a pictorial representation of the simulation model for the Roll Shop under analysis.

**Figure 11.** The Roll Shop System representation

The Roll Shop is composed by the following elements (Figure 11).

 *Buffer areas* where the rolls are kept while waiting to be transferred to the Roll Shop or the Milling system. The buffer areas can be:

	- *Stand-by Area*, represents the entry/exit point of the Roll Shop and only the overhead crane can access it. The batches of rolls coming from the Roll Milling System and the batches of grinded rolls to send back to the Roll Milling System are placed here.

Distributed Modeling of Discrete Event Systems 29


[min/roll]

[min/roll]

[min/roll]


Component Statistics Unit of Measurement

Shop system for every roll

System Time for every roll

Grinding time for every roll

Transfer time for every roll

Number of processed rolls - Number of rolls in the buffer

Idle time [min] Utilization [%]

In this section the simulator of the Roll Milling System developed for the case of interest will

The Roll Milling System simulator has been developed in C++ language. The C++ based

A generic milling system can be represented as shown in Figure 12, whereas the simulation

simulator emulates a COTS, following the approach showed by Wang et al. [67], [68].

model realized for the case presented is graphically represented in Figure 13.

Waiting time in the buffer [min/roll]

Processing Time [min/roll]

Roll Number of rolls in the Roll

type

type

type

type

area

Buffer Number of rolls in the buffer -

**Table 4.** Roll Shop Statistics

be explained in detail.

**Figure 12.** Roll Milling System [87]

**4.2. The Roll Milling System simulator** 

Conveyor Transfer time [min/transfer]

Machine Utilization [%]

	- *Grinder Machine Work* (Grinder WR), i.e. grinding machine dedicated to *work roll* type.
	- *Grinder Machine Intermediate* (Grinder IMR), i.e. grinding machine dedicated to *intermediate roll* type.
	- *Electro-Discharge Texturing Machine* (EDT), i.e. machine executing a surface finishing process on the rolls.
	- *Loader*, i.e. an automatic conveyor that transfers rolls to the grinding machines
	- *Overhead crane*, i.e. a semi-automatic handling system to transfer rolls from the arrival point in the roll shop (Stand by area) to the exchange areas of the system

All the workstations, but the EDT, have two buffer positions (grey rectangles in Figure 11) within the working area.

The main parameters needed to configure the simulator are:


Table 4 summarizes the output statistics that can be gathered from the Roll Shop simulator. The minimum, maximum average values are supplied for every statistic, and the variance is computed as well.



*Workstations* where the rolls are reconditioned;

The main parameters needed to configure the simulator are:

*intermediate roll* type.

process on the rolls.

Two types of conveyors:

within the working area.

*Size* of the roll batches

workstations).

computed as well.

quantity as a function of the path.

type.

placed here.

type.

 *Stand-by Area*, represents the entry/exit point of the Roll Shop and only the overhead crane can access it. The batches of rolls coming from the Roll Milling System and the batches of grinded rolls to send back to the Roll Milling System are

*Exchange Area 1*, represents the interface between the part of the system managed

*Exchange Area 2*, represents the interface between the exit of the grinding system

*Grinder Machine Work* (Grinder WR), i.e. grinding machine dedicated to *work roll*

*Grinder Machine Intermediate* (Grinder IMR), i.e. grinding machine dedicated to

*Electro-Discharge Texturing Machine* (EDT), i.e. machine executing a surface finishing

 *Loader*, i.e. an automatic conveyor that transfers rolls to the grinding machines *Overhead crane*, i.e. a semi-automatic handling system to transfer rolls from the arrival point in the roll shop (Stand by area) to the exchange areas of the system

All the workstations, but the EDT, have two buffer positions (grey rectangles in Figure 11)

 *Process Sequence* for every roll, i.e. the process plan and the assignment of operations to the production resources in the Roll Shop. The process sequence depends on the roll

 *Processing time* of each operation. Those processing time have been considered deterministic for the experiments presented in this section (i.e. no failures affect the

*Transfer time, i.e.* the time to move the roll within the Roll Shop. This is a deterministic

 *Number of workstations* of each type. We have dedicated machines, in particular we have one grinding machine dedicated to WR type (i.e. WR1 and WR2) and one grinding machine dedicated to IMR type roll (i.e. IMR1 and IMR2). We have one EDT machine (processing only WR type rolls). In case the WR type roll finds the dedicated grinding machine occupied and the IMR machine is idle, then it can be processed also on the

Table 4 summarizes the output statistics that can be gathered from the Roll Shop simulator. The minimum, maximum average values are supplied for every statistic, and the variance is

*Type* of rolls (e.g. WR1 and WR2, IMR1 and IMR2 as described in Section 4)

IMR machine. However the time required for the process increases.

by the overhead-crane and the grinding system served by the loader

and the exit from the roll shop system managed by the overhead-crane.

#### **4.2. The Roll Milling System simulator**

In this section the simulator of the Roll Milling System developed for the case of interest will be explained in detail.

The Roll Milling System simulator has been developed in C++ language. The C++ based simulator emulates a COTS, following the approach showed by Wang et al. [67], [68].

A generic milling system can be represented as shown in Figure 12, whereas the simulation model realized for the case presented is graphically represented in Figure 13.

**Figure 12.** Roll Milling System [87]

Distributed Modeling of Discrete Event Systems 31



Component Statistics Unit of Measurement

Milling system for every roll

Number of rolls in the buffer

The Service Level (SL) is the typical key performance indicator (KPI) for analyzing the Roll Milling System. It is defined as the time the Roll Milling System produces sheet metal over the time it is operative (the total simulation time in the case of the computer experiment). It must be highlighted that the system cannot produce if all the required rolls are not present at each station. For this reason every station will have the same "Busy Time", i.e. the same time period during which it produces. The service level of a system can be increased managing the plant in a way such that we always have an available batch of rolls to change

It is then clear that the Service Level would be reduced if the Roll Milling System had to

The Roll Shop designer may choose two possible approaches to estimate the effects of the

 *Approach A*. The Roll Milling System is represented by a simplified model inside the detailed simulation model of the Roll Shop (Section 4.1). This simplified model roughly reproduces the Roll Milling System by generating the arrival of worn out rolls and

 *Approach B*. The performance of the whole factory is evaluated by adopting the HLAbased Infrastructure integrating the simulators of the Roll Milling System and of the

The simulator of the whole system is realized introducing within the detailed simulation model of the Roll Shop a simplified model of the Roll Milling System. Also this simulator is

Interchange Time [min]

Waiting Time [min] Busy Time [min]

Waiting time in the buffer [min/roll]

Roll Number of rolls in the Roll

type

Station Utilization [%]

area

Service Level Busy Time/Simulation Time [%]

Buffer Number of rolls in the buffer -

**Table 5.** Roll Milling System Simulator Statistics

accepting the reconditioned ones.

Roll Shop.

*4.3.1. Approach A* 

developed in C++ language.

the worn out ones that have to be sent to the Roll Shop.

wait too long for reconditioned rolls coming from the Roll Shop.

**4.3. Sheet metal production system analysis: Approach A, approach B**

Roll Shop design choices over the performance of the Roll Milling System:

**Figure 13.** Milling System Simulation Model

The Roll Milling System considered for the industrial case is composed of 5 stands (refer to Figure 12 for the definition of *stand*) each characterized by a milling station and a buffer for rolls (Figure 13). It is important to highlight that more rolls are needed with respect to those in use to process the metal sheet in order to minimize the waiting time of the milling system when the rolls are changed.

Once the rolls worn out, the Roll Milling System interrupts the process and the rolls must be replaced. The rolls are then replaced with the rolls of the same type available at the rolls buffer (Figure 13) close to the station. In the case the rolls required are not available the Roll Milling System stops producing. The worn out rolls are sent to the Roll Shop System as soon as a batch of rolls is ready (the size of the batch is usually fixed and represents a parameter of the simulation model). The batches are then transferred to the Roll Shop by means of a special conveyor, the *Transfer Car*.

The interval between roll changes (*interchange time*) is the interval between two consecutive sending of the same roll to the Roll Shop. This interval is fundamental to the correct sizing of the Roll Shop Plant (i.e. the number of grinding machines for every type, the size of the buffer areas) and of the number of rolls, for every type, populating the system.

This interval is mainly related to the life duration and the roll changing policies adopted within the system. For the industrial case considered, these policies can be brought back to two main criteria:


The main parameters to configure the simulation model are:


The Roll Milling System simulator supplies the following output statistics (Table 5):


**Table 5.** Roll Milling System Simulator Statistics

30 Discrete Event Simulations – Development and Applications

**Figure 13.** Milling System Simulation Model

when the rolls are changed.

special conveyor, the *Transfer Car*.

grinding process to avoid multiple sending.

Capacity of the buffers at every stand

 Number of rolls in the system Size of roll batches and types of rolls Life duration for each roll type

The main parameters to configure the simulation model are:

two main criteria:

on life.

Number of stands

The Roll Milling System considered for the industrial case is composed of 5 stands (refer to Figure 12 for the definition of *stand*) each characterized by a milling station and a buffer for rolls (Figure 13). It is important to highlight that more rolls are needed with respect to those in use to process the metal sheet in order to minimize the waiting time of the milling system

Once the rolls worn out, the Roll Milling System interrupts the process and the rolls must be replaced. The rolls are then replaced with the rolls of the same type available at the rolls buffer (Figure 13) close to the station. In the case the rolls required are not available the Roll Milling System stops producing. The worn out rolls are sent to the Roll Shop System as soon as a batch of rolls is ready (the size of the batch is usually fixed and represents a parameter of the simulation model). The batches are then transferred to the Roll Shop by means of a

The interval between roll changes (*interchange time*) is the interval between two consecutive sending of the same roll to the Roll Shop. This interval is fundamental to the correct sizing of the Roll Shop Plant (i.e. the number of grinding machines for every type, the size of the

This interval is mainly related to the life duration and the roll changing policies adopted within the system. For the industrial case considered, these policies can be brought back to

1. If an IMR roll has to be changed, since this requires high setup time, also the WR rolls from the same station are sent to the Roll Shop even if they have not reached their end

2. If more batches have almost reached their life duration, they are sent together to

The Roll Milling System simulator supplies the following output statistics (Table 5):

buffer areas) and of the number of rolls, for every type, populating the system.

The Service Level (SL) is the typical key performance indicator (KPI) for analyzing the Roll Milling System. It is defined as the time the Roll Milling System produces sheet metal over the time it is operative (the total simulation time in the case of the computer experiment). It must be highlighted that the system cannot produce if all the required rolls are not present at each station. For this reason every station will have the same "Busy Time", i.e. the same time period during which it produces. The service level of a system can be increased managing the plant in a way such that we always have an available batch of rolls to change the worn out ones that have to be sent to the Roll Shop.

It is then clear that the Service Level would be reduced if the Roll Milling System had to wait too long for reconditioned rolls coming from the Roll Shop.

#### **4.3. Sheet metal production system analysis: Approach A, approach B**

The Roll Shop designer may choose two possible approaches to estimate the effects of the Roll Shop design choices over the performance of the Roll Milling System:


#### *4.3.1. Approach A*

The simulator of the whole system is realized introducing within the detailed simulation model of the Roll Shop a simplified model of the Roll Milling System. Also this simulator is developed in C++ language.

More specifically, the Roll Milling System is modeled as an oracle sending (receiving) batches of worn out (reconditioned) rolls to the Roll Shop based on the information on the wearing out time of every roll type.

Distributed Modeling of Discrete Event Systems 33

dynamics of the Roll Milling System. A possible idea to increase the accuracy of Approach A could be to replace data in Table 6 with the estimated information in Table 7, thus taking

However, the missing feature of this approach is that in any case, even updating the life duration, we will not be able to reproduce the *dynamics* over time of the Roll Milling System which is what really affects the estimation of the performance of the whole system more

In addition, the only statistics of the Roll Milling System that can be gathered if Approach A

 Service Level. In particular, the utilization time (Section 4.2) is estimated as the difference between the total simulation time and the computed sum of time intervals during which the number of rolls for at least one type are equal to 0. If this condition is verified then the Roll Milling System cannot produce and has to wait for reconditioned rolls. As defined in Section 4.2, the ratio between this time and the total simulation time

Summarizing, Approach A supplies an approximate estimate of the whole system

 The real behavior of the Roll Milling System cannot be precisely modeled since it is reduced to a black box sending and receiving rolls (e.g. the roll changing policies are

The performance of the Roll Milling System cannot be evaluated in detail (e.g. mean

In Approach B, the detailed models of the two subsystems are directly adopted. Indeed the two simulators described in Section 4.1 and 4.2 are linked together thanks to the HLA-based

 MAK-RTI 3.3.2 (www.mak.com) was used as the RTI component implementation. The middleware was developed in C++ language following the specifications defined in

The *FederateAmbassador* and *RTIAmbassador* were provided by MAK-RTI as C++ classes and were linked to the *SimulationMiddleware*. Further extensions were needed to implement the proposed modification to ETS (Section 3.2) and the Simulation Messages (Section 3.3.1). The former required a modification to *FederateAmbassador* class, whereas the latter led to the development of a new C++ class. The *SimulationMiddleware* was implemented to manage the

starvation time for every station, mean level of roll buffers, etc).

The HLA-based architecture was implemented as follows:

Section 3.3 and was named *SimulationMiddleware*.

information contained in Simulation Messages.

into account, at least on average, the behavior of the Roll Milling System.

than the average behavior captured by data in Table 7.

 Average Number of Rolls at every stand Average Waiting Time for grinded rolls

gives the estimate of the service level.

is adopted are:

performance:

not modeled).

developed infrastructure (Figure 15).

*4.3.2. Approach B* 

It must be stressed that the simulator of Approach A cannot be considered as a proper monolithic simulator of the whole factory, since the Roll Milling System is only poorly modeled.

The input parameters needed to initialize the oracle are:


Table 6 defines the life duration of the rolls (for every type) given as input for Approach A..


**Table 6.** Rolls parameters

Table 7 reports the results in terms of average intervals between rolls change, estimated running the simulation model of the Roll Milling System (Section 4.2) as standalone. The life duration given as input to the Roll Milling System simulator was the same given in Table 6. Although the life duration used is the same the resulting intervals between rolls change are different because of the effect of the roll changing policies which are not taken into account in the simplified model of the Roll Milling System (Section 4.2).


**Table 7.** Rolls interarrival time estimated running the MS simulation model

In particular the intervals between rolls change result decreased because the WR2 type, characterized by the shorter life duration, draws the change of all other roll types. This result represents a first motivation towards the distributed approach. Indeed it shows that the aggregated information related to the interchange time is not enough to represent the dynamics of the Roll Milling System. A possible idea to increase the accuracy of Approach A could be to replace data in Table 6 with the estimated information in Table 7, thus taking into account, at least on average, the behavior of the Roll Milling System.

However, the missing feature of this approach is that in any case, even updating the life duration, we will not be able to reproduce the *dynamics* over time of the Roll Milling System which is what really affects the estimation of the performance of the whole system more than the average behavior captured by data in Table 7.

In addition, the only statistics of the Roll Milling System that can be gathered if Approach A is adopted are:


Summarizing, Approach A supplies an approximate estimate of the whole system performance:


#### *4.3.2. Approach B*

32 Discrete Event Simulations – Development and Applications

The input parameters needed to initialize the oracle are:

System, whereas the Roll Shop is empty.

in the simplified model of the Roll Milling System (Section 4.2).

Time to change rolls batch

Life duration of every roll type.

**Table 6.** Rolls parameters

in Section 4.1 and Section 4.2).

wearing out time of every roll type.

modeled.

More specifically, the Roll Milling System is modeled as an oracle sending (receiving) batches of worn out (reconditioned) rolls to the Roll Shop based on the information on the

It must be stressed that the simulator of Approach A cannot be considered as a proper monolithic simulator of the whole factory, since the Roll Milling System is only poorly

 Number of rolls present in the system, when the simulation starts, for every type. In this case the number of rolls represents the total number of rolls in the Roll Milling System. The only initial condition we can set using this simulation model is that all the rolls at the simulation start are at the beginning of their life and are all in the Roll Milling

Size of the batch of rolls for every roll type (the rolls are moved in batches as explained

Table 6 defines the life duration of the rolls (for every type) given as input for Approach A..

Roll Type WR1 WR2 IMR1 IMR2 Batch Size [#Rolls/batch] 8 2 8 2 Duration [min/batch] 438 288 3228 1452

Table 7 reports the results in terms of average intervals between rolls change, estimated running the simulation model of the Roll Milling System (Section 4.2) as standalone. The life duration given as input to the Roll Milling System simulator was the same given in Table 6. Although the life duration used is the same the resulting intervals between rolls change are different because of the effect of the roll changing policies which are not taken into account

[min/batch] Detailed Model

In particular the intervals between rolls change result decreased because the WR2 type, characterized by the shorter life duration, draws the change of all other roll types. This result represents a first motivation towards the distributed approach. Indeed it shows that the aggregated information related to the interchange time is not enough to represent the

WR1 432 WR2 288 IMR1 3024.04 IMR2 1440

**Table 7.** Rolls interarrival time estimated running the MS simulation model

In Approach B, the detailed models of the two subsystems are directly adopted. Indeed the two simulators described in Section 4.1 and 4.2 are linked together thanks to the HLA-based developed infrastructure (Figure 15).

The HLA-based architecture was implemented as follows:


The *FederateAmbassador* and *RTIAmbassador* were provided by MAK-RTI as C++ classes and were linked to the *SimulationMiddleware*. Further extensions were needed to implement the proposed modification to ETS (Section 3.2) and the Simulation Messages (Section 3.3.1). The former required a modification to *FederateAmbassador* class, whereas the latter led to the development of a new C++ class. The *SimulationMiddleware* was implemented to manage the information contained in Simulation Messages.

The interaction tables (Section 3.2) developed for this case are shown in Tables 10,11,12 and 13.

Distributed Modeling of Discrete Event Systems 35

*difference* 

average statistics related to the Roll Milling System from the simulation of the High level

High Level 0.995 0.872 12.3 Medium Level 0.946 0.682 27.7 Low Level 0.308 0.273 3.5

> **Output Statistic Value** Busy Time 5.2[months] Waiting Time 0.8[months] Utilization 87.2[%]

Number of Rolls

**Table 9.** Roll Milling System simulator. Output Statistics

effectiveness as supporting tool for decision making.

Service Level 87.2[%]

estimation due to this additional information that characterizes the model.

The number of rolls in the system together with the roll changing policy have a strong impact on the workload conditions of the Roll Shop. This aspect can only be taken into account under Approach B and the results show the dramatic difference in performance

Approach A is overestimating the performance of the system, thus decreasing its

Based on the analysis carried out so far, it was decided to design further experiments to analyze the behavior of the whole system with different starting workload conditions, i.e.

*Approach A Approach B Percentage* 

Stand 1: 1.177 WR 1.470 IMR Stand 2: 1.177 WR 1.470 IMR Stand 3: 1.177 WR 1.470 IMR Stand 4: 1.177 WR 1.470 IMR Stand 5: 3.401 WR 1.740 IMR

condition in Table 8.

**Table 8.** Service Level Results

Experimental Conditions

**Figure 14.** Distributed Simulation of Roll Shop and Roll Milling System

#### *4.3.3. Comparison between approach A and approach B*

The two approaches have been compared by designing a set of experiments characterized as follows:


The results of the experiments are shown in Table 8. Approach A and Approach B are compared in terms of the estimated Service Level (refer to Section 4.2 for the definition). The results show that the difference between the two approaches is larger for the High and Medium level conditions.

When the level of rolls is Low the roll changing policy does not affect the overall performance of the production system because the Roll Milling System is frequently starved and therefore the estimations are similar (consider also that the Low level condition has no industrial meaning, but was considered to study the service level response). In case of Medium and High level conditions the workload of the rolls in the Roll Shop can be strongly influenced by the roll changing policy, thus generating a higher difference in the estimation between the two approaches.

Approach B generates more accurate estimates of the whole system performance because the Roll Milling System is represented with high level of detail and the roll changing policies are modeled. In addition to this if Approach B is used, detailed information related to the Roll Milling System performance are available. Table 9 shows an example of output for the average statistics related to the Roll Milling System from the simulation of the High level condition in Table 8.


**Table 8.** Service Level Results

34 Discrete Event Simulations – Development and Applications

**Figure 14.** Distributed Simulation of Roll Shop and Roll Milling System

*4.3.3. Comparison between approach A and approach B* 

13.

follows:

High level.

the experimentation.

Medium level conditions.

between the two approaches.

The interaction tables (Section 3.2) developed for this case are shown in Tables 10,11,12 and

The two approaches have been compared by designing a set of experiments characterized as

 Three experimental conditions were designed with reference to the total number of rolls circulating in the whole system. These three conditions were defined as Low, Medium,

 The simulation run length was set to six months (4 weeks of transitory period). The roll changing policy adopted for the Approach B simulator has been kept fixed throughout

The results of the experiments are shown in Table 8. Approach A and Approach B are compared in terms of the estimated Service Level (refer to Section 4.2 for the definition). The results show that the difference between the two approaches is larger for the High and

When the level of rolls is Low the roll changing policy does not affect the overall performance of the production system because the Roll Milling System is frequently starved and therefore the estimations are similar (consider also that the Low level condition has no industrial meaning, but was considered to study the service level response). In case of Medium and High level conditions the workload of the rolls in the Roll Shop can be strongly influenced by the roll changing policy, thus generating a higher difference in the estimation

Approach B generates more accurate estimates of the whole system performance because the Roll Milling System is represented with high level of detail and the roll changing policies are modeled. In addition to this if Approach B is used, detailed information related to the Roll Milling System performance are available. Table 9 shows an example of output for the


**Table 9.** Roll Milling System simulator. Output Statistics

The number of rolls in the system together with the roll changing policy have a strong impact on the workload conditions of the Roll Shop. This aspect can only be taken into account under Approach B and the results show the dramatic difference in performance estimation due to this additional information that characterizes the model.

Approach A is overestimating the performance of the system, thus decreasing its effectiveness as supporting tool for decision making.

Based on the analysis carried out so far, it was decided to design further experiments to analyze the behavior of the whole system with different starting workload conditions, i.e.

the number of rolls that are present in the Roll Milling System when the simulation starts. These experiments can be useful to analyze the ramp-up period and select the roll changing policy that avoids the arising of critical workload conditions. These additional experiments can be carried out only adopting *Approach B*, since the starting workload conditions cannot be modeled with *Approach A* (Section 4.3.1). Indeed, the simplified model generates rolls for the Roll Shop independently from the starting workload conditions. Therefore, the simplified model would generate roll arrivals even if all the rolls are already located in the Roll Shop, thus incorrectly increasing the number of rolls in the whole system.

Distributed Modeling of Discrete Event Systems 37

HLAreliable TimeStamp

HLAreliable TimeStamp

TransferEntityFromFedSourceExMillingtoFe

dDestEnRollShop(P/S)

RollEntity RollType HLAreliable TimeStamp

ReceivedRollEntity RollType HLAreliable TimeStamp

RollShopBuffer HLAInteger32BE HLAreliable TimeStamp

SourcePriority HLAInteger32BE HLAreliable TimeStamp

EntityTransferTime HLAFloat32BE HLAreliable TimeStamp

RollEntity RollType HLAreliable TimeStamp ReceivedRollEntity RollType HLAreliable TimeStamp

EntityTransferTime HLAFloat32BE HLAreliable TimeStamp

Interaction Parameter DataType Transportation Order

Interaction Parameter DataType Transportation Order

MillingBuffer HLAInteger32B

SourcePriority HLAInteger32B

E

E

HLAinter action Root(N)

Entity

TransferEntityFromF edSourceExMillingto FedDestEnRollShop(

TransferEntityFromFe dSourceExRollMillingt oFedDestEnRollShop

P/S)

(P/S)

Transfer Entity (N/S)

(N/S)

**Table 10.** Interaction Class Table

**Table 11.** Parameter Table for Roll Shop SOM

**Table 12.** Parameter Table for Milling System SOM

RollType Roll HLAASCIIString

BatchSize HLAInteger32BE

Record Name Field Encoding

Name Type HLAfixedRecord

**Table 13.** Fixed Record Datatype Table (SOM) for Roll Shop and Milling Systems

TransferEntityFromFedSource

TransferEntityFromFedSource

Milling(N/S)

RollShop(N/S)

The second set of experiments was designed as follows:


Fig. 15 shows the main effects plot for the *Service Level* evaluated by simulating the 27 resulting experimental conditions with *Approach B*. The plot suggests a significant influence of the factor *Starting Workload for WR*. This roll type assumes a key role because of its short roll life (Section 4.3.1). If the *Starting Workload For WR* is *Low*, the Roll Shop can hardly follow the frequent roll requests of *WR* from the Roll Milling System and low values of *SL* are observed. This phenomenon occurs in all conditions of the simulation lengths, however it mitigates when the simulation length increases. Indeed the *SL* tends to a stationary value that is independent from the starting conditions. Nonetheless this analysis can be useful for the Roll Milling System owner that can individuate critical conditions, thus designing roll changing policies that avoid the occurrence of these situations during the ramp-up period.

**Figure 15.** Main Effects Plot of Service Lovel, Approach B


**Table 10.** Interaction Class Table

36 Discrete Event Simulations – Development and Applications

The second set of experiments was designed as follows:

The roll changing policy is fixed for all experiments

roll life than *WR*

Milling System are considered

fixed for all the experiments.

0,9 0,8 0,7 0,6 0,5

0,9 0,8 0,7 0,6 0,5

**Figure 15.** Main Effects Plot of Service Lovel, Approach B

**Mean Service Levels**

the number of rolls that are present in the Roll Milling System when the simulation starts. These experiments can be useful to analyze the ramp-up period and select the roll changing policy that avoids the arising of critical workload conditions. These additional experiments can be carried out only adopting *Approach B*, since the starting workload conditions cannot be modeled with *Approach A* (Section 4.3.1). Indeed, the simplified model generates rolls for the Roll Shop independently from the starting workload conditions. Therefore, the simplified model would generate roll arrivals even if all the rolls are already located in the

Two *types of roll* circulate in the factory (*WR* and *IMR*). The roll of type *IMR* has a longer

For each type of roll three levels of the *Starting Workload* (i.e. number of rolls) in the Roll

total number of rolls is equal to the *High* level of the previous experimentation and is

Fig. 15 shows the main effects plot for the *Service Level* evaluated by simulating the 27 resulting experimental conditions with *Approach B*. The plot suggests a significant influence of the factor *Starting Workload for WR*. This roll type assumes a key role because of its short roll life (Section 4.3.1). If the *Starting Workload For WR* is *Low*, the Roll Shop can hardly follow the frequent roll requests of *WR* from the Roll Milling System and low values of *SL* are observed. This phenomenon occurs in all conditions of the simulation lengths, however it mitigates when the simulation length increases. Indeed the *SL* tends to a stationary value that is independent from the starting conditions. Nonetheless this analysis can be useful for the Roll Milling System owner that can individuate critical conditions, thus designing roll changing policies that avoid the occurrence of these situations during the ramp-up period.

**Main Effects Plot for Service Level**

Low Medium High

WR

Three simulation run lengths are considered, i.e. 1 week, 2 weeks and 4 weeks

Low Medium High

Simulation Run Length

IMR

1[week] 2[weeks] 4[weeks]

Roll Shop, thus incorrectly increasing the number of rolls in the whole system.


**Table 11.** Parameter Table for Roll Shop SOM


**Table 12.** Parameter Table for Milling System SOM


**Table 13.** Fixed Record Datatype Table (SOM) for Roll Shop and Milling Systems

## **5. Conclusions**

This chapter has presented an overview of distributed simulation and the contemporary innovations in the use of distributed modeling to support the analysis of complex systems. The attention has been focused on CSP-based distributed simulation in civilian applications and especially in manufacturing domain. The literature review showed the need of a general standard solution to the distributed simulation of systems within the civilian domain to increase the use of distributed techniques for the analysis of complex systems. The Interoperability Reference Model standard released by SISO CSPI-PDG has been analyzed. Furthermore, the latest advancements in ETS standard proposal and Communication Protocol between federates within HLA-based distributed environment have been shown. In the end the industrial application of an HLA-based infrastructure proved the benefits of the distributed approach to effectively analyze the behavior of complex industrial systems.

Distributed Modeling of Discrete Event Systems 39

interoperability between digital tools. Future research will evaluate if the VFF approach can

*Istituto Tecnologie Industriali e Automazione (ITIA), Consiglio Nazionale delle Ricerche (CNR),* 

The research reported in this chapter has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement No: NMP2 2010- 228595, Virtual Factory Framework (VFF). Sections of this chapter are based on the paper: Pedrielli, G., Sacco, M., Terkaj, W., Tolio, T., Simulation of complex manufacturing systems via HLA-based infrastructure. *Journal Of Simulation, to be published*. Authors would like to acknowledge Professor Taylor S.J.E and Professor Starssburger S. for their ongoing support and their work the in developing current range of standards for distributed simulation. Authors acknowledge Tenova Pomini Company for providing the support to build the

[1] Ieee standard for modeling and simulation (m amp;s) high level architecture (hla) -

[2] Khaldoon Al-Zoubi and Gabriel A. Wainer. Performing distributed simulation with

[3] James G. Anderson. Simulation in the health services and biomedicine, pages 275{293.

[4] J. Banks, J.C. Hugan, P. Lendermann, C. McLean, E.H. Page, C.D. Pegden, O. Ulgen, and J.R. Wilson. The future of the simulation industry. In Simulation Conference, 2003.

[6] Lee A. Belfore, II. Simulation in environmental and ecological systems, pages 295-314.

[7] S. Bergmann, S. Stelzer, and S Strassburger. Initialization of simulation models using cmsd. In Simulation Conference (WSC), Proceedings of the 2011 Winter, WSC '11, pages

[8] C.A. Boer. Distributed Simulation in Industry. PhD thesis, Erasmus University

Proceedings of the 2003 Winter, volume 2, pages 2033-2043 vol.2, dec. 2003. [5] J. Banks and B.L. Nelson. Discrete-event system simulation. Prentice Hall, 2010.

framework and rules. IEEE Std. 1516-2000, pages i-22, 2000.

Kluwer Academic Publishers, Norwell, MA, USA, 2003.

Kluwer Academic Publishers, Norwell, MA, USA, 2003.

Rotterdam, Rotterdam, The Netherlands, October 2005.

*Politecnico di Milano, Dipartimento di Ingegneria Meccanica, Milano (MI), Italy* 

be exploited by distributed simulation applications.

**Author details** 

*Milano (MI) Italy* 

industrial case.

**6. References** 

restful web services approach.

594-602. IEEE Press, 2011.

**Acknowledgement** 

Giulia Pedrielli and Tullio Tolio

Walter Terkaj and Marco Sacco

The distribute simulation is acquiring more and more interest also because the need to interoperate several simulators leads to the need to improve the methodologies and the tools developed so far for the simulation of Discrete Event Systems. In other words the issues arising when trying to make several simulators interoperate are those issues the simulation community has been dealing with in the last years.

Further effort is needed in the formalization of the IRMs. In particular the presence of shared resources and the modeling of the control policies characterizing the system represent challenging issues not yet solved.

The communication protocols need to be enhanced and the zero-lookahead issue represents one of the main bottlenecks against the increase of efficiency of distributed simulations. Research effort is necessary to come up with new algorithms enabling the avoidance of zerolookahead. In this area the research on protocols that do not force to send interaction at every time unit to communicate the state of the federates, but enable the interaction depending on the system state (Adaptive Communication Protocols) look promising.

The need of a shared data model ([26], [7]) and of a common definition of the objects which are input and output of the simulation and a common simulation language, are all scientific and technical challenging topics that make research in distributed simulation always up to date.

In particular the definition of a common reference model to describe information generated by each simulator while it is running will be a key factor for the success of the distributed simulation technique. A fundamental contribution in this field was given by the Core Manufacturing Simulation Data CMSD [7], but the interoperability between simulators is still far from being reached. The data modeling research topic is wider than what just stated and covers several areas and simulation is just one of those. For example proposals to standardize the information modeling for manufacturing systems have been done in [83-84].

The European project Virtual Factory Framework VFF ([38] [85-86]) represents one of the latest proposals from the research community in terms of framework supporting the interoperability between digital tools. Future research will evaluate if the VFF approach can be exploited by distributed simulation applications.

## **Author details**

38 Discrete Event Simulations – Development and Applications

This chapter has presented an overview of distributed simulation and the contemporary innovations in the use of distributed modeling to support the analysis of complex systems. The attention has been focused on CSP-based distributed simulation in civilian applications and especially in manufacturing domain. The literature review showed the need of a general standard solution to the distributed simulation of systems within the civilian domain to increase the use of distributed techniques for the analysis of complex systems. The Interoperability Reference Model standard released by SISO CSPI-PDG has been analyzed. Furthermore, the latest advancements in ETS standard proposal and Communication Protocol between federates within HLA-based distributed environment have been shown. In the end the industrial application of an HLA-based infrastructure proved the benefits of the distributed approach to effectively analyze the behavior of complex industrial systems.

The distribute simulation is acquiring more and more interest also because the need to interoperate several simulators leads to the need to improve the methodologies and the tools developed so far for the simulation of Discrete Event Systems. In other words the issues arising when trying to make several simulators interoperate are those issues the

Further effort is needed in the formalization of the IRMs. In particular the presence of shared resources and the modeling of the control policies characterizing the system

The communication protocols need to be enhanced and the zero-lookahead issue represents one of the main bottlenecks against the increase of efficiency of distributed simulations. Research effort is necessary to come up with new algorithms enabling the avoidance of zerolookahead. In this area the research on protocols that do not force to send interaction at every time unit to communicate the state of the federates, but enable the interaction

The need of a shared data model ([26], [7]) and of a common definition of the objects which are input and output of the simulation and a common simulation language, are all scientific and technical challenging topics that make research in distributed simulation always up to

In particular the definition of a common reference model to describe information generated by each simulator while it is running will be a key factor for the success of the distributed simulation technique. A fundamental contribution in this field was given by the Core Manufacturing Simulation Data CMSD [7], but the interoperability between simulators is still far from being reached. The data modeling research topic is wider than what just stated and covers several areas and simulation is just one of those. For example proposals to standardize the information modeling for manufacturing systems have been done in [83-84]. The European project Virtual Factory Framework VFF ([38] [85-86]) represents one of the latest proposals from the research community in terms of framework supporting the

depending on the system state (Adaptive Communication Protocols) look promising.

simulation community has been dealing with in the last years.

represent challenging issues not yet solved.

date.

**5. Conclusions** 

Giulia Pedrielli and Tullio Tolio *Politecnico di Milano, Dipartimento di Ingegneria Meccanica, Milano (MI), Italy* 

Walter Terkaj and Marco Sacco *Istituto Tecnologie Industriali e Automazione (ITIA), Consiglio Nazionale delle Ricerche (CNR), Milano (MI) Italy* 

## **Acknowledgement**

The research reported in this chapter has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement No: NMP2 2010- 228595, Virtual Factory Framework (VFF). Sections of this chapter are based on the paper: Pedrielli, G., Sacco, M., Terkaj, W., Tolio, T., Simulation of complex manufacturing systems via HLA-based infrastructure. *Journal Of Simulation, to be published*. Authors would like to acknowledge Professor Taylor S.J.E and Professor Starssburger S. for their ongoing support and their work the in developing current range of standards for distributed simulation. Authors acknowledge Tenova Pomini Company for providing the support to build the industrial case.

### **6. References**

	- [9] C.A. Boer, A. de Bruin, and A. Verbraeck. Distributed simulation in industry a survey part 1 - the cots vendors. volume 0, pages 1053-1060, Los Alamitos, CA, USA, 2006. IEEE Computer Society.

Distributed Modeling of Discrete Event Systems 41

distributed decision support system for a semiconductor supply chain. SimTech

[23] Sumit Ghosh and Tony Lee. Modeling and Asynchronous Distributed Simulation

[24] Strong D.R. Richards N. Goel N.C. Goel, S. A simulation-based method for the process to allow continuous tracking of quality, cost, and time. Simulation, 78(5):330-337, 2001. [25] Hironori Hibino, Yoshiyuki Yura, Yoshiro Fukuda, Keiji Mitsuyuki, and Kiyoshi Kaneda. Manufacturing modeling architectures: manufacturing adapter of distributed simulation systems using hla. In Proceedings of the 34th conference on Winter simulation: exploring new frontiers, WSC '02, pages 1099-1107. Winter Simulation

[26] Marcus Johansson, Bjorn Johansson, Anders Skoogh, Swee Leong, Frank Riddick, Y. Tina Lee, Guodong Shao, and Par Klingstam. A test implementation of the core manufacturing simulation data specification. In Proceedings of the 39th conference on Winter simulation: 40 years! The best is yet to come, WSC '07, pages 1673-1681,

[27] Frederick Kuhl, Richard Weatherly, and Judith Dahmann. Creating computer simulation systems: an introduction to the high level architecture. Prentice Hall PTR,

[28] P. Lendermann. About the need for distributed simulation technology for the resolution of realworld manufacturing and logistics problems. In Simulation Conference, 2006.

[29] P. Lendermann, M.U. Heinicke, L.F. McGinnis, C. McLean, S. Strassburger, and S.J.E. Taylor. Panel: distributed simulation in industry - a real-world necessity or ivory tower

[30] Richard J. Linn, Chin-Sheng Chen, and Jorge A. Lozan. Manufacturing supply chain applications 2: development of distributed simulation model for the transporter entity in a supply chain process. In Proceedings of the 34th conference on Winter simulation: exploring new frontiers, WSC '02, pages 1319-1326. Winter Simulation

[31] A. W. Malik. An optimistic parallel simulation for cloud computing environments. SCS

[32] Charles McLean and Swee Leong. The expanding role of simulation in future manufacturing. In Proceedings of the 33nd conference on Winter simulation, WSC '01,

[33] Charles McLean, Swee Leong, Charley Harrell, Philomena M. Zimmerman, and Roberto F. Lu. Simulation standards: current status, needs, future directions, panel: simulation standards: current status, needs, and future directions. In Proceedings of the 35th conference on Winter simulation: driving innovation, WSC '03, pages 2019-2026. Winter

pages 1478-1486, Washington, DC, USA, 2001. IEEE Computer Society.

fancy? In Simulation Conference, 2007 Winter, pages 1053-1062, dec. 2007.

WSC 06. Proceedings of the Winter, pages 1119 -1128, dec. 2006.

Analyzing Complex Systems. Wiley-IEEE Press, 1st edition, 2000.

technical reports, 7(4):220-226, 2007.

Piscataway, NJ, USA, 2007. IEEE Press.

Upper Saddle River, NJ, USA, 1999.

Conference, 2002.

Conference, 2002.

M&S Magazine, 6:1-9, 2010.

Simulation Conference, 2003.


distributed decision support system for a semiconductor supply chain. SimTech technical reports, 7(4):220-226, 2007.

[23] Sumit Ghosh and Tony Lee. Modeling and Asynchronous Distributed Simulation Analyzing Complex Systems. Wiley-IEEE Press, 1st edition, 2000.

40 Discrete Event Simulations – Development and Applications

IEEE Computer Society.

pages 1094-1102, dec. 2008.

Simulation Conference, 2003.

SIMULATION, 71(6):378- 387, 1998.

Publishers, Norwell, MA, USA, 2003.

distributed computing. Wiley, 2000.

pages 1412 -1420 vol.2, 2001.

Conference, 2006.

[9] C.A. Boer, A. de Bruin, and A. Verbraeck. Distributed simulation in industry - a survey part 1 - the cots vendors. volume 0, pages 1053-1060, Los Alamitos, CA, USA, 2006.

[10] C.A. Boer, A. de Bruin, and A. Verbraeck. Distributed simulation in industry - a survey part 3 - the hla standard in industry. In Simulation Conference, 2008. WSC 2008. Winter,

[11] Csaba Attila Boer, Arie de Bruin, and Alexander Verbraeck. Distributed simulation in industry - a survey: part 2 - experts on distributed simulation. In Proceedings of the 38th conference on Winter simulation, WSC '06, pages 1061-1068. Winter Simulation

[12] Csaba Attila Boer and Alexander Verbraeck. Distributed simulation and manufacturing: distributed simulation with cots simulation packages. In Proceedings of the 35th conference on Winter simulation: driving innovation, WSC '03, pages 829-837. Winter

[13] Vesna Bosilj-Vuksic, Mojca Indihar Stemberger, Jurij Jaklic, and Andrej Kovacic. Assessment of e-business transformation using simulation modeling. Simulation. [14] Agostino G. Bruzzone. Preface to modeling and simulation methodologies for logistics

[15] Judith S. Dahmann, Frederick Kuhl, and Richard Weatherly. Standards for simulation: As simple as possible but not simpler the high level architecture for simulation.

[16] Wilhelm Dangelmaier and Bengt Mueck. Simulation in business administration and management, pages 391-406. Kluwer Academic Publishers, Norwell, MA, USA, 2003. [17] Paul K. Davis. Military applications of simulation, pages 407-435. Kluwer Academic

[18] T. Eldabi and R.J. Paul. A proposed approach for modeling healthcare systems for understanding. In Simulation Conference, 2001. Proceedings of the Winter, volume 2,

[19] Richard Fujimoto, Michael Hunter, Jason Sirichoke, Mahesh Palekar, Hoe Kim, and Wonho Suh. Ad hoc distributed simulations. In Proceedings of the 21st International Workshop on Principles of Advanced and Distributed Simulation, PADS '07, pages 15-

[20] R.M. Fujimoto. Parallel and distributed simulation systems. Wiley series on parallel and

[21] Boon Ping Gan, Peter Lendermann, Malcolm Yoke Hean Low, Stephen J. Turner, Xiaoguang Wang, and Simon J. E. Taylor. Interoperating autosched ap using the high level architecture. In Proceedings of the 37th conference on Winter simulation, WSC '05,

[22] Boon Ping Gan, Peter Lendermann, Malcolm Yoke Hean Low, Stephen J. Turner, Xiaoguang Wang, and Simon J. E. Taylor. Architecture and performance of an hla-based

and manufacturing optimization. Simulation, 80(3):119-120, 2004.

24, Washington, DC, USA, 2007. IEEE Computer Society.

pages 394-401. Winter Simulation Conference, 2005.

	- [34] Navonil Mustafee, Simon J.E. Taylor, Korina Katsaliaki, and Sally Brailsford. Facilitating the analysis of a uk national blood service supply chain using distributed simulation. SIMULATION, 85(2):113-128, 2009.

Distributed Modeling of Discrete Event Systems 43

[48] S. Straburger, G. Schmidgall, and S. Haasis. Distributed manufacturing simulation as an enabling technology for the digital factory. Journal of Advanced ManufacturingSystem,

[49] S. Strassburger. Distributed Simulation Based on the High Level Architecture in Civilian Application Domains. PhD thesis, Computer Science Faculty, University Otto

[50] Steffen Strassburger. The road to cots-interoperability: from generic hla-interfaces towards plug and play capabilities. In Proceedings of the 38th conference on Winter

[51] Taylor, Strassburger, S.J. Turner, M.Y.H. Low, Xiaoguang Wang, and J. Ladbrook. Developing interoperability standards for distributed simulaton and cots simulation packages with the cspi pdg. In Simulation Conference, 2006. WSC 06. Proceedings of the

[52] S J E Taylor and N Mustafee. An analysis of internal/external event ordering strategies for cots distributed simulation. pages 193-198. Proceedings of the 15th European

[53] Simon J. E. Taylor, Navonil Mustafee, Steffen Strassburger, Stephen J. Turner, Malcolm Y. H. Low, and John Ladbrook. The siso cspi pdg standard for commercial off-the-shelf simulation package interoperability reference models. In Proceedings of the 39th conference on Winter simulation: 40 years! The best is yet to come, WSC '07, pages 594-

[54] Simon J. E. Taylor, Stephen J. Turner, and Steffen Strassburger. Guidelines for commercial-off-the-shelf simulation package interoperability. In Proceedings of the 40th Conference on Winter Simulation, WSC '08, pages 193-204. Winter Simulation

[55] Simon J.E. Taylor. A proposal for an entity transfer specification standard for cots simulation package interoperation. In Proceedings of the 2004 European Simulation

[56] Simon J.E. Taylor. Distributed Modeling, pages 9:1-9:19. Chapman & Hall/CRC, 2007. [57] S.J.E. Taylor. Realising parallel and distributed simulation in industry: A roadmap. In Principles of Advanced and Distributed Simulation (PADS), 2011 IEEE Workshop on,

[58] S.J.E. Taylor, A. Bruzzone, R. Fujimoto, Boon Ping Gan, S. Strassburger, and R.J. Paul. Distributed simulation and industry: potentials and pitfalls. In Simulation Conference,

[59] S.J.E Taylor, M. Ghorbani, N. Mustafee, S. J. Turner, T. Kiss, D. Farkas, S. Kite, and S. Strassburger. Distributed computing and modeling & simulation: speeding up simulations and creating large models. In Proceedings of the 2011 Winter Simulation

[60] S.J.E. Taylor, M. Ghorbani, N. Mustafee, S.J. Turner, T. Kiss, D. Farkas, S. Kite, and S. Strassburger. Distributed computing and modeling amp; simulation: Speeding up

2002. Proceedings of the Winter, volume 1, pages 688- 694 vol.1, dec. 2002.

simulation, WSC '06, pages 1111-1118. Winter Simulation Conference, 2006.

2(1):111 -126, 2003.

Conference, 2008.

page 1, june 2011.

Interoperability Workshop, 2004.

Conference, pages 1-15, December 2011.

von Guericke, Magdeburg, 2001.

Winter, pages 1101 -1110, dec. 2006.

Simulation Symposium (ESS2003), Delft, 2003.

602, Piscataway, NJ, USA, 2007. IEEE Press.


[48] S. Straburger, G. Schmidgall, and S. Haasis. Distributed manufacturing simulation as an enabling technology for the digital factory. Journal of Advanced ManufacturingSystem, 2(1):111 -126, 2003.

42 Discrete Event Simulations – Development and Applications

International, San Diego, Calif. 1988.

for discrete event simulations, 1998.

2005.

USA, 2003.

81(1):5-13, 2005.

73(5):296-303, 1999.

simulation. SIMULATION, 85(2):113-128, 2009.

[34] Navonil Mustafee, Simon J.E. Taylor, Korina Katsaliaki, and Sally Brailsford. Facilitating the analysis of a uk national blood service supply chain using distributed

[35] Brian Unger, and David Jefferson. Distributed simulation, 1988 : proceedings of the SCS Multiconference on Distributed Simulation, 3-5 February, 1988, San Diego, California / edited by Brian Unger and David Jefierson. Society for Computer Simulation

[36] Ernest H. Page and Roger Smith. Introduction to military training simulation: A guide

[37] Jaebok Park, R. Moraga, L. Rabelo, J. Dawson, M.N. Marin, and J. Sepulveda. Addressing complexity using distributed simulation: a case study in spaceport modeling. In Simulation Conference, 2005 Proceedings of the Winter, page 9 pp., dec.

[38] G. Pedrielli, P. Scavardone, T. Tolio, M. Sacco, and W. Terkaj. Simulation of complex manufacturing systems via hla-based infrastructure. In Principles of Advanced and

[40] A. R. Pritchett, M. M. van Paassen, F. P. Wieland, and E. N. Johnson. Aerospace vehicle and air traffic simulation, pages 365-389. Kluwer Academic Publishers, Norwell, MA,

[41] M. Raab, S. Masik, and T. Schulze. Support system for distributed hla simulations in industrial applications. In Principles of Advanced and Distributed Simulation (PADS),

[42] M. Rabe, F.W. Jkel, and H. Weinaug. Supply chain demonstrator based on federated

[43] Stewart Robinson. Distributed simulation and simulation practice. SIMULATION,

[44] Marco Sacco, Giovanni Dal Maso, Ferdinando Milella, Paolo Pedrazzoli, Diego Rovere, and Walter Terkaj. Virtual Factory Manager, volume 6774 of Lecture Notes in

[45] B. Sadoun. Simulation in city planning and engineering, pages 315{341. Kluwer

[46] Thomas Schulze, Steffen Strassburger, and Ulrich Klein. Migration of hla into civil domains: Solutions and prototypes for transportation applications. SIMULATION,

[47] S. Straburger. Overview about the high level architecture for modeling and simulation

Computer Science, pages 397- 406. Springer Berlin, Heidelberg, 2011.

and recent developments. Simulation News Europe, 16(2):5-14, 2006.

Distributed Simulation (PADS), 2011 IEEE Workshop on, pages 1 -9, june 2011. [39] P. Peschlow and P. Martini. Eficient analysis of simultaneous events in distributed simulation. In Distributed Simulation and Real-Time Applications, 2007. DS-RT 2007.

11th IEEE International Symposium, pages 244-251, oct. 2007.

2011 IEEE Workshop on, pages 1 -7, june 2011.

Academic Publishers, Norwell, MA, USA, 2003.

models and hla application. 2006.


simulations and creating large models. In Simulation Conference (WSC), Proceedings of the 2011 Winter, pages 161-175, dec. 2011.

Distributed Modeling of Discrete Event Systems 45

[72] Y. Zhang and L. M. Zhang. The Viewable Distributed Simulation Linkage Development Tool Based on Factory Mechanism. Applied Mechanics and Materials, 58:1813-1818,

[73] Cosby, L.N. SIMNET: An Insider's Perspective. Ida documenT D-1661, Institute for Defense Analyses 1801 N. Beauregard Street, Alexandria, Virginia 22311-1772. pp: 1-19.

[75] SISO COTS Simulation Package Interoperability Product Development Group, (2010). Standard for Commercial - off - the - shelf Simulation Package Interoperability

[76] Kubat, C., Uygun O. (2007). HLA Based Distributed Simulation Model for Integrated Maintenance and Production Scheduling System in Textile Industry. *P.T. Pham, E.E. Eldukhri, A. Soroka (Eds.), Proceedings of 3rd I\*PROMS Virtual International Conference, 2–*

[77] Pedrielli, G., Sacco, M, Terkaj, W., Tolio, T. (2012). Simulation of complex manufacturing systems via HLA-based infrastructure. *Journal Of Simulation*. To be

[78] Kewley, R., Cook, J., Goerger, N., Henderson, D., Teague, E., (2008). "Federated simulations for systems of systems integration," *Simulation Conference, 2008. WSC 2008.* 

[79] Bruccoleri M, Capello C, Costa A, Nucci F, Terkaj W, Valente A (2009) Testing. In: Tolio T (ed) Design of Flexible Production Systems. Springer: 239-293. ISBN 978-3-540-85413-5.

[82] Ke Pan; Turner, S.J.; Wentong Cai; Zengxiang Li; , "Implementation of Data Distribution Management services in a Service Oriented HLA RTI," *Simulation Conference (WSC), Proceedings of the 2009 Winter* , vol., no., pp.1027-1038, 13-16 Dec. 2009

[83] Colledani M, Terkaj W, Tolio T, Tomasella M (2008) Development of a Conceptual Reference Framework to manage manufacturing knowledge related to Products, Processes and Production Systems. In Bernard A, Tichkiewitch S (eds) Methods and Tools for Effective Knowledge Life-Cycle-Management. Springer: 259-284. ISBN 978-3-

[84] Colledani M, Terkaj W, Tolio T (2009) Product-Process-System Information Formalization. In: Tolio T (ed) Design of Flexible Production Systems. Springer: 63-86.

[85] Pedrazzoli, P, Sacco, M, Jönsson, A, Boër, C (2007) Virtual Factory Framework: Key Enabler For Future Manufacturing. In Cunha, PF, Maropoulos, PG (eds) Digital

June 2011.

March 1995

[74] SISO CSPI-PDG www.sisostds.org.

*13 July 2007*: 413–418.

published.

[80] http://www.pitch.se/ [81] MAK RTI, www.mak.com

540-78430-2.

ISBN 978-3-540-85413-5.

doi: 10.1109/WSC.2009.5429557

Referenve Models. SISO-STD-006-2010.

*Winter* , vol., no., pp.1121-1129, 7-10 Dec. 2008.

Enterprise Technology, Springer US, pp 83-90


the 2011 Winter, pages 161-175, dec. 2011.

Transactions on, 36(1):109- 122, jan. 2006.

Transactions on, 36(1):109- 122, jan 2006.

Computers in Industry, 53(1):3-16, 2004.

2008.

2005.

2008.

1709-1713, aug. 2007.

simulations and creating large models. In Simulation Conference (WSC), Proceedings of

[61] S.J.E. Taylor, Xiaoguang Wang, S.J. Turner, and M.Y.H. Low. Integrating heterogeneous distributed cots discrete-event simulation packages: an emerging standards-based approach. Systems, Man and Cybernetics, Part A: Systems and Humans, IEEE

[62] S.J.E. Taylor, Xiaoguang Wang, S.J. Turner, and M.Y.H. Low. Integrating heterogeneous distributed cots discrete-event simulation packages: an emerging standards-based approach. Systems, Man and Cybernetics, Part A: Systems and Humans, IEEE

[63] Sergio Terzi and Sergio Cavalieri. Simulation in the supply chain context: a survey.

[64] J. Vancza, P. Egri, and L. Monostori. A coordination mechanism for rolling horizon planning in supply networks. CIRP Annals - Manufacturing Technology, 57(1):455- 458,

[65] Gabriel A. Wainer, Olivier Khaldoon, Al-Zoub Dalle, David R.C. Hill Hill, Saurabh Mittal, Jos L. Risco Martn, Hessam Sarjoughian, Luc Touraille, Mamadou K. Traor, and

[66] Gang Wang, Chun Jin, and Peng Gao. Adapting arena into hla: Approach and experiment. In Automation and Logistics, 2007 IEEE International Conference on, pages

[67] Xiaoguang Wang, Stephen John Turner, and Simon J. E. Taylor. Cots simulation package (csp) interoperability -a solution to synchronous entity passing. In Proceedings of the 20th Workshop on Principles of Advanced and Distributed Simulation, PADS '06,

[68] Xiaoguang Wang, Stephen John Turner, Simon J. E. Taylor, Malcolm Yoke Hean Low, and Boon Ping Gan. A cots simulation package emulator (cspe) for investigating cots simulation package interoperability. In Proceedings of the 37th conference on Winter

[69] Xiaoguang Wang, Sthephen John Turner, Malcolm Yoke Hean Low, and Boon Ping Gan. A generic architecture for the integration of cots packages with the hla. In Proceedings of the 2004 Operational Research Society Simulation Workshop, Special Interest Group for Simulation, pages 225-233. Association for Computing Machinery's,

[70] Gregory Zacharewicz, Claudia Frydman, and Norbert Giambiasi. G-devs/hla environment for distributed simulations of workows. Simulation, 84:197-213, May

[71] zer Uygun, Ercan ztemel, and Cemalettin Kubat. Scenario based distributed manufacturing simulation using hla technologies. Information Sciences, 179(10):1533- 1541, 2009. <ce:title>Including Special Issue on Artificial Imune Systems</ce:title>.

Zeigler Bernard P. Standardizing DEVS model representation. 2010.

pages 201-210, Washington, DC, USA, 2006. IEEE Computer Society.

simulation, WSC '05, pages 402-411. Winter Simulation Conference, 2005.

	- [86] Sacco M, Pedrazzoli P, Terkaj W (2010) VFF: Virtual Factory Framework. Proceedings of 16th International Conference on Concurrent Enterprising, Lugano, Switzerland, 21-23 June 2010.

**Chapter 0**

**Chapter 2**

**The Speedup of Discrete Event Simulations**

Discrete event simulation (DES) has been widely used for both theoretical and application researches in various fields. The low-cost numerical experiments by DES can be carried out repeatedly to investigate the collective behavior or dynamics of a complex system which usually consists of a huge number of elements or is too expensive to be established realistically [3]. As for engineering applications, new algorithms, measures, or even standards, often require a great deal of DES experiments before actual deploying, especially in the field of communication networks. However, the feasibility of DES method relies not only on the correctness of computational model of the target system but also the spend in the time of computation. An experiment that costs over days of computer computation will limit the

For a DES system designed for communication network simulation, which is implemented by the event-driven method, computation tasks waiting for processing are usually represented by their corresponding events. These events need to be maintained in-sequence by a systematic scheduler. The causality that a scheduler must preserve is usually guaranteed by sorting and dispatching according to the time of event. An event is to be suspended unless its time value is smaller than or equal to the present time. When a great number of tasks or events are in pending, the cost of event scheduling becomes the top contributor to the computational time. For a large scale communication network simulation, for example, scheduling can consume over 40% of the total computational time [12]. To accelerate scheduling, several practical algorithms have been proposed and adopted for some widely used simulators such as NS2 [9] and OMNet++ [14]. These algorithms can be categorized, according to the data structure used, into linear list, partitioning list, heap and tree. Among them, the partitioning list has

The partitioning list algorithm proposed by R. Brown is called Calendar Queue (CQ) [4], where a data structure, named bucket, is used to partition events into a sequence of sub-lists.

and reproduction in any medium, provided the original work is properly cited.

©2012 Wang and Yang, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly

© 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution,

**by Utilizing CPU Caching**

Additional information is available at the end of the chapter

Wennai Wang and Yi Yang

http://dx.doi.org/10.5772/50397

application area of DES considerably.

attracted much attention.

cited.

**1. Introduction**

[87] Kalpakjian Serope; Schmid Steven R. Book title: Tecnologia Meccanica, Editor: Pearon Education Year: 2008 The Figure 12 picture was inspired by the book

## **The Speedup of Discrete Event Simulations by Utilizing CPU Caching**

Wennai Wang and Yi Yang

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50397

### **1. Introduction**

46 Discrete Event Simulations – Development and Applications

June 2010.

[86] Sacco M, Pedrazzoli P, Terkaj W (2010) VFF: Virtual Factory Framework. Proceedings of 16th International Conference on Concurrent Enterprising, Lugano, Switzerland, 21-23

[87] Kalpakjian Serope; Schmid Steven R. Book title: Tecnologia Meccanica, Editor: Pearon

Education Year: 2008 The Figure 12 picture was inspired by the book

Discrete event simulation (DES) has been widely used for both theoretical and application researches in various fields. The low-cost numerical experiments by DES can be carried out repeatedly to investigate the collective behavior or dynamics of a complex system which usually consists of a huge number of elements or is too expensive to be established realistically [3]. As for engineering applications, new algorithms, measures, or even standards, often require a great deal of DES experiments before actual deploying, especially in the field of communication networks. However, the feasibility of DES method relies not only on the correctness of computational model of the target system but also the spend in the time of computation. An experiment that costs over days of computer computation will limit the application area of DES considerably.

For a DES system designed for communication network simulation, which is implemented by the event-driven method, computation tasks waiting for processing are usually represented by their corresponding events. These events need to be maintained in-sequence by a systematic scheduler. The causality that a scheduler must preserve is usually guaranteed by sorting and dispatching according to the time of event. An event is to be suspended unless its time value is smaller than or equal to the present time. When a great number of tasks or events are in pending, the cost of event scheduling becomes the top contributor to the computational time. For a large scale communication network simulation, for example, scheduling can consume over 40% of the total computational time [12]. To accelerate scheduling, several practical algorithms have been proposed and adopted for some widely used simulators such as NS2 [9] and OMNet++ [14]. These algorithms can be categorized, according to the data structure used, into linear list, partitioning list, heap and tree. Among them, the partitioning list has attracted much attention.

The partitioning list algorithm proposed by R. Brown is called Calendar Queue (CQ) [4], where a data structure, named bucket, is used to partition events into a sequence of sub-lists.

©2012 Wang and Yang, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 2 Will-be-set-by-IN-TECH 48 Discrete Event Simulations – Development and Applications The Speedup of Discrete Event Simulations by Utilizing CPU Caching <sup>3</sup>

CQ can reduce the time for event searching and decrease the complexity down to O(1). Since then, some improved algorithms, such as Dynamic CQ (DCQ) [1], SNOOPy CQ [13], and sluggish CQ [15], have been proposed in order to improve the adaptability of CQ to the event distribution in the time domain. However, the way of event enqueue and dequeue remains unchanged. For these algorithms, the complexity of O(1) is only valid for bucket searching and it doesn't take into account event locating within a bucket. In addition, the overhead for bucket managing inevitably results in some negative effects on the performance of computation. K. Chung, J. Sang and V. Rego compared some earlier-day scheduling algorithms and found that, for a token-ring network simulation, both CQ and DCQ are no better than those based on heap or tree structure [6].

**Scheduler**

**Application**

Aiming at network simulating at packet level, NS2 has been developed as a platform consisting of a comprehensive network elements by the object-oriented programming method. Four classes of fundamental objects, named Node, Link, Agent, and Application, are modeled for packets transporting, generating and terminating. While, functionalities such as configuration helping, routes computing, experiment managing and so on, are put into a

Unlike Agent and Application, both Node and Link are compound classes which consist of some other elemental objects. For example, a simplest Node contains two routing objects, one for addressing/routing at network layer and the other for multiplexing/demultiplexing at transport layer, both derived from the same parent class Classifier. The physical link is modeled by Link which consists in series of three objects: Queue, LinkDelay, and TTLChecker. The elemental classes including Agent, Queue, LinkDelay, and TTLChecker are all derived from the parent class Connector. Together with Classifier,

In the class NsObject, two abstract functions named recv() and send() are defined for simulating packet receiving and sending, respectively. The class Handler encapsulates a function named handle() for specific event handling. These functions are overridden in the derived classes for different purposes. For instance, the implementation of recv() in Classifier, ie. Classifier::recv(), is to forward the received packet to the next object predefined by Simulator according to routing logics. Details of the implementation

During network simulating, a packet handling causes usually one or more elemental objects to change their state. In most cases, the change results in at least one new event which relates to another function Handler::handle() of the next specified object. This event-driven invocation is managed systematically by a global sole object of Scheduler, which is instantiated and kept by another global object of Simulator. The derived classes based on

Connector is derived from NsObject and Handler in turn, as depicted in Fig.1.

**Process**

**NsObject**

**Handler TclObject**

**Queue TTLChecker Agent**

**Figure 1.** The hierarchical structure of main classes of NS2 in UML schema

**LinkDelay**

global container class and named after Simulator.

mechanism of NS2 can be referred to [10].

**2.1. Event handling mechanism**

**Classifier Connector**

**1**

The Speedup of Discrete Event Simulations by Utilizing CPU Caching 49

**CalendarScheduler**

**1**

**Simulator**

**ListScheduler**

In the past decades, researchers have put their interesting into parallel and distributed approaches [7] to speedup simulation by a cluster of computers or processors. Parallel and distributed event scheduling need to cope with challenges such as time synchronization and task balance and remain a gap to extensive adoption. An alternative way to speedup DES simulation, we think, is to exploit the potentials offered by the modern processors. As an instance of attempts, recently H. Park and P.A. Fishwick proposed a graphics processing units (GPU) based algorithm which shows 10-fold speedup [11]. However, until now the caching mechanism of processor has not been considered in existing event scheduling algorithms. Herein we provide a cache aware algorithm and verify its improvement on the performance by extending the conventional NS2.

The motivation to utilize caching of processor or CPU is straightforward. Both cache based method and cache aware method have been applied successfully in some heavy data-retrieving algorithms, including IP routing table lookup [5] and Radix tree implementation [2], where a huge amount of data are frequently accessed. This situation happens for a large scale DES as well, and CPU caching is benefit to speed up event scheduling.

The chapter is organized as follows. Section 2 gives a analysis of event driven mechanism of the NS2 simulator, an estimation on the number of event of the typical application, and a description of the conventional CQ algorithm. A cache aware algorithm is presented in Section 3, followed by a complexity analysis on enqueue and dequeue operations. In section 4, experiments to verify CPU cache awareness are provided, aiming at the evaluation of performance and its relationship with the size of event queue. Section 5 summarizes the chapter.

#### **2. NS2 and Calendar Queue algorithm**

For the typical event-driven DES system, simulation events waiting for processing are buffered in a queue and sorted according to their simulation time or timestamp. A simple and natural choice is to introduce a linear list to manage the queue of event. The number and distribution of events in the time domain hence dominate the performance of queue management. Such performance depends on not only simulation complexity, but also the implementation of a simulator. In order to make the issue tractable, the following analyses are developed on the widely used network simulator NS2 [9].

**Figure 1.** The hierarchical structure of main classes of NS2 in UML schema

#### **2.1. Event handling mechanism**

2 Will-be-set-by-IN-TECH

CQ can reduce the time for event searching and decrease the complexity down to O(1). Since then, some improved algorithms, such as Dynamic CQ (DCQ) [1], SNOOPy CQ [13], and sluggish CQ [15], have been proposed in order to improve the adaptability of CQ to the event distribution in the time domain. However, the way of event enqueue and dequeue remains unchanged. For these algorithms, the complexity of O(1) is only valid for bucket searching and it doesn't take into account event locating within a bucket. In addition, the overhead for bucket managing inevitably results in some negative effects on the performance of computation. K. Chung, J. Sang and V. Rego compared some earlier-day scheduling algorithms and found that, for a token-ring network simulation, both CQ and DCQ are no

In the past decades, researchers have put their interesting into parallel and distributed approaches [7] to speedup simulation by a cluster of computers or processors. Parallel and distributed event scheduling need to cope with challenges such as time synchronization and task balance and remain a gap to extensive adoption. An alternative way to speedup DES simulation, we think, is to exploit the potentials offered by the modern processors. As an instance of attempts, recently H. Park and P.A. Fishwick proposed a graphics processing units (GPU) based algorithm which shows 10-fold speedup [11]. However, until now the caching mechanism of processor has not been considered in existing event scheduling algorithms. Herein we provide a cache aware algorithm and verify its improvement on the performance

The motivation to utilize caching of processor or CPU is straightforward. Both cache based method and cache aware method have been applied successfully in some heavy data-retrieving algorithms, including IP routing table lookup [5] and Radix tree implementation [2], where a huge amount of data are frequently accessed. This situation happens for a large scale DES as well, and CPU caching is benefit to speed up event

The chapter is organized as follows. Section 2 gives a analysis of event driven mechanism of the NS2 simulator, an estimation on the number of event of the typical application, and a description of the conventional CQ algorithm. A cache aware algorithm is presented in Section 3, followed by a complexity analysis on enqueue and dequeue operations. In section 4, experiments to verify CPU cache awareness are provided, aiming at the evaluation of performance and its relationship with the size of event queue. Section 5 summarizes the

For the typical event-driven DES system, simulation events waiting for processing are buffered in a queue and sorted according to their simulation time or timestamp. A simple and natural choice is to introduce a linear list to manage the queue of event. The number and distribution of events in the time domain hence dominate the performance of queue management. Such performance depends on not only simulation complexity, but also the implementation of a simulator. In order to make the issue tractable, the following analyses are

better than those based on heap or tree structure [6].

by extending the conventional NS2.

**2. NS2 and Calendar Queue algorithm**

developed on the widely used network simulator NS2 [9].

scheduling.

chapter.

Aiming at network simulating at packet level, NS2 has been developed as a platform consisting of a comprehensive network elements by the object-oriented programming method. Four classes of fundamental objects, named Node, Link, Agent, and Application, are modeled for packets transporting, generating and terminating. While, functionalities such as configuration helping, routes computing, experiment managing and so on, are put into a global container class and named after Simulator.

Unlike Agent and Application, both Node and Link are compound classes which consist of some other elemental objects. For example, a simplest Node contains two routing objects, one for addressing/routing at network layer and the other for multiplexing/demultiplexing at transport layer, both derived from the same parent class Classifier. The physical link is modeled by Link which consists in series of three objects: Queue, LinkDelay, and TTLChecker. The elemental classes including Agent, Queue, LinkDelay, and TTLChecker are all derived from the parent class Connector. Together with Classifier, Connector is derived from NsObject and Handler in turn, as depicted in Fig.1.

In the class NsObject, two abstract functions named recv() and send() are defined for simulating packet receiving and sending, respectively. The class Handler encapsulates a function named handle() for specific event handling. These functions are overridden in the derived classes for different purposes. For instance, the implementation of recv() in Classifier, ie. Classifier::recv(), is to forward the received packet to the next object predefined by Simulator according to routing logics. Details of the implementation mechanism of NS2 can be referred to [10].

During network simulating, a packet handling causes usually one or more elemental objects to change their state. In most cases, the change results in at least one new event which relates to another function Handler::handle() of the next specified object. This event-driven invocation is managed systematically by a global sole object of Scheduler, which is instantiated and kept by another global object of Simulator. The derived classes based on

4 Will-be-set-by-IN-TECH 50 Discrete Event Simulations – Development and Applications The Speedup of Discrete Event Simulations by Utilizing CPU Caching <sup>5</sup>

transmitting. Therefore, the number *N* of pending events is as follows,

*N* =

network with 100 nodes and a full mesh typed traffic pattern, *n* equals 9900.

It is easy to figure out that the number of events is approximately,

can be seen taken roughly as the low bound for event scheduling.

speed of event scheduling becomes a critical factor for fast simulation.

are located in bucket *B*(0), *E*<sup>3</sup> in *B*(1), and *E*<sup>4</sup> in *B*(2) of the third year.

time, 16, of *E*<sup>4</sup> gives a value 2, indexing the bucket *B*(2).

**2.3. The Calendar Queue scheduling**

stored in the same bucket as shown in Fig.3.

ms. Then,

*n* ∑ *k*=1

where, *k* ∈ [1..*n*], is also the index of Applications of *n* concurrent traffic flows. For a

For a typical traffic with bandwidth demand 1 Mbps and packet size 1000 bytes, the packet generation rate *r* is 125 packet/s. Assuming the propagation delay along each link is fixed as 2 ms and each path consists of 6 hops in average, we have *hk* = 6 and OWD is greater than 12

Here, for simplicity, we neglect the effects of packet transmitting and queuing at an output port. As one knows, transmitting delay is inverse proportion to the bandwidth of output link, and it contributes an increment to OWD. Queuing can bring about packet drops and introduce a decrement to *gk* during network congestion. However, as it is the result of overloading, network congestion implies much more pending events. Therefore, the above estimated *N*

If *N* events are buffered in single one linear list, the average times of memory accessing for sorted insertion is about *N*/2 (105). Such heavy volume of accessing can lead to much too high time overhead. For networks with burst pattern traffics and higher traffic demands, the

The default event scheduler of NS2 is based on Calendar Queue (CQ) algorithm, which works similar to a desktop calendar [4]. In CQ, a bucket is used to stand for the day of year and store a sub-list of events. Events happening on the same day, despite the deference of year, are

For simplicity, 5 events, *E*<sup>0</sup> to *E*4, are illustrated in Fig.3, whose timestamps are defined as 0, 1, 1, 3 and 16 in second. Given that the depth or size of bucket, *TB*, is 2 seconds in time and one year consists of 3 buckets, the length of year, *Ty*, is 6 seconds. Therefore, events *E*<sup>0</sup> to *E*<sup>2</sup>

where, *te* is the timestamp of an event. *E*0, *E*<sup>1</sup> and *E*<sup>2</sup> are inserted in the bucket *B*(0) since computations according to Eq.(4) result in the same value 0. Again, the computation of the

The computation of Eq.(4) is independent on the number of events in the queue, the complexity is therefore of O(1) [4]. However, event enqueue in a single bucket works as the

The bucket location or index is determined by the following arithmetic computation,

*gk* = *r* × *OWD* = 125(*packet*/*s*) × 12(*ms*) = 1.5. (2)

*<sup>N</sup>* <sup>=</sup> <sup>9900</sup> <sup>×</sup> 1.5 <sup>×</sup> <sup>13</sup> <sup>≈</sup> <sup>2</sup> <sup>×</sup> 105. (3)

*nB* = �(*te* mod *Ty*)/*TB*�, (4)

*gk* × (2*hk* + 1), (1)

The Speedup of Discrete Event Simulations by Utilizing CPU Caching 51

**Figure 2.** The configuration node and link of NS2 and typical events written in italic

the Scheduler include ListScheduler, which implements the linear list algorithm, and CalendarScheduler, which implements CQ and DCQ.

The event generation and consumption can be illustrated by tracing a packet transportation along a simulation link. As it can be imaged, the end of packet transferring at a output port of node will trigger a next serving event which will be sent back to the output queue, i.e. Queue, in the future time of simulation. While time delay due to propagation through the LinkDelay will create a packet arriving event which will be forwarded the input interface, i.e. TTLChecker, of the next node later. These two kinds of events lead to invocations of the class Queue and TTLChecker are buffered into the event queue of Scheduler, as showed in Fig.2.

Packet generating is controlled by the object Application, which has been developed and divided into two types which represent the stochastic generator and application simulator, respectively. The former creates an event that will trigger packet generating of Application recursively. The latter generates packet according to the mechanism of the simulated application. The following discussion will focus on the former for simplicity.

#### **2.2. Population of standing events**

As mentioned above, there are 3 types of events in NS2, including packet generating, timer and "at-event". Packet is used to model packet's transmission, timer to build protocol state-machine and packet's generator, and "at-event"to control computation. Application and its derivations are designed to generate a flow of packets, where one or more timer(s) is/are used to trigger the generation in a recursive manner. Queue and LinkDelay are designed for packet forwarding and transmitting through a link connecting a pair of neighbor nodes, respectively. Two timers are generated when a LinkDelay handles a received packet, one for triggering packet receiving at the next node later on, the other for calling Queue back when the LinkDelay is available for transmitting next packet buffered in the Queue. Hopping of a packet along a link, hence, induce two consequential events.

Let *hk* denotes the number of hops or links over an end-to-end path indexed by *k*. The number, represented by *gk*, of packets are generated by the corresponding Application during One Way Delay (OWD) of the path. The total number of Application's instances, *n*, is assumed as the same as the number of paths. Then, in the context of NS2, one packet-typed event induces 2*hk* timer-typed events if the packet does not drop during forwarding and transmitting. Therefore, the number *N* of pending events is as follows,

4 Will-be-set-by-IN-TECH

the Scheduler include ListScheduler, which implements the linear list algorithm, and

The event generation and consumption can be illustrated by tracing a packet transportation along a simulation link. As it can be imaged, the end of packet transferring at a output port of node will trigger a next serving event which will be sent back to the output queue, i.e. Queue, in the future time of simulation. While time delay due to propagation through the LinkDelay will create a packet arriving event which will be forwarded the input interface, i.e. TTLChecker, of the next node later. These two kinds of events lead to invocations of the class Queue and TTLChecker are buffered into the event queue of Scheduler, as showed

Packet generating is controlled by the object Application, which has been developed and divided into two types which represent the stochastic generator and application simulator, respectively. The former creates an event that will trigger packet generating of Application recursively. The latter generates packet according to the mechanism of the simulated

As mentioned above, there are 3 types of events in NS2, including packet generating, timer and "at-event". Packet is used to model packet's transmission, timer to build protocol state-machine and packet's generator, and "at-event"to control computation. Application and its derivations are designed to generate a flow of packets, where one or more timer(s) is/are used to trigger the generation in a recursive manner. Queue and LinkDelay are designed for packet forwarding and transmitting through a link connecting a pair of neighbor nodes, respectively. Two timers are generated when a LinkDelay handles a received packet, one for triggering packet receiving at the next node later on, the other for calling Queue back when the LinkDelay is available for transmitting next packet buffered in the Queue.

Let *hk* denotes the number of hops or links over an end-to-end path indexed by *k*. The number, represented by *gk*, of packets are generated by the corresponding Application during One Way Delay (OWD) of the path. The total number of Application's instances, *n*, is assumed as the same as the number of paths. Then, in the context of NS2, one packet-typed event induces 2*hk* timer-typed events if the packet does not drop during forwarding and

**Figure 2.** The configuration node and link of NS2 and typical events written in italic

application. The following discussion will focus on the former for simplicity.

Hopping of a packet along a link, hence, induce two consequential events.

CalendarScheduler, which implements CQ and DCQ.

**2.2. Population of standing events**

in Fig.2.

$$N = \sum\_{k=1}^{n} g\_k \times (2h\_k + 1)\_\prime \tag{1}$$

where, *k* ∈ [1..*n*], is also the index of Applications of *n* concurrent traffic flows. For a network with 100 nodes and a full mesh typed traffic pattern, *n* equals 9900.

For a typical traffic with bandwidth demand 1 Mbps and packet size 1000 bytes, the packet generation rate *r* is 125 packet/s. Assuming the propagation delay along each link is fixed as 2 ms and each path consists of 6 hops in average, we have *hk* = 6 and OWD is greater than 12 ms. Then,

$$g\_k = r \times \text{OWD} = 125(packet/s) \times 12(ms) = 1.5.\tag{2}$$

It is easy to figure out that the number of events is approximately,

$$N = 9900 \times 1.5 \times 13 \approx 2 \times 10^5. \tag{3}$$

Here, for simplicity, we neglect the effects of packet transmitting and queuing at an output port. As one knows, transmitting delay is inverse proportion to the bandwidth of output link, and it contributes an increment to OWD. Queuing can bring about packet drops and introduce a decrement to *gk* during network congestion. However, as it is the result of overloading, network congestion implies much more pending events. Therefore, the above estimated *N* can be seen taken roughly as the low bound for event scheduling.

If *N* events are buffered in single one linear list, the average times of memory accessing for sorted insertion is about *N*/2 (105). Such heavy volume of accessing can lead to much too high time overhead. For networks with burst pattern traffics and higher traffic demands, the speed of event scheduling becomes a critical factor for fast simulation.

#### **2.3. The Calendar Queue scheduling**

The default event scheduler of NS2 is based on Calendar Queue (CQ) algorithm, which works similar to a desktop calendar [4]. In CQ, a bucket is used to stand for the day of year and store a sub-list of events. Events happening on the same day, despite the deference of year, are stored in the same bucket as shown in Fig.3.

For simplicity, 5 events, *E*<sup>0</sup> to *E*4, are illustrated in Fig.3, whose timestamps are defined as 0, 1, 1, 3 and 16 in second. Given that the depth or size of bucket, *TB*, is 2 seconds in time and one year consists of 3 buckets, the length of year, *Ty*, is 6 seconds. Therefore, events *E*<sup>0</sup> to *E*<sup>2</sup> are located in bucket *B*(0), *E*<sup>3</sup> in *B*(1), and *E*<sup>4</sup> in *B*(2) of the third year.

The bucket location or index is determined by the following arithmetic computation,

$$m\_B = \lfloor (t\_\ell \bmod T\_\circ) / T\_B \rfloor,\tag{4}$$

where, *te* is the timestamp of an event. *E*0, *E*<sup>1</sup> and *E*<sup>2</sup> are inserted in the bucket *B*(0) since computations according to Eq.(4) result in the same value 0. Again, the computation of the time, 16, of *E*<sup>4</sup> gives a value 2, indexing the bucket *B*(2).

The computation of Eq.(4) is independent on the number of events in the queue, the complexity is therefore of O(1) [4]. However, event enqueue in a single bucket works as the

**Figure 4.** Event digests and their relationship with DES events

adaptive to event distribution in the time domain.

**3.2. Enqueue and dequeue operations**

including,

sub-list (*Ei*, ..., *Ek*−1). The time delimiter of digest is defined as follows,

as depicted in Fig.4.

ring-typed buffer. The element of *Qd*, called digest, is used to index the tail of a sub-list of *Qe*,

A digest *Dj* contains a time delimiter *tm* and an event pointer *ptr* to the corresponding DES event of a sub-list. In Fig.2, *D*<sup>0</sup> is the digest of the sub-list (*E*0, *E*1), and *Dj* is the digest of the

where the coefficient *g* is used to control the size of sub-list, and *ts* is the timestamp of the corresponding event. If *g* = 1, in unit of timestamp, a sub-list consists of only those events with the same timestamp. For *g* > 1, events happen within *g* interval are indexed by a digest pointing to the latest of them. In effect, the configurable coefficient *g* is analog to bucket size of CQ. The difference is that two successive digests need not have their *tm* contiguous. In other words, *Dj*<sup>+</sup><sup>1</sup> → *tm* - *Dj* → *tm* is allowed to be greater than 1×*g*. This makes our sub-list more

The value of *ptr* of a digest is the address of corresponding DES event. The function of this address is two-fold, one for heading of the next sub-list, the other for tailing of the current sub-list. The reason why the *ptr* points to the tail of a sublist is to avoid repeating comparisons for concurrent events with the same time-stamp. This situation happens more frequently in case simulation parameters, ex. bandwidth or link delay, are configured with the same value.

As described in Algorithm 1, inserting a new DES event, or ENQUEUE, is a little more complex operation than a linear list. There are 6 significant steps in ENQUEUE operation,

• to calculate the time maker, *x*, of the new DES event, according to Eq.(4);

• to create a digest and insert it into digest ring buffer by InsertDigQueue().

• to find the target sub-list by comparing *x* with each digest of *Qe*; • to sort and insert the new event into the sub-list by InsertList(); • to replace the digest if the new event is inserted at the tail, otherwise;

• to update the digest if digest queue is full, otherwise;

*tm* = �*g* × *ts*� (5)

The Speedup of Discrete Event Simulations by Utilizing CPU Caching 53

**Figure 3.** The structural relationship between event and bucket in CQ

same way as a linear list. In the case as showed in Fig.3, the insertion of *E*<sup>2</sup> will require two extra steps of comparing. The final complexity is dominated by the sorted insertion in the bucket unless the number of bucket is greater than the total number of events and the distribution of event in time is uniform.

Actually, bucket plays as a container and its size determines events partitioning. The fixed structure of bucket becomes an obstacle for performance improvement. It has been shown that a skewed nonuniform distribution of events in the time domain can degrade the performance of CQ [12]. Since then, several improved scheduling algorithms [1, 13, 15] have been proposed to solve the problem of adaptability to the distribution.

In addition to the structure of bucket, there are two issues that CQ and the related improved algorithms can not cope with readily. One is concurrent events insertion and the other is the earliest event fetching. The former can be seen from the insertion of *E*<sup>2</sup> after *E*<sup>1</sup> as shown in Fig.3. Although they have equal timestamp, logically *E*<sup>2</sup> occurs later than *E*<sup>1</sup> and one more comparison is needed to insert *E*2. If the number of concurrent events is huge, say 9900 as cited in Eq.(3), a large number of comparisons are needed. The latter can be illustrated via *E*4's fetching. After *E*<sup>3</sup> departing, a CQ scheduler will carry out 6 times of bucket access across two years, then reach *E*4. This is more complex than a linear list, in which the earliest event is always located at the head.

As is seen, CQ and the related improved algorithms do not take the modern structure of processor into considerations, especially high speed CPU caches. It has been demonstrated that a cache aware algorithm can speedup greatly an application that involves heavy data-retrieving [2, 5]. In the following section, we provide a fast cache aware scheduling algorithm to accelerate the simulation of large scale communication networks.

### **3. The cache aware scheduling**

#### **3.1. Data structure for event partitioning**

Similar to CQ, the algorithm with cache awareness belongs to the category of partitioning list. Two lists work in a correlated manner, one for DES events and named event queue (*Qe*), and the other for indexing the sub-list of DES events and named digest queue (*Qd*). *Qe* is organized the same as a linear list, while *Qd* is implemented in an array structure acts as a

**Figure 4.** Event digests and their relationship with DES events

6 Will-be-set-by-IN-TECH

same way as a linear list. In the case as showed in Fig.3, the insertion of *E*<sup>2</sup> will require two extra steps of comparing. The final complexity is dominated by the sorted insertion in the bucket unless the number of bucket is greater than the total number of events and the

Actually, bucket plays as a container and its size determines events partitioning. The fixed structure of bucket becomes an obstacle for performance improvement. It has been shown that a skewed nonuniform distribution of events in the time domain can degrade the performance of CQ [12]. Since then, several improved scheduling algorithms [1, 13, 15] have been proposed

In addition to the structure of bucket, there are two issues that CQ and the related improved algorithms can not cope with readily. One is concurrent events insertion and the other is the earliest event fetching. The former can be seen from the insertion of *E*<sup>2</sup> after *E*<sup>1</sup> as shown in Fig.3. Although they have equal timestamp, logically *E*<sup>2</sup> occurs later than *E*<sup>1</sup> and one more comparison is needed to insert *E*2. If the number of concurrent events is huge, say 9900 as cited in Eq.(3), a large number of comparisons are needed. The latter can be illustrated via *E*4's fetching. After *E*<sup>3</sup> departing, a CQ scheduler will carry out 6 times of bucket access across two years, then reach *E*4. This is more complex than a linear list, in which the earliest

As is seen, CQ and the related improved algorithms do not take the modern structure of processor into considerations, especially high speed CPU caches. It has been demonstrated that a cache aware algorithm can speedup greatly an application that involves heavy data-retrieving [2, 5]. In the following section, we provide a fast cache aware scheduling

Similar to CQ, the algorithm with cache awareness belongs to the category of partitioning list. Two lists work in a correlated manner, one for DES events and named event queue (*Qe*), and the other for indexing the sub-list of DES events and named digest queue (*Qd*). *Qe* is organized the same as a linear list, while *Qd* is implemented in an array structure acts as a

algorithm to accelerate the simulation of large scale communication networks.

**Figure 3.** The structural relationship between event and bucket in CQ

to solve the problem of adaptability to the distribution.

distribution of event in time is uniform.

event is always located at the head.

**3. The cache aware scheduling**

**3.1. Data structure for event partitioning**

ring-typed buffer. The element of *Qd*, called digest, is used to index the tail of a sub-list of *Qe*, as depicted in Fig.4.

A digest *Dj* contains a time delimiter *tm* and an event pointer *ptr* to the corresponding DES event of a sub-list. In Fig.2, *D*<sup>0</sup> is the digest of the sub-list (*E*0, *E*1), and *Dj* is the digest of the sub-list (*Ei*, ..., *Ek*−1). The time delimiter of digest is defined as follows,

$$tm = \lfloor \lg \times t \text{s} \rfloor \tag{5}$$

where the coefficient *g* is used to control the size of sub-list, and *ts* is the timestamp of the corresponding event. If *g* = 1, in unit of timestamp, a sub-list consists of only those events with the same timestamp. For *g* > 1, events happen within *g* interval are indexed by a digest pointing to the latest of them. In effect, the configurable coefficient *g* is analog to bucket size of CQ. The difference is that two successive digests need not have their *tm* contiguous. In other words, *Dj*<sup>+</sup><sup>1</sup> → *tm* - *Dj* → *tm* is allowed to be greater than 1×*g*. This makes our sub-list more adaptive to event distribution in the time domain.

The value of *ptr* of a digest is the address of corresponding DES event. The function of this address is two-fold, one for heading of the next sub-list, the other for tailing of the current sub-list. The reason why the *ptr* points to the tail of a sublist is to avoid repeating comparisons for concurrent events with the same time-stamp. This situation happens more frequently in case simulation parameters, ex. bandwidth or link delay, are configured with the same value.

#### **3.2. Enqueue and dequeue operations**

As described in Algorithm 1, inserting a new DES event, or ENQUEUE, is a little more complex operation than a linear list. There are 6 significant steps in ENQUEUE operation, including,


#### **Algorithm 1:** ENQUEUE operation over a non-empty digest queue

**input** : A new event *new*, a DES event list headed by *H*, and a digest ring headed by *h* and tailed by *t*. **output**: Updated event list and digest ring. *x* ←− �*g* × *new* → *ts*�; *list* ←− *H*; **for** *j* ←− *h* **to** *t* **do** /\*To find the target sub-list\*/ **if** *D*[*j*] → *tm* > *x* **then break; end end if** *j != h* **then** *j* ←− *j* − 1; *list* ←− *D*[*j*] → *ptr*; **end** InsertList (*list*, *new*); **if** *D*[*j*] → *tm* == *x* **then** /\*To check the time maker of the digest\*/ **if** *list* → *ts* <= *new* → *ts* **then** /\*To replace the end of sub-list\*/ (*D*[*j*] → *ptr*) ←− *new*; **return; end else if** *t* + 1 ==*h* **then** /\*To update the digest\*/ (*D*[*j*] → *ptr*) ←− *new*; (*D*[*j*] → *tm*) ←− *x*; **else** /\*To create a digest and insert it at (*j*+1)-th of the digest ring \*/ InsertDigQueue (*new*, *x*, *j* + 1); **end**

**3.3. Cache awareness and computational complexity**

inserting (InsertDigQueue()) execute at first over the array.

of sub-list searching is equivalent to O( *<sup>m</sup>*

conclusion, the algorithm is better than CQ.

**4.1. Description of implementation**

<sup>2</sup>*<sup>β</sup>* ) + O( *<sup>n</sup>*

Consequently, the lowest complexity of ENQUEUE is O(*<sup>n</sup>*

**4. Implementation and performance evaluations**

size of each sub-list bing *<sup>n</sup>*

the complexity of O( *<sup>m</sup>*

follows,

O(*<sup>n</sup> β* ).

latter is cache awareness and has a better compatibility for general purposes.

The difference between a cache aware algorithm and CQ is the way to partition events into sub-lists. For the cache aware algorithm, an array of digest is designed to index sub-lists and works as a ring buffer. In the ENQUEUE operation, both sub-list searching and digest

Generally speaking, a fixed array of data is with great possibility allocated by the memory management of an OS to be continuously distributed in the memory space. This locality is feasible for fast accessing by a processor with caches. Caches are fast memories in which a block of data around the request are loaded from the main memory. The caching operation is benefit to the future requests if they locate in the loaded block. The ratio, *β*, of access latency of main memory to cache, is usually greater than 10. Therefore, an optimum data structure and its accessing can be designed to utilize the benefit. Two types of approach have been used in designing, one is to access cache directly by instructions based on the hardware structure, the other is to design a deliberate data structure and access them as locally as possible. The former promises better performance but less applicable for different hardware structures, while the

Buckets of CQ/DCQ embedded in NS2 are allocated by an array and can assure locality. But this locality can not make use of the benefit of caching, because the bucket determination for event searching is computed directly and has no relationship with buckets' allocation.

For the cache aware algorithm, the target sub-list is searched via the digest queue *Qd*. If whole of *Qd* can be loaded into caches, the searching can speed up *β* times. Therefore, the complexity

with *m* = *βn* buckets. However, for the case of a sub-list consists of concurrent events with the same timestamp, the cache aware algorithm has an advantage over CQ since that the tail of sub-list is indexed. Sublist searching of the algorithm is of O(1) while that of CQ goes up of

As for DEQUEUE operation, the cache aware algorithm is equivalent to remove the head of a list. The complexity is O(1) and it is independent on the distribution of event in time. As a

The cache aware algorithm has been implemented within an extended class CacheScheduler which is derived from the existing ListScheduler. Two functions, insert() and deque(), are overridden, and some helper functions appended. Modifications are based on the NS with version 2.33 and should be compatible with

*<sup>m</sup>* , the complexity of InsertList() is then O( *<sup>n</sup>*

<sup>2</sup>*<sup>β</sup>* ). For an average distribution of event in time, the

The Speedup of Discrete Event Simulations by Utilizing CPU Caching 55

*m* = *βn*. (6)

*<sup>β</sup>* ), twice as that of CQ configured

<sup>2</sup>*<sup>m</sup>* ) for an ENQUEUE operation. The optimum condition is as

<sup>2</sup>*<sup>m</sup>* ). This leads to

Since an array structure is used for the digest queue, the conventional ring-buffer access method can be adopted to implement InsertDigQueue(). The function InsertList() is, however, the same as that for a linear list.

The event departure or DEQUEUE operation is simply as a linear list, as shown in Algorithm 2. The function FetchList() is to fetch the head of DES event list, and the following block of codes to update the variable *h* is to remove the first of digest queue. The extra processing on the head is used to keep the correctness of digesting.

#### **Algorithm 2:** DEQUEUE operation when digest queue is not empty

```
input : An event queue Qe, and a digest ring headed by h.
output: The earliest event evt.
evt ←− FetchList(Qe);
if D[h] → ptr == evt then
   h ←− h + 1;
end
```
#### **3.3. Cache awareness and computational complexity**

8 Will-be-set-by-IN-TECH

**input** : A new event *new*, a DES event list headed by *H*, and a digest ring headed by *h*

**Algorithm 1:** ENQUEUE operation over a non-empty digest queue

**if** *D*[*j*] → *tm* == *x* **then** /\*To check the time maker of the digest\*/ **if** *list* → *ts* <= *new* → *ts* **then** /\*To replace the end of sub-list\*/

**else** /\*To create a digest and insert it at (*j*+1)-th of the digest ring \*/

Since an array structure is used for the digest queue, the conventional ring-buffer access method can be adopted to implement InsertDigQueue(). The function InsertList()

The event departure or DEQUEUE operation is simply as a linear list, as shown in Algorithm 2. The function FetchList() is to fetch the head of DES event list, and the following block of codes to update the variable *h* is to remove the first of digest queue. The extra processing

and tailed by *t*.

**if** *D*[*j*] → *tm* > *x* **then break;**

*list* ←− *D*[*j*] → *ptr*;

(*D*[*j*] → *ptr*) ←− *new*;

InsertDigQueue (*new*, *x*, *j* + 1);

is, however, the same as that for a linear list.

on the head is used to keep the correctness of digesting.

**Algorithm 2:** DEQUEUE operation when digest queue is not empty **input** : An event queue *Qe*, and a digest ring headed by *h*.

**else if** *t* + 1 ==*h* **then** /\*To update the digest\*/

InsertList (*list*, *new*);

**return;**

(*D*[*j*] → *ptr*) ←− *new*; (*D*[*j*] → *tm*) ←− *x*;

**output**: The earliest event *evt*. *evt* ←− FetchList(*Qe*); **if** *D*[*h*] → *ptr* == *evt* **then** *h* ←− *h* + 1;

*x* ←− �*g* × *new* → *ts*�;

*list* ←− *H*;

**end**

**if** *j != h* **then** *j* ←− *j* − 1;

**end**

**end**

**end**

**end**

**end**

**output**: Updated event list and digest ring.

**for** *j* ←− *h* **to** *t* **do** /\*To find the target sub-list\*/

The difference between a cache aware algorithm and CQ is the way to partition events into sub-lists. For the cache aware algorithm, an array of digest is designed to index sub-lists and works as a ring buffer. In the ENQUEUE operation, both sub-list searching and digest inserting (InsertDigQueue()) execute at first over the array.

Generally speaking, a fixed array of data is with great possibility allocated by the memory management of an OS to be continuously distributed in the memory space. This locality is feasible for fast accessing by a processor with caches. Caches are fast memories in which a block of data around the request are loaded from the main memory. The caching operation is benefit to the future requests if they locate in the loaded block. The ratio, *β*, of access latency of main memory to cache, is usually greater than 10. Therefore, an optimum data structure and its accessing can be designed to utilize the benefit. Two types of approach have been used in designing, one is to access cache directly by instructions based on the hardware structure, the other is to design a deliberate data structure and access them as locally as possible. The former promises better performance but less applicable for different hardware structures, while the latter is cache awareness and has a better compatibility for general purposes.

Buckets of CQ/DCQ embedded in NS2 are allocated by an array and can assure locality. But this locality can not make use of the benefit of caching, because the bucket determination for event searching is computed directly and has no relationship with buckets' allocation.

For the cache aware algorithm, the target sub-list is searched via the digest queue *Qd*. If whole of *Qd* can be loaded into caches, the searching can speed up *β* times. Therefore, the complexity of sub-list searching is equivalent to O( *<sup>m</sup>* <sup>2</sup>*<sup>β</sup>* ). For an average distribution of event in time, the size of each sub-list bing *<sup>n</sup> <sup>m</sup>* , the complexity of InsertList() is then O( *<sup>n</sup>* <sup>2</sup>*<sup>m</sup>* ). This leads to the complexity of O( *<sup>m</sup>* <sup>2</sup>*<sup>β</sup>* ) + O( *<sup>n</sup>* <sup>2</sup>*<sup>m</sup>* ) for an ENQUEUE operation. The optimum condition is as follows,

$$m = \sqrt{\beta n}.\tag{6}$$

Consequently, the lowest complexity of ENQUEUE is O(*<sup>n</sup> <sup>β</sup>* ), twice as that of CQ configured with *m* = *βn* buckets. However, for the case of a sub-list consists of concurrent events with the same timestamp, the cache aware algorithm has an advantage over CQ since that the tail of sub-list is indexed. Sublist searching of the algorithm is of O(1) while that of CQ goes up of O(*<sup>n</sup> β* ).

As for DEQUEUE operation, the cache aware algorithm is equivalent to remove the head of a list. The complexity is O(1) and it is independent on the distribution of event in time. As a conclusion, the algorithm is better than CQ.

#### **4. Implementation and performance evaluations**

#### **4.1. Description of implementation**

The cache aware algorithm has been implemented within an extended class CacheScheduler which is derived from the existing ListScheduler. Two functions, insert() and deque(), are overridden, and some helper functions appended. Modifications are based on the NS with version 2.33 and should be compatible with the most of other versions because no change is required except a config modification for Simulator is required to select CacheScheduler rather than CalendarScheduler in default.

which the new should be inserted, and the index idx of the sublist in the digest queue. The

The Speedup of Discrete Event Simulations by Utilizing CPU Caching 57

Codes in the above from Line 05 to 06 are to find out the sublist that the event timestamped with *key* should be contained. Codes from Line 08 to 10 are to handle the special case of the first sublist, and codes from Line 13 to 17 are introduced to compensate tail pointing mechanism used in sublist digesting, considering that the head of a sublist is equivalent to

The insertion function insertDigest() is little more complex as listed in the following,

05 if (tmp == head\_) { // buffer is full

07 } else if (head\_ == tail\_) { // buffer is empty

09 } else if (key\_[idx] == key) { // replace only

12 } else if (!e->next\_) { // append at tail

01 void CacheScheduler::insertDigest(double t, unsigned int key, \\

03 unsigned int tmp = (tail\_+1) % size\_; // target of tail moving

01 unsigned int CacheScheduler::findDigest(double t, \\ 02 unsigned int key, Event \*\*& p) {

05 for (; idx != tail\_; idx++, idx %= size\_)

13 if (idx == tail\_ || event\_[idx]->time\_ > t) { 14 // go to the tail of previous sub-list

c++ codes of findDigest() is defined as,

03 unsigned int idx = head\_;

15 idx = idx + size\_ - 1;

19 p = (Event\*\*)event\_[idx];

02 unsigned int idx, Event \*e) {

10 if (event\_[idx]->time\_ <= t)

11 event\_[idx] = e;

13 key\_[tail\_] = key; 14 event\_[tail\_] = e;

09 p = &queue\_; 10 return idx;

16 idx %= size\_;

20 return idx;

the tail of its previous one.

06 ...

08 ...

06 if (key\_[idx] >= key) break;

08 if (idx == head\_) { // is it at head

04

07

12

18

04

11 }

17 }

21 }

Two fixed array are defined in the class CacheScheduler, one named key\_ with type of unsigned int for digest's time, the other named event\_ with type of Event \* for digest's pointer. Two variable members, head\_ and tail\_, are defined to index the ring-typed buffer. The coefficient *g* in Eq.(5) is represented by the third variable member pricision\_, and the size of digest by the last size\_.

The overridden function deque() is defined in c++ as following,

```
01 Event* CacheScheduler::deque() {
02 Event *e = queue_;
03 if (e)} {
04 queue_ = e->next_;
05 if (event_[head_] == e) {
06 event_[head_] = 0;
07 key_[head_] = 0;
08 head_ = ++head_ % size_;
09 }
10 }
11 return (e);
12 }
```
where, the variable member queue\_ points to the head of DES event list. The condition (event\_[head\_] == e) is for checking whether it should remove the head of digest queue or not. The c++ codes of the function insert() is showed as following,

```
01 Event* CacheScheduler::insert(Event *e) {
02 Event **p;
03 unsigned int idx, key;
04 double t = e->time_;
05
06 key = t * pricision_;
07 idx = findDigest(t, key, p);
08
09 for( ;* p != 0; p = &(*p)->next_ )
10 if ( t < (*p)->time_ ) break;
11
12 e->next_ = *p;
13 *p = e;
14
15 insertDigest(t, key, idx, *p);
16 }
```
where, the functions findDigest() and insertDigest() operate over the array key\_ and event\_. The function findDigest() results the pointer p of sublist of DES event after which the new should be inserted, and the index idx of the sublist in the digest queue. The c++ codes of findDigest() is defined as,

10 Will-be-set-by-IN-TECH

the most of other versions because no change is required except a config modification for Simulator is required to select CacheScheduler rather than CalendarScheduler in

Two fixed array are defined in the class CacheScheduler, one named key\_ with type of unsigned int for digest's time, the other named event\_ with type of Event \* for digest's pointer. Two variable members, head\_ and tail\_, are defined to index the ring-typed buffer. The coefficient *g* in Eq.(5) is represented by the third variable member pricision\_, and the

where, the variable member queue\_ points to the head of DES event list. The condition (event\_[head\_] == e) is for checking whether it should remove the head of digest queue

where, the functions findDigest() and insertDigest() operate over the array key\_ and event\_. The function findDigest() results the pointer p of sublist of DES event after

or not. The c++ codes of the function insert() is showed as following,

01 Event\* CacheScheduler::insert(Event \*e) {

The overridden function deque() is defined in c++ as following,

01 Event\* CacheScheduler::deque() {

05 if (event\_[head\_] == e) { 06 event\_[head\_] = 0; 07 key\_[head\_] = 0;

08 head\_ = ++head\_ % size\_;

default.

size of digest by the last size\_.

03 if (e)} {

09 } 10 }

12 }

05

08

11

14

16 }

11 return (e);

02 Event \*\*p;

03 unsigned int idx, key; 04 double t = e->time\_;

06 key = t \* pricision\_;

12 e->next\_ = \*p;

13 \*p = e;

07 idx = findDigest(t, key, p);

15 insertDigest(t, key, idx, \*p);

09 for( ;\* p != 0; p = &(\*p)->next\_ ) 10 if ( t < (\*p)->time\_ ) break;

02 Event \*e = queue\_;

04 queue\_ = e->next\_;

```
01 unsigned int CacheScheduler::findDigest(double t, \\
02 unsigned int key, Event **& p) {
03 unsigned int idx = head_;
04
05 for (; idx != tail_; idx++, idx %= size_)
06 if (key_[idx] >= key) break;
07
08 if (idx == head_) { // is it at head
09 p = &queue_;
10 return idx;
11 }
12
13 if (idx == tail_ || event_[idx]->time_ > t) {
14 // go to the tail of previous sub-list
15 idx = idx + size_ - 1;
16 idx %= size_;
17 }
18
19 p = (Event**)event_[idx];
20 return idx;
21 }
```
Codes in the above from Line 05 to 06 are to find out the sublist that the event timestamped with *key* should be contained. Codes from Line 08 to 10 are to handle the special case of the first sublist, and codes from Line 13 to 17 are introduced to compensate tail pointing mechanism used in sublist digesting, considering that the head of a sublist is equivalent to the tail of its previous one.

The insertion function insertDigest() is little more complex as listed in the following,

```
01 void CacheScheduler::insertDigest(double t, unsigned int key, \\
02 unsigned int idx, Event *e) {
03 unsigned int tmp = (tail_+1) % size_; // target of tail moving
04
05 if (tmp == head_) { // buffer is full
06 ...
07 } else if (head_ == tail_) { // buffer is empty
08 ...
09 } else if (key_[idx] == key) { // replace only
10 if (event_[idx]->time_ <= t)
11 event_[idx] = e;
12 } else if (!e->next_) { // append at tail
13 key_[tail_] = key;
14 event_[tail_] = e;
```
12 Will-be-set-by-IN-TECH 58 Discrete Event Simulations – Development and Applications The Speedup of Discrete Event Simulations by Utilizing CPU Caching <sup>13</sup>

```
15 tail_ = tmp;
16 } else { // in the middle
17 idx = ++idx % size_;
18 if (key_[idx] > key) { // insert a new element
19 if (tail_ > idx) { // right moving
20 for (tmp = tail_; tmp > idx; tmp--) {
21 key_[tmp] = key_[tmp-1];
22 event_[tmp] = event_[tmp-1];
23 }
24 tail_ = ++tail_ % size_;
25 } else { // left moving
26 head_ --;
27 idx --;
28 for (tmp = head_; tmp < idx; tmp++) {
29 key_[tmp] = key_[tmp+1];
30 event_[tmp] = event_[tmp+1];
31 }
32 }
33 } // end of moving
34 key_[tmp] = key;
35 event_[tmp] = e;
36 }
37 }
```
computations, the simulated network is examined and replaced by a regenerated topology until it is fully connected. In order to make the experiment more generalized, the simulated network topology is built randomly. The number of network nodes and the pattern of traffics

The Speedup of Discrete Event Simulations by Utilizing CPU Caching 59

**Figure 5.** Computational time and speedup factor over CQ VS traffics (in unit of 99 flows)

Num. of List Calendar Cache Speedup traffic flows (ms) (ms) (ms) factor ×99 2,273 1,517 1,086 1.40 ×99 238,797 19,245 14,455 1.33 ×99 2,523,570 203,784 34,991 5.82 ×99 6,916,930 1,144,140 63,792 17.94 ×99 11,862,000 621,804 99,360 6.26 ×99 17,934,800 1,692,130 140,996 12.00 ×99 26,835,100 1,159,060 192,093 6.03 ×99 38,388,600 3,712,680 257,185 14.44 ×99 48,637,600 2,880,920 330,402 8.72 ×99 - 10,459,600 410,934 25.45 ×99 - 18,573,400 501,229 37.06

**Table 1.** Computational times of three scheduling algorithms and speedups of the proposed (Cache) to

Table 1 and Fig.5 show the computational time spending for simulations *T* versus the number of traffic flows, for CQ (Calendar), the cache aware (Cache) and the linear list (List) schedulers. Experiments are carried out on a personal computer with a Intel(R) Pentium(R) 4 typed CPU with 2.93 GHz in frequency and 1024 kB in cache size. The computational time is evaluated by invoking TCL predefined command clock at 0.0 and 2.0 seconds in the simulation time. It can be seen from Fig.5 that the cache aware algorithm is always faster than CQ and the

are adjustable.

CQ

maximum speedup factor reaches 37.

There are 5 conditional branches in the function insertDigest(), the 1st (omitted at Line 06) and 2nd (omitted at Line 08) handle the buffer of digest with full and empty, respectively. The 3rd case is that the new digest has the same digest time as an existing one, os replacing is required. The 4th is identified in order to avoid the complex operations as defined in the 5th case. The last case involves element movings over an array and can lead to much more memory accesses. Such processing is, however, executing on a memory area that can be allocated continually or locally. The time overhead of insertDigest() can be hence reduced by the benefit of CPU caching, as discussed in the section 3.

The variables precision\_ and size\_ dominate the dynamics of cache aware scheduling. The bigger the precision\_ is, the longer the sub-list of DES event tends to be. The optimal value depends on the population of standing events and distribution in the time domain. However, size\_ can be determined according to the size of CPU caches.

#### **4.2. Experiment environment and results**

For simplicity in performance evaluation, simulation experiments are carried out over a random network with 100 nodes, each node connects to 6 others being selected randomly. Simulation configurations are coded in a Tools Command Language (TCL) script. Every link of the network is configured with fixed bandwidth 155 Mbps and fixed propagation delay 2 ms. A Constant Bit Rate (CBR) generator with demand bandwidth 10 Mbps and packet size 1000 bytes is assigned for each traffic flow and kept active during simulation time from 0.01 to 2.0 seconds. The number of flows varies from 99, i.e. one node is chosen to generate packets to the rest, to 9900, i.e. a full mesh-typed flow pattern is arranged. Before putting into

computations, the simulated network is examined and replaced by a regenerated topology until it is fully connected. In order to make the experiment more generalized, the simulated network topology is built randomly. The number of network nodes and the pattern of traffics are adjustable.

12 Will-be-set-by-IN-TECH

16 } else { // in the middle

19 if (tail\_ > idx) { // right moving 20 for (tmp = tail\_; tmp > idx; tmp--) {

25 } else { // left moving

28 for (tmp = head\_; tmp < idx; tmp++) {

33 } // end of moving

reduced by the benefit of CPU caching, as discussed in the section 3.

However, size\_ can be determined according to the size of CPU caches.

**4.2. Experiment environment and results**

There are 5 conditional branches in the function insertDigest(), the 1st (omitted at Line 06) and 2nd (omitted at Line 08) handle the buffer of digest with full and empty, respectively. The 3rd case is that the new digest has the same digest time as an existing one, os replacing is required. The 4th is identified in order to avoid the complex operations as defined in the 5th case. The last case involves element movings over an array and can lead to much more memory accesses. Such processing is, however, executing on a memory area that can be allocated continually or locally. The time overhead of insertDigest() can be hence

The variables precision\_ and size\_ dominate the dynamics of cache aware scheduling. The bigger the precision\_ is, the longer the sub-list of DES event tends to be. The optimal value depends on the population of standing events and distribution in the time domain.

For simplicity in performance evaluation, simulation experiments are carried out over a random network with 100 nodes, each node connects to 6 others being selected randomly. Simulation configurations are coded in a Tools Command Language (TCL) script. Every link of the network is configured with fixed bandwidth 155 Mbps and fixed propagation delay 2 ms. A Constant Bit Rate (CBR) generator with demand bandwidth 10 Mbps and packet size 1000 bytes is assigned for each traffic flow and kept active during simulation time from 0.01 to 2.0 seconds. The number of flows varies from 99, i.e. one node is chosen to generate packets to the rest, to 9900, i.e. a full mesh-typed flow pattern is arranged. Before putting into

21 key\_[tmp] = key\_[tmp-1]; 22 event\_[tmp] = event\_[tmp-1];

29 key\_[tmp] = key\_[tmp+1]; 30 event\_[tmp] = event\_[tmp+1];

24 tail\_ = ++tail\_ % size\_;

18 if (key\_[idx] > key) { // insert a new element

15 tail\_ = tmp;

23 }

31 } 32 }

36 } 37 }

26 head\_ --; 27 idx --;

34 key\_[tmp] = key; 35 event\_[tmp] = e;

17 idx = ++idx % size\_;

**Figure 5.** Computational time and speedup factor over CQ VS traffics (in unit of 99 flows)


**Table 1.** Computational times of three scheduling algorithms and speedups of the proposed (Cache) to CQ

Table 1 and Fig.5 show the computational time spending for simulations *T* versus the number of traffic flows, for CQ (Calendar), the cache aware (Cache) and the linear list (List) schedulers. Experiments are carried out on a personal computer with a Intel(R) Pentium(R) 4 typed CPU with 2.93 GHz in frequency and 1024 kB in cache size. The computational time is evaluated by invoking TCL predefined command clock at 0.0 and 2.0 seconds in the simulation time. It can be seen from Fig.5 that the cache aware algorithm is always faster than CQ and the maximum speedup factor reaches 37.

Fig.6 shows the computational time varying with the size of digest queue, for the case configured with a cache aware scheduler and 30 sources, i.e. 99 × 30 = 2970 flows. Since a ring-buffer management is used in the cache aware algorithm, the actual capacity of the queue is one smaller than its literal size. The condition with the queue sized 1 means to disable the digest queue, and the cache aware algorithm is degraded and equivalent to the linear list. The condition with the queue sized 2 allows only one digest which always points to the tail of DES event list. Experiment results on the size of digest queue are also listed in Tab.2.

The algorithm provided in this chapter is not the best that can fully utilize the capability of CPU resources since that the design of digest queue and operations have no knowledge of CPU's microstructure. However, a further improvement that depends on the specific CPU products will bring a serious issue in compatibility. The trade-off between performance and compatibility needs extensive investigations and relates to the application field of discrete

The Speedup of Discrete Event Simulations by Utilizing CPU Caching 61

Discrete event simulation method has been employed in various research fields for exploring a complex or large system numerically. The applicability of this kind of simulation method relies in large on the computation speed. In this chapter, fast event scheduling approaches for simulations of large scale networks is discussed. The typical simulation procedures are described according to NS2 platform, followed by a brief analysis of the conventional Calendar Queue algorithm. Based on the list partitioning method, a digest queue on the fixed array structure is introduced to partition and index the sub-list of DES events. The double-list

Details of enqueue and dequeue algorithm are given and the complexity analysis is presented. Also, developments based on the open-source NS2 software are showed for the algorithm implementation. In order to verify the benefit of cache awareness to speedup simulation, we report computational experiments over an on-shelf personal computer. It is shown that the performance of a cache aware algorithm is considerably better than the conventional Calendar Queue. Experiment results show that the cache aware algorithm makes simulation faster by up to a factor of 37. The improvement in computation efficiency is applicable for any ordinary network topology and traffic pattern since the randomized topologies and patterns

For simplicity, the algorithm's design and implementation are focused on the network simulation with all traffics assumed over UDP/IP protocol stacks. As for the simulation involves in TCP connection, events used for protocol machine controlling, such as those for time-out processing, bring numbers of canceling operations into event scheduling. The event cancelation in the queue based on the ring-typed buffer should result in more computational complex than that of convectional algorithms. A better solution, we think, is to put these "at-events" being significant for the local network node into a separate object from the global event scheduler. This separation is not only benefit in speedup to the scheduler but also to

On the other hand, the advancement of multi-core architecture of CPU deserves to be verified whether the proposed algorithm can utilize caches of multi-core chip as well. Moreover, it's

*Key Lab of Broadband Wireless Communication and Sensor Network Technology, Nanjing University*

valuable to study multi-core technology to achieve parallel speedup in the meantime.

structure can utilize the benefit of caching mechanism of modern CPU chips.

event simulations.

**5. Conclusion**

are employed in our experiments.

parallel enabler.

**Author details**

Wennai Wang and Yi Yang

*of Posts and Telecommunications, China*


**Table 2.** Computational time (T) varies with the size of digest queue

From Fig.6 and Tab.2, it can be concluded that a cache aware algorithm makes use of the benefit of CPU caches when the size of digest queue is greater than 16. It is can be seen also that the speedup efficiency relies on the size of cache-line rather than the size and topology of the simulated network. The detailed analysis is given in the next section.

**Figure 6.** Computation time VS digest queue size

#### **4.3. Effectiveness of CPU caching**

The result of the experiment indicates that the cache aware algorithm can speedup DES event scheduling for large scale networks. Continuously allocated digest queue in memory is benefit to fully use the feature of caching. Meanwhile, the digest queue acts as a classifier to partition DES event list into sub-lists with a shorter length, and hence reduce the searching time of ENQUEUE operations.

The digest in the cache aware algorithm is managed by two variables, key\_ and Event\_, each 4 bytes in length. Hence, a digest queue with size 16 occupies 128 bytes in the main memory. This is exactly equivalent to 128-byte sized L2 cache line of the processor [8] used in experiments. The dropping approaching the size 16 shown in Fig.6 coincides with the CPU L2 caching capability. Therefore, the algorithm is cache-line awareness for the CPU chip employed in our experiments. This means that computation time reduction remains nearly constant when the size of digest queue increases over that of cache line, i.e., 8 elements or 64 bytes, roughly 10 as showed in experiments. In this condition, missing in single a cache line causes a regular main memory access so that the contribution of CPU caching is saturated.

The algorithm provided in this chapter is not the best that can fully utilize the capability of CPU resources since that the design of digest queue and operations have no knowledge of CPU's microstructure. However, a further improvement that depends on the specific CPU products will bring a serious issue in compatibility. The trade-off between performance and compatibility needs extensive investigations and relates to the application field of discrete event simulations.

## **5. Conclusion**

14 Will-be-set-by-IN-TECH

Fig.6 shows the computational time varying with the size of digest queue, for the case configured with a cache aware scheduler and 30 sources, i.e. 99 × 30 = 2970 flows. Since a ring-buffer management is used in the cache aware algorithm, the actual capacity of the queue is one smaller than its literal size. The condition with the queue sized 1 means to disable the digest queue, and the cache aware algorithm is degraded and equivalent to the linear list. The condition with the queue sized 2 allows only one digest which always points to the tail of DES

> size 1 2 3 4 8 16 32 64 T(ms) 4,052,043 4,131,210 1,929,621 281,398 64,113 63,186 64,140 63,109

From Fig.6 and Tab.2, it can be concluded that a cache aware algorithm makes use of the benefit of CPU caches when the size of digest queue is greater than 16. It is can be seen also that the speedup efficiency relies on the size of cache-line rather than the size and topology of

The result of the experiment indicates that the cache aware algorithm can speedup DES event scheduling for large scale networks. Continuously allocated digest queue in memory is benefit to fully use the feature of caching. Meanwhile, the digest queue acts as a classifier to partition DES event list into sub-lists with a shorter length, and hence reduce the searching time of

The digest in the cache aware algorithm is managed by two variables, key\_ and Event\_, each 4 bytes in length. Hence, a digest queue with size 16 occupies 128 bytes in the main memory. This is exactly equivalent to 128-byte sized L2 cache line of the processor [8] used in experiments. The dropping approaching the size 16 shown in Fig.6 coincides with the CPU L2 caching capability. Therefore, the algorithm is cache-line awareness for the CPU chip employed in our experiments. This means that computation time reduction remains nearly constant when the size of digest queue increases over that of cache line, i.e., 8 elements or 64 bytes, roughly 10 as showed in experiments. In this condition, missing in single a cache line causes a regular main memory access so that the contribution of CPU caching is saturated.

event list. Experiment results on the size of digest queue are also listed in Tab.2.

the simulated network. The detailed analysis is given in the next section.

**Table 2.** Computational time (T) varies with the size of digest queue

**Figure 6.** Computation time VS digest queue size

**4.3. Effectiveness of CPU caching**

ENQUEUE operations.

Discrete event simulation method has been employed in various research fields for exploring a complex or large system numerically. The applicability of this kind of simulation method relies in large on the computation speed. In this chapter, fast event scheduling approaches for simulations of large scale networks is discussed. The typical simulation procedures are described according to NS2 platform, followed by a brief analysis of the conventional Calendar Queue algorithm. Based on the list partitioning method, a digest queue on the fixed array structure is introduced to partition and index the sub-list of DES events. The double-list structure can utilize the benefit of caching mechanism of modern CPU chips.

Details of enqueue and dequeue algorithm are given and the complexity analysis is presented. Also, developments based on the open-source NS2 software are showed for the algorithm implementation. In order to verify the benefit of cache awareness to speedup simulation, we report computational experiments over an on-shelf personal computer. It is shown that the performance of a cache aware algorithm is considerably better than the conventional Calendar Queue. Experiment results show that the cache aware algorithm makes simulation faster by up to a factor of 37. The improvement in computation efficiency is applicable for any ordinary network topology and traffic pattern since the randomized topologies and patterns are employed in our experiments.

For simplicity, the algorithm's design and implementation are focused on the network simulation with all traffics assumed over UDP/IP protocol stacks. As for the simulation involves in TCP connection, events used for protocol machine controlling, such as those for time-out processing, bring numbers of canceling operations into event scheduling. The event cancelation in the queue based on the ring-typed buffer should result in more computational complex than that of convectional algorithms. A better solution, we think, is to put these "at-events" being significant for the local network node into a separate object from the global event scheduler. This separation is not only benefit in speedup to the scheduler but also to parallel enabler.

On the other hand, the advancement of multi-core architecture of CPU deserves to be verified whether the proposed algorithm can utilize caches of multi-core chip as well. Moreover, it's valuable to study multi-core technology to achieve parallel speedup in the meantime.

## **Author details**

Wennai Wang and Yi Yang *Key Lab of Broadband Wireless Communication and Sensor Network Technology, Nanjing University of Posts and Telecommunications, China*

16 Will-be-set-by-IN-TECH <sup>62</sup> Discrete Event Simulations – Development and Applications **Chapter 3** 

#### **6. References**


URL: *http://doi.acm.org/10.1145/324138.324142*


© 2012 Montevechi et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution,

**Sensitivity Analysis in Discrete-Event** 

Rafael de Carvalho Miranda and Jonathan Daniel Friend

many sectors due to its versatility, flexibility and analysis potential [5, 6].

Additional information is available at the end of the chapter

most desirable output result for simulation models [8].

possess integrated optimization routines.

José Arnaldo Barra Montevechi,

http://dx.doi.org/10.5772/50196

**1. Introduction** 

**Simulation Using Design of Experiments** 

The use of discrete-event simulation as an aid in decision-making has grown over recent decades [1, 2, 3, 4]. It is already used as one of the most utilized research techniques for

However, one of simulation's greatest disadvantages is that, on its own, it does not serve as an optimization technique [7]. This forces simulation practitioners to simulate multiple system configurations and choose the one which presents the best system performance. Computational development has helped alter this scenario due to the increasing availability of faster computers and ever improving search and heuristic optimization techniques.

Simulation optimization can be defined as the process of testing different combinations of values for controllable values, aiming to find the combination of values which offers the

In support of this claim, [1, 9, 10, 11] assert that using optimization along with simulation has been continuously increasing due to the emergence of simulation packages which

The overarching idea of including these routines is to search for improved definitions for the system parameters in relation to its performance. However, according to [10], at the end of

Despite the fact that simulation has been around for more than half a century, until quite recently, the scientific community was reluctant to use optimization tools in simulation. The first time the subject emerged in two renowned simulation books, [12] and [13], was at the close of the 20th century [9]. This resistance has begun to diminish with the convent of meta-

and reproduction in any medium, provided the original work is properly cited.

optimization, the user has no way of knowing if an optimal point was truly reached.

heuristic research, along with strides being made in statistic analysis [14].

## **Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments**

José Arnaldo Barra Montevechi, Rafael de Carvalho Miranda and Jonathan Daniel Friend

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50196

## **1. Introduction**

16 Will-be-set-by-IN-TECH

[1] Ahn, J. & Seunghyun, O. [1999]. Dynamic calendar queue, *Annual Simulation Symposium*,

[2] Askitis, N. & Sinha, R. [2007]. Hat-trie: a cache-conscious trie-based data structure for strings, *Proc. of the thirtieth Australasian conference on Computer science - Volume 62*,

[3] Banks, J. [1999]. Introduction to simulation, *Proceedings of the 31st conference on Winter simulation: Simulation—a bridge to the future - Volume 1*, WSC '99, ACM, New York, NY,

[4] Brown, R. [1988]. Calendar queues: A fast O(1) priority queue implementation for the

[5] Chiueh, T. & Pradhan, P. [1999]. High performance ip routing table lookup using cpu

[6] Chung, K., Sang, J. & Rego, V. [1993]. A performance comparison of event calendar

[7] Fujimoto, R. [1999]. Parallel and distributed simulation., *Winter Simulation Conference'99*,

[8] Hinton, G., Sager, D., Upton, M., Boggs, D., Carmean, D., Kyker, A. & Roussel, P. [2001]. The microarchitecture of the pentium 4 processor, *Intel Technology Journal* 5(1): 1–13. [9] Issariyakul, T. & Hossain, E. [2008]. *Introduction to Network Simulator NS2*, 1 edn, Springer

[11] Park, H. & Fishwick, P. A. [2011]. An analysis of queuing network simulation using gpu-based hardware acceleration, *ACM Trans. Model. Comput. Simul.* 21(18): 1–22. [12] Siangsukone, T., Aswakul, C. & Wuttisittikulkij, L. [2003]. Study of optimised bucket widths in calendar queue for discrete event simulator, *Proc. of Thailand's Electrical*

[13] Tan, K. L. & Thng, L. [2000]. Snoopy calendar queue, *Winter Simulation Conference*,

[14] Varga, A. & Hornig, R. [2008]. An overview of the omnet++ simulation environment, *Simutools '08: Proceedings of the 1st international conference on Simulation tools and techniques for communications, networks and systems & workshops*, ICST, Brussels, Belgium, pp. 1–10. [15] Yan, G. & Eidenbenz, S. [2006]. Sluggish calendar queues for network simulation, *Proc. of the 14th IEEE International Symposium on Modeling, Analysis, and Simulation of Computer*

*and Telecommunication Systems (MASCOTS'06)*, Monterey, CA, USA, pp. 127–136.

simulation event set problem, *Communications of the ACM* 31(10): 1220–1227.

algorithms: an empirical approach, *Softw. Pract. Exper.* 23: 1107–1138.

Australian Computer Society, Inc., Darlinghurst, Australia, pp. 97–105.

URL: *http://doi.acm.org/10.1145/324138.324142*

caching, *INFOCOM*, pp. 1421–1428.

Publishing Company, Incorporated. [10] Fall, K. & Varadhan, K. [2009]. The ns manual.

URL: *http://www.isi.edu/nsnam/ns/ns-documentation.html*

*Engineering Conference (EECON-26)*, Phetchaburi, pp. 6–7.

**6. References**

pp. 20–25.

USA, pp. 7–13.

pp. 122–131.

pp. 487–495.

The use of discrete-event simulation as an aid in decision-making has grown over recent decades [1, 2, 3, 4]. It is already used as one of the most utilized research techniques for many sectors due to its versatility, flexibility and analysis potential [5, 6].

However, one of simulation's greatest disadvantages is that, on its own, it does not serve as an optimization technique [7]. This forces simulation practitioners to simulate multiple system configurations and choose the one which presents the best system performance. Computational development has helped alter this scenario due to the increasing availability of faster computers and ever improving search and heuristic optimization techniques.

Simulation optimization can be defined as the process of testing different combinations of values for controllable values, aiming to find the combination of values which offers the most desirable output result for simulation models [8].

In support of this claim, [1, 9, 10, 11] assert that using optimization along with simulation has been continuously increasing due to the emergence of simulation packages which possess integrated optimization routines.

The overarching idea of including these routines is to search for improved definitions for the system parameters in relation to its performance. However, according to [10], at the end of optimization, the user has no way of knowing if an optimal point was truly reached.

Despite the fact that simulation has been around for more than half a century, until quite recently, the scientific community was reluctant to use optimization tools in simulation. The first time the subject emerged in two renowned simulation books, [12] and [13], was at the close of the 20th century [9]. This resistance has begun to diminish with the convent of metaheuristic research, along with strides being made in statistic analysis [14].

© 2012 Montevechi et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

According to [15], verification of system performance for a determined set of system parameters with reasonable precision using simulation, demands a considerable amount of computational potential. In order to find the optimal or near-optimal solution, one needs to verify a large number of parameter values, thus, optimization via simulation is normally exhaustive from a computational standpoint.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 65

**y1 y2 ym ...**

݂ሺߠሻ (1)

Optimization helps respond to the following questions: What are the optimal adjustments to the input variables (x) which maximize (or minimize) a given simulation model output? The objective is to find an optimal value which maximizes or minimizes a determined

**Simulation Model**

**Input Output**

According to [21], simulation optimization is one of the most important technologies to come about in recent years. These authors recall that previous methodologies demanded carrying out complex changes to the simulation model, thus consuming time and computational potential and, in many cases, not even being economically viable for real

A traditional simulation optimization problem (minimization of a single objective) is given

߆ א ߠݐ Ǥݏ Where ݂ሺߠሻ ൌ ܧሾ߰ሺߠǡ ߱ሻሿ is the system's expected performance, estimated for ݂መሺߠሻ, which is obtained using the simulation model samples ߰ሺߠǡ ߱ሻ, observed according to the discrete or

According to [23], the optimization method which serves for the problem presented in Eq. (1) depends on whether the simulated variables are discrete or continuous. There are many methods for resolving this problem in the literature, like the one presented in Eq. (1), and unfortunately depending on the model being optimized, some methods cannot guarantee

Table 1 shows the main optimization software packages which are both on the market and cited in academic literature, as well as the simulation packages with which they are sold.

As shown in Table 1, different optimization software packages utilize different search methods, such as: Evolutionary Algorithms [25], Genetic Algorithms [26], Scatter Search

According to [31] and [32], the simulation optimization's greatest limitation is the number of variables being manipulated, as the software's performance is considerably reduced in models with a great number of variables. Thus, [33] asserts that convergence time is the most significant restriction in reaching computational efficiency for optimization algorithms.

continuous input parameters, restricted by θ within a viability set ؿ߆ॆௗ.

The optimization techniques utilized in each software package is also shown.

[27], Taboo Search [28], Neural Networks [29] and Simulated Annealing [30].

**Figure 1.** Simulation Model [20]

performance indicator [11].

in Eq. 1 [22]:

cases due to the large number of decision variables.

**x1 x2 xn ...**

that an optimal solution is found [24].

Having highlighted the computational strains, [8] states that despite the evolution of optimization software, a common criticism made about such commercial packages is that, when more than one variable is manipulated at a time, the software becomes very slow.

Considering that not all decision variables are of equal importance in respect to the response variable that one desires to optimize [16, 17], a sensitivity analysis may be carried out on the simulation model in order to select the variables which will compose the optimization search space in order to limit the number of variables and, in turn, make the search faster.

Thus, in order to proceed to the variable selection, screening can be done in order to separate the most important variables from those which may be eliminated from consideration [16, 17]. The same author presents some examples of experimental design utilized in screening experiments:


The current chapter presents an application of Design of Experiments (DOE), specifically fractional factorial design, in order to select the significant input variables in a simulation model, and thus accelerate the optimization process. For information about experimental design, the reader can consult [1, 4, 18, 19].

Fractional factorial design is a DOE technique in which only a fraction of the total number of experiments is executed, thus realizing fewer experiments than full factorial design. Throughout this chapter, it is shown that the use of such a design serves to reduce the search space in the optimization phase of simulation studies.

In this chapter, real examples of how to conduct sensitivity analysis with factorial design are given. In order to reach this goal, two study objects are presented, comparing the optimization carried out without previous investigation of input variable significance, with the optimization carried out in reduced search space. Finally, a comparison is made of the results of the optimization, with and without the sensitivity analysis.

### **2. Simulation optimization**

A simulation model generally includes *n* input variables (x1, x2,...,xn) and *m* output variables (y1, y2, ..., ym) (Figure 1). The optimization of this simulation method implies finding the optimal configuration of input variables; that is, the values of x1, x2, …, xn which optimize the response variable(s) [20].

**Figure 1.** Simulation Model [20]

exhaustive from a computational standpoint.

utilized in screening experiments:

 2n-p fractional factorial design; Supersaturated designs; Groups screening designs.

**2. Simulation optimization** 

the response variable(s) [20].

design, the reader can consult [1, 4, 18, 19].

search space in the optimization phase of simulation studies.

results of the optimization, with and without the sensitivity analysis.

2n factorial design;

According to [15], verification of system performance for a determined set of system parameters with reasonable precision using simulation, demands a considerable amount of computational potential. In order to find the optimal or near-optimal solution, one needs to verify a large number of parameter values, thus, optimization via simulation is normally

Having highlighted the computational strains, [8] states that despite the evolution of optimization software, a common criticism made about such commercial packages is that, when more than one variable is manipulated at a time, the software becomes very slow.

Considering that not all decision variables are of equal importance in respect to the response variable that one desires to optimize [16, 17], a sensitivity analysis may be carried out on the simulation model in order to select the variables which will compose the optimization search space in order to limit the number of variables and, in turn, make the search faster.

Thus, in order to proceed to the variable selection, screening can be done in order to separate the most important variables from those which may be eliminated from consideration [16, 17]. The same author presents some examples of experimental design

The current chapter presents an application of Design of Experiments (DOE), specifically fractional factorial design, in order to select the significant input variables in a simulation model, and thus accelerate the optimization process. For information about experimental

Fractional factorial design is a DOE technique in which only a fraction of the total number of experiments is executed, thus realizing fewer experiments than full factorial design. Throughout this chapter, it is shown that the use of such a design serves to reduce the

In this chapter, real examples of how to conduct sensitivity analysis with factorial design are given. In order to reach this goal, two study objects are presented, comparing the optimization carried out without previous investigation of input variable significance, with the optimization carried out in reduced search space. Finally, a comparison is made of the

A simulation model generally includes *n* input variables (x1, x2,...,xn) and *m* output variables (y1, y2, ..., ym) (Figure 1). The optimization of this simulation method implies finding the optimal configuration of input variables; that is, the values of x1, x2, …, xn which optimize Optimization helps respond to the following questions: What are the optimal adjustments to the input variables (x) which maximize (or minimize) a given simulation model output? The objective is to find an optimal value which maximizes or minimizes a determined performance indicator [11].

According to [21], simulation optimization is one of the most important technologies to come about in recent years. These authors recall that previous methodologies demanded carrying out complex changes to the simulation model, thus consuming time and computational potential and, in many cases, not even being economically viable for real cases due to the large number of decision variables.

A traditional simulation optimization problem (minimization of a single objective) is given in Eq. 1 [22]:

$$\min f(\theta)$$
 
$$\text{s.t.} \; \theta \in \Theta$$

Where ݂ሺߠሻ ൌ ܧሾ߰ሺߠǡ ߱ሻሿ is the system's expected performance, estimated for ݂መሺߠሻ, which is obtained using the simulation model samples ߰ሺߠǡ ߱ሻ, observed according to the discrete or continuous input parameters, restricted by θ within a viability set ؿ߆ॆௗ.

According to [23], the optimization method which serves for the problem presented in Eq. (1) depends on whether the simulated variables are discrete or continuous. There are many methods for resolving this problem in the literature, like the one presented in Eq. (1), and unfortunately depending on the model being optimized, some methods cannot guarantee that an optimal solution is found [24].

Table 1 shows the main optimization software packages which are both on the market and cited in academic literature, as well as the simulation packages with which they are sold. The optimization techniques utilized in each software package is also shown.

As shown in Table 1, different optimization software packages utilize different search methods, such as: Evolutionary Algorithms [25], Genetic Algorithms [26], Scatter Search [27], Taboo Search [28], Neural Networks [29] and Simulated Annealing [30].

According to [31] and [32], the simulation optimization's greatest limitation is the number of variables being manipulated, as the software's performance is considerably reduced in models with a great number of variables. Thus, [33] asserts that convergence time is the most significant restriction in reaching computational efficiency for optimization algorithms.


Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 67

Determine the values of x (significant variables) in order that the effect of the non-

The experimentation strategy is the method of designing and conducting experiments [18]. According to this author, there are many methods which can be used for the realization of

**Best-guess**: This strategy is based on the specialists' technical or theoretical knowledge, which alters the value of one or two variables for the test in function of the previous result. This procedure presents at least two disadvantages: The first disadvantage occurs when the initial configuration does not produce the desired result and then the analyst must search for another input variable configuration. These attempts may continue indefinitely and certainly take a long time without guaranteeing success. The second disadvantage is that, supposing an acceptable initial configuration, the analyst will be tempted to stop testing,

**One factor at a time**: This strategy involves selecting the starting configuration for each input variable and then successively varying each variable within a given range, while simultaneously maintaining the other variables constant. The greatest disadvantage of this strategy is its inability to detect interaction between variables; nonetheless, many analysts

**Factorial Design**: According to [18], when an experiment involves the study of two or more factors, the most effective strategy is factorial design. In using this strategy, factors are altered at the same time, instead of one at a time. That is, in each complete attempt or experimental replica, all possible combinations are investigated [35]. This strategy is more efficient than the one previously mentioned, as it allows for the effects of a single factor to be estimated across various factor levels, thus leading to valid conclusions within the experimental conditions, [18] and is the only way to discover interaction between factors [18, 35], avoiding incorrect conclusions when there is interaction between factors. The main problem with factorial design is the exponentially increasing number of combinations with

even though there is no guarantee that the best results has been obtained.

controllable variable effects are minimized.

**Figure 2.** General process model [34]

experiments. Some examples are listed below:

disregard this fact and it is often used [18].

each increase in the number of factors [19].

**Table 1.** Optimization software packages [4, 7, 9]

In order to ease this process, the use of fractional factorial design can be used to conduct sensitivity analysis on a simulation model in order to select the input variables which truly impact the response variable and enable the elimination of variables which are not statistically significant. In terms of the simulation, sensitivity analysis may be interpreted as a systematic investigation of the model's outputs, in accordance with the model's input variables [19].

By using DOE techniques, it is possible to reduce the number of experiments executed, determine which independent variables affect the dependent variable, and identify the amplitude or intensity of this effect. For optimization purposes, identification of the most significant variables is important, as the greater the number of variables in the search space, the longer the optimization process will take.

Thus, by using sensitivity analysis in simulation optimization problems, one can work with those input variables which actually have a significant effect over the determined response variable, thus reducing the number of experiments necessary and the computational potential involved in this process.

## **3. Experimentation strategies**

An experiment can be defined as a test or series of tests in which purposeful changes are made to input variables, with the objective of observing and identifying the form in which the system responses are affected, in function of the changes carried out on the input variables [18].

According to [34], there are two types of process variables (Figure 2): controllable variables (x1, x2, …, xp), and non-controllable variables (z1, z2, …, zq), which are many times called "sound". The same author states that the experiment's objectives can be:


 Determine the values of x (significant variables) in order that the effect of the noncontrollable variable effects are minimized.

**Figure 2.** General process model [34]

66 Discrete Event Simulations – Development and Applications

AutoStat® AutoMod®,

**Table 1.** Optimization software packages [4, 7, 9]

the longer the optimization process will take.

potential involved in this process.

**3. Experimentation strategies** 

**Software Simulation Package Optimization Technique** 

OptQuest® Arena® Scatter Search, Taboo and Neural

Optimizer® Witness® Simulated Annealing e Taboo Search SimRunner® ProModel® Evolutionary and Genetic Algorithms

In order to ease this process, the use of fractional factorial design can be used to conduct sensitivity analysis on a simulation model in order to select the input variables which truly impact the response variable and enable the elimination of variables which are not statistically significant. In terms of the simulation, sensitivity analysis may be interpreted as a systematic investigation of the model's outputs, in accordance with the model's input

By using DOE techniques, it is possible to reduce the number of experiments executed, determine which independent variables affect the dependent variable, and identify the amplitude or intensity of this effect. For optimization purposes, identification of the most significant variables is important, as the greater the number of variables in the search space,

Thus, by using sensitivity analysis in simulation optimization problems, one can work with those input variables which actually have a significant effect over the determined response variable, thus reducing the number of experiments necessary and the computational

An experiment can be defined as a test or series of tests in which purposeful changes are made to input variables, with the objective of observing and identifying the form in which the system responses are affected, in function of the changes carried out on the input

According to [34], there are two types of process variables (Figure 2): controllable variables (x1, x2, …, xp), and non-controllable variables (z1, z2, …, zq), which are many times called

Determine the values of x (significant variables) in order that the response is close to the

Determine the values of x (significant variables) in order that the variability in y is

"sound". The same author states that the experiment's objectives can be:

Determine the variables which have the most influence over the response (y);

Optimiz® Simul8® Neural Networks

AutoSched® Evolutionary and Genetic Algorithms

Networks

**Optimization** 

variables [19].

variables [18].

small;

nominal demand;

The experimentation strategy is the method of designing and conducting experiments [18]. According to this author, there are many methods which can be used for the realization of experiments. Some examples are listed below:

**Best-guess**: This strategy is based on the specialists' technical or theoretical knowledge, which alters the value of one or two variables for the test in function of the previous result. This procedure presents at least two disadvantages: The first disadvantage occurs when the initial configuration does not produce the desired result and then the analyst must search for another input variable configuration. These attempts may continue indefinitely and certainly take a long time without guaranteeing success. The second disadvantage is that, supposing an acceptable initial configuration, the analyst will be tempted to stop testing, even though there is no guarantee that the best results has been obtained.

**One factor at a time**: This strategy involves selecting the starting configuration for each input variable and then successively varying each variable within a given range, while simultaneously maintaining the other variables constant. The greatest disadvantage of this strategy is its inability to detect interaction between variables; nonetheless, many analysts disregard this fact and it is often used [18].

**Factorial Design**: According to [18], when an experiment involves the study of two or more factors, the most effective strategy is factorial design. In using this strategy, factors are altered at the same time, instead of one at a time. That is, in each complete attempt or experimental replica, all possible combinations are investigated [35]. This strategy is more efficient than the one previously mentioned, as it allows for the effects of a single factor to be estimated across various factor levels, thus leading to valid conclusions within the experimental conditions, [18] and is the only way to discover interaction between factors [18, 35], avoiding incorrect conclusions when there is interaction between factors. The main problem with factorial design is the exponentially increasing number of combinations with each increase in the number of factors [19].

**Response Surface Methodology (MSR)**: This method consists in a set of mathematical and statistical techniques that are used for modeling and analysis, in which the response of interest is influenced by multiple variables, and the objective is to optimize this response [35]. According to these authors, the relation between the independent and dependent response variables is unknown in most problems. [35] states that the first step of MSR is finding an accurate approximation for the true relation between the response (y) and the independent variables. In general, polynomials with a low degree are used to model a given region of the independent variables.

## **4. Simulated experiment design**

According to [36], although classic experimental design methods were developed for real world experiments, they are perfectly applicable to simulated experiments. In fact, according to the same author, simulated experiment design presents many opportunities for improvements which are difficult or impossible to carry out using actual experiments.

[37] asserts that research related to experimental design are frequently found in specialized publications, but they are rarely read by simulation practitioners. According to these same authors, most simulation practitioners can get more from their analyses by using DOE theory developed specifically to experiment with computational models.

The benefits of DOE enable the possibility of improving simulation process performance by avoiding trial-and-error searches for solutions [38]. More specifically, the use of factorial design can minimize or even eliminate the disadvantages brought about by experimenting with simulated systems instead of the real system.

According to [36], in order to facilitate the understanding of simulation's role in experimental execution, it is necessary to imagine a response value (Y) or a dependent value variable can be represented in the following equation:

$$\mathbf{Y} = \mathbf{f} \text{ (хі, хz, ..., х\_n)}\tag{2}$$

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 69

Experimentation using simulation presents some special advantages over using physical or

By using simulation, it is possible to control factors that, in reality, are uncontrollable,

By using simulation, it is possible to control the basic origin of variation, which is

Another experimental design characteristic is that commercial simulators come with random number generators and therefore, from an experimental point of view, the trouble of randomizing the experimental replicas is eliminated. Randomization is a problem with

According to [18], DOE can be defined as the process of designing experiments in order that the appropriate data are collected and analyzed by statistical methods, thus resulting in valid and objective conclusions. Any experimental problem must contain two elements:

DOE techniques are seen with a broad range of application in many knowledge areas, thus showing itself as a set of tools of great importance for process and product development.

Those involved in the research should have a previous idea of the experiment's objective, which factors will be studied, how the experiment will be conducted and a comprehension

1. **Problem recognition and definition:** Completely develop all ideas about the problem and the objectives to be attained through the experiment, thus contributing to greater

2. **Choice of factors and working levels:** Choose the factors to undergo alterations, the intervals of these factors and the specific levels for each run to be carried out; 3. **Selection of the response variables:** Determine the response variables which really supply useful information about the performance of the process under study; 4. **Selection of the experimental design:** Consider the sample size (number of replications), selection of the correct order of runs for the experimental attempts, or the

5. **Realization of experiments:** Monitor the process to guarantee that everything is being completed according to the design – errors in this stage can destroy the experiment's

6. **Statistical data analysis:** Analyze the data using statistical methods, given that results and conclusions are objective and not the outcome of opinions – residuals analysis and

different from physical experiments, thus avoiding use of blocks.

industrial systems [4]:

such as client arrival rate;

physical experimentation [36].

**5. Design and analysis of experiments** 

experimental design and statistical data analysis.

According to [18], DOE should consider the following stages:

comprehension of the process and eventual problem solution;

formation of blocks and other randomization restrictions involved;

verification of model validity are important to this phase;

of how the data will be analyzed [34].

validity;

Where:


[39] declares that simulation is a black box which transforms inputs variables into simulated outputs, which imitate the real system's perspective output. For each scenario, the analyst carries out one or more runs and registers the average output values.

In simulation models, the levels chosen for each factor must enable the effects to be programmed in the model. In order to exemplify this question, the following situation is proposed: a determined factor which is desired to be optimized corresponds to the possibility of using an experienced employee (upper level) or a new hire (lower level), thus verifying, what the impact would be on daily throughput. In simulation models, the modeler must be familiar with each variable to be affected by the change in levels. Thus, the modeler must decide which distribution to use for each variable time in the operation.

Experimentation using simulation presents some special advantages over using physical or industrial systems [4]:


Another experimental design characteristic is that commercial simulators come with random number generators and therefore, from an experimental point of view, the trouble of randomizing the experimental replicas is eliminated. Randomization is a problem with physical experimentation [36].

## **5. Design and analysis of experiments**

68 Discrete Event Simulations – Development and Applications

region of the independent variables.

**4. Simulated experiment design** 

with simulated systems instead of the real system.

variable can be represented in the following equation:

Where:

**Response Surface Methodology (MSR)**: This method consists in a set of mathematical and statistical techniques that are used for modeling and analysis, in which the response of interest is influenced by multiple variables, and the objective is to optimize this response [35]. According to these authors, the relation between the independent and dependent response variables is unknown in most problems. [35] states that the first step of MSR is finding an accurate approximation for the true relation between the response (y) and the independent variables. In general, polynomials with a low degree are used to model a given

According to [36], although classic experimental design methods were developed for real world experiments, they are perfectly applicable to simulated experiments. In fact, according to the same author, simulated experiment design presents many opportunities for improvements which are difficult or impossible to carry out using actual experiments.

[37] asserts that research related to experimental design are frequently found in specialized publications, but they are rarely read by simulation practitioners. According to these same authors, most simulation practitioners can get more from their analyses by using DOE

The benefits of DOE enable the possibility of improving simulation process performance by avoiding trial-and-error searches for solutions [38]. More specifically, the use of factorial design can minimize or even eliminate the disadvantages brought about by experimenting

According to [36], in order to facilitate the understanding of simulation's role in experimental execution, it is necessary to imagine a response value (Y) or a dependent value

[39] declares that simulation is a black box which transforms inputs variables into simulated outputs, which imitate the real system's perspective output. For each scenario, the analyst

In simulation models, the levels chosen for each factor must enable the effects to be programmed in the model. In order to exemplify this question, the following situation is proposed: a determined factor which is desired to be optimized corresponds to the possibility of using an experienced employee (upper level) or a new hire (lower level), thus verifying, what the impact would be on daily throughput. In simulation models, the modeler must be familiar with each variable to be affected by the change in levels. Thus, the modeler must decide which distribution to use for each variable time in the operation.

Y = f (x1, x2, ..., xn) (2)

theory developed specifically to experiment with computational models.

x1, x2, xn represent the input variables, factors or independent variables;

*f* represents the simulation model's transformation function.

carries out one or more runs and registers the average output values.

According to [18], DOE can be defined as the process of designing experiments in order that the appropriate data are collected and analyzed by statistical methods, thus resulting in valid and objective conclusions. Any experimental problem must contain two elements: experimental design and statistical data analysis.

DOE techniques are seen with a broad range of application in many knowledge areas, thus showing itself as a set of tools of great importance for process and product development.

Those involved in the research should have a previous idea of the experiment's objective, which factors will be studied, how the experiment will be conducted and a comprehension of how the data will be analyzed [34].

According to [18], DOE should consider the following stages:

	- 7. **Conclusions and recommendations:** Provide practical conclusions based on the results and recommend a plan of action. Accompanying sequences and confirmation tests must be conducted in order to validate the experiment's conclusions.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 71

According to [35], in order to test if the alteration in one of the levels or interaction is significant, a hypothesis test for the average can be used. In the case of DOE, this test can be conducted using ANOVA. The statistical test ANOVA is utilized to accept or reject hypotheses investigated with DOE. Its objective is to analyze the average variation of results and demonstrate which factors actually produce significant events over the system response

However, according to [18], it is not advisable to trust solely in ANOVA since the validity of its suppositions may be unreliable. Problems with these results may be identified using

Residual analysis is an important procedure to guarantee that the models developed by means of experimentation adequately represent the responses of interest. [18] defines residuals as the difference between the predicted value and the observed experimental value; the same author also asserts that residues should be normal, random and non-

Full factorial design with two levels or factorial 2k is a type of design in which two levels are defined for each factor, an upper and lower level, and combinations of factors are tested [8]. 2k factorial design is one of the most important types of factorial design, according to [35], and can be particularly useful in the initial phases of experimental work, especially when many factors are being investigated. Full factorial design offers the fewest executions for the

In full factorial design, the number of experiments is equal to the number of experimental levels, elevated to the number of factors. In the case of factorials with two levels, the number of experiments (N) in order to evaluate k factors is given by N = 2k. These designs possess a

In using this strategy, the factors are altered simultaneously and not just one at a time, which indicates that, for each run or complete replica, all possible combinations are investigated [35]. For example, if there are *a* levels for factor *A* and *b* levels for factor *B*, then

One aspect to be considered is that, as there are only two levels per factor, it must be supposed that the response is approximately linear within the range of levels chosen [35]. Another important aspect is that, for experiments with a great number of factors, full factorial design results in an extremely large number of combinations. In this situation, fractional factorial planning is used in order to select a subset of combinations within the full factorial design, aiming to identify the significant factors in system performance [8].

According to [39], many studies in operational research use full factorial design due to its simplicity and because the technique allows the analyst to identify interactions between

simplified analysis and form the base of many other experimental designs [34].

variables [42].

residual analysis.

correlated.

**5.2. Full factorial design** 

k factors to be studied.

each replica will contain *ab* combinations [18].

factors as well as their main effects.

Stages 1 – 3 are commonly called the pre-experimental design and, for the experiment's success, it is important that these steps are carried out in the most appropriate manner possible [34].

## **5.1 DOE: Main concepts**

There are three basic principles to DOE [18]:


Now that the basic principles of DOE have been defined, the following list presents some of the fundamental terms which are used when dealing with DOE techniques:


Aside from these commonly utilized experimental design terms, two further important concepts should be presented: Analysis of variance (ANOVA) and residuals analysis.

According to [35], in order to test if the alteration in one of the levels or interaction is significant, a hypothesis test for the average can be used. In the case of DOE, this test can be conducted using ANOVA. The statistical test ANOVA is utilized to accept or reject hypotheses investigated with DOE. Its objective is to analyze the average variation of results and demonstrate which factors actually produce significant events over the system response variables [42].

However, according to [18], it is not advisable to trust solely in ANOVA since the validity of its suppositions may be unreliable. Problems with these results may be identified using residual analysis.

Residual analysis is an important procedure to guarantee that the models developed by means of experimentation adequately represent the responses of interest. [18] defines residuals as the difference between the predicted value and the observed experimental value; the same author also asserts that residues should be normal, random and noncorrelated.

## **5.2. Full factorial design**

70 Discrete Event Simulations – Development and Applications

There are three basic principles to DOE [18]:

of statistical methods in experimental design;

homogeneity of experimental conditions.

**Levels:** The variations possible for each factor [41];

quantitative or qualitative;

inferior to superior level;

input factors [8];

possible [34].

effect.

**5.1 DOE: Main concepts** 

7. **Conclusions and recommendations:** Provide practical conclusions based on the results and recommend a plan of action. Accompanying sequences and confirmation tests must

Stages 1 – 3 are commonly called the pre-experimental design and, for the experiment's success, it is important that these steps are carried out in the most appropriate manner

 **Randomization:** Execution of experiments in a random order in order that the phenomenon's unknown effects are distributed among the factors, thus increasing the investigation's validity. According to the author, randomization is the base for the use

 **Replication:** Repetition of the same test multiple times, thus creating a variation in the response variable which is utilized to evaluate experimental error. With the use of each replication, it is possible to obtain an estimate of experimental error, allowing for the determination of whether the differences observed in the data are statistically different, as well as obtaining a more accurate estimate of an experimental factor's

 **Blocking:** Design technique used to increase precision with the comparisons between the factors of interest. It is frequently utilized to reduce or eliminate variability transmitted by factors of sound. It should be utilized when it is not possible to maintain

Now that the basic principles of DOE have been defined, the following list presents some of

 **Factor:** According to [37], factors are input parameters and the structural considerations which compose an experiment. Factors are altered during experimental conduction. According to [40], a factor may assume at least two values during an experiment, being

 **Main effect:** According to [36], the main effect for a factor may be defined as the average of the differences in the response variable, when the factor changes from an

 **Response variable:** The response variable is the performance measure for the DOE. The response variables describe how the system responds under a certain configuration of

 **Interaction:** There is interaction between the factors when the response difference between the levels of a given factor is not the same as for the rest of the factors.

Aside from these commonly utilized experimental design terms, two further important

concepts should be presented: Analysis of variance (ANOVA) and residuals analysis.

the fundamental terms which are used when dealing with DOE techniques:

be conducted in order to validate the experiment's conclusions.

Full factorial design with two levels or factorial 2k is a type of design in which two levels are defined for each factor, an upper and lower level, and combinations of factors are tested [8]. 2k factorial design is one of the most important types of factorial design, according to [35], and can be particularly useful in the initial phases of experimental work, especially when many factors are being investigated. Full factorial design offers the fewest executions for the k factors to be studied.

In full factorial design, the number of experiments is equal to the number of experimental levels, elevated to the number of factors. In the case of factorials with two levels, the number of experiments (N) in order to evaluate k factors is given by N = 2k. These designs possess a simplified analysis and form the base of many other experimental designs [34].

In using this strategy, the factors are altered simultaneously and not just one at a time, which indicates that, for each run or complete replica, all possible combinations are investigated [35]. For example, if there are *a* levels for factor *A* and *b* levels for factor *B*, then each replica will contain *ab* combinations [18].

One aspect to be considered is that, as there are only two levels per factor, it must be supposed that the response is approximately linear within the range of levels chosen [35]. Another important aspect is that, for experiments with a great number of factors, full factorial design results in an extremely large number of combinations. In this situation, fractional factorial planning is used in order to select a subset of combinations within the full factorial design, aiming to identify the significant factors in system performance [8].

According to [39], many studies in operational research use full factorial design due to its simplicity and because the technique allows the analyst to identify interactions between factors as well as their main effects.

Factorial designs are more efficient than the one-at-a-time approach, as they allow for the factors' effects to be estimated via the other factors' levels, thus leading to valid conclusions about the experimental scope; they are also the only manner to discover interactions among the variables, thus avoiding erroneous conclusions when interactions between the factors are present [18].

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 73

 **Design resolution V:** These are the designs in which no main effect or second order interaction is associated with any other main effect or second order interaction, but the

As a way of simplifying DOE application in sensitivity analysis of simulation models, the

The flowchart proposed in Figure 3 presents the necessary steps for conducting sensitivity analysis in discrete-event simulation models. Four stages are defined for the proposed method:

In the first step, the analyst should define the optimization objectives as well as verify if the simulation model is verified and validated. In doing so, the model will be apt to proceed to the following step, fractional factorial design. In this phase, the analyst should determine the

Once these initial steps have been completed, fractional factorial design can be applied. During execution of the experiments, the analyst should return to the simulation step and carry out the experiments in the simulation package. With the experiments done, the data should be analyzed using statistical means, determining the factors' significance level as well as its lower order interactions. At the end of this phase, the non-significant factors can

The third stage defined in the sensitivity analysis phase may be omitted by the analyst, depending on the degree of precision desired or for the cases in which the simulation models demand a lot of computational time in order to be processed. In this stage, a full factorial design is generated with experimental data, and only factors which show to be significant in the previous steps are tested. Depending on the necessity for executing more experiments in order to comply with full factorial design, the analyst will have to return to the simulation phase again to execute new experiments. In this stage, the residues should be analyzed in order to validate the results. Statistical analysis should be conducted once again

In the following stage, optimization simulation is utilized again. The significant factors after full factorial design application are utilized to configure the optimization tool, which is then executed. Many different configurations are tested for the input parameters until the optimizer converges on a solution. It is up to the analyst to evaluate the results and generate

In order to demonstrate the utilization of this method in sensitivity analysis of discrete-event simulation models, two simulation models will be used as study objects in this chapter.

model's input factors, their levels and select the response variables for analysis.

second order interactions are associated with third order interactions.

**6. Sensitivity analysis development stages** 

following sequence of steps is proposed.

2. Fractional Factorial Design; 3. Full Factorial Design;

be removed from analysis.

in order to finalize this stage.

his or her conclusions and recommendations.

1. Simulation;

4. Optimization.

## **5.3. Fractional factorial designs**

When there is little interest in interaction behavior among the factors which compose the system, this can be disregarded [35]. Instead, fractional factorial design can be used.

For example, consider a factorial design of 25. In this planning, five degrees of freedom correspond to the main effects, 10 degrees of freedom correspond to second order interactions and 16 degrees of freedom correspond to the highest order of interactions. In initial system or project studies, there is little interest in the highest level of interactions [35].

If interactions can be disregarded, fractional factorial design involving fewer executions for a complete set of 2k executions can be used in order to obtain information about the main effects and lower order interactions [35].

Thus, fractional factorial design provides a means by which to obtain estimates of main effects and, perhaps, second order interactions, with a fraction of the computational force required for full factorial design [4].

According to [18], the greatest application of fractional factorial designs is in screening experiments, where many factors are present in a system and the objective is to identify which factors indeed exercise a significant effect over the given response of interest. For the factors identified as significant through the use of fractional designs, the author recommends a more careful analysis with the use of other designs, such as full factorial design.

In fractional factorial design, a subset of 2k-p is constructed from a set of all of the possible points for a 2k design, and a simulation is executed for only the chosen points [4].

For this type of factorial design, the analyst must be attentive to its resolution. According to [35], the concept of resolution design is the form in which the fractional factorial designs are related in accordance with the associative standards which they produce. A design's resolution is represented by a Roman subscript number; for example, 2��� �����represents a factorial design with resolution III, with half of the experiments used in full factorial design [18]. The designs with resolution of III, IV and V are particularly important; they are listed in detail below [35].


 **Design resolution V:** These are the designs in which no main effect or second order interaction is associated with any other main effect or second order interaction, but the second order interactions are associated with third order interactions.

## **6. Sensitivity analysis development stages**

As a way of simplifying DOE application in sensitivity analysis of simulation models, the following sequence of steps is proposed.

The flowchart proposed in Figure 3 presents the necessary steps for conducting sensitivity analysis in discrete-event simulation models. Four stages are defined for the proposed method:

1. Simulation;

72 Discrete Event Simulations – Development and Applications

**5.3. Fractional factorial designs** 

effects and lower order interactions [35].

required for full factorial design [4].

are present [18].

design.

in detail below [35].

Factorial designs are more efficient than the one-at-a-time approach, as they allow for the factors' effects to be estimated via the other factors' levels, thus leading to valid conclusions about the experimental scope; they are also the only manner to discover interactions among the variables, thus avoiding erroneous conclusions when interactions between the factors

When there is little interest in interaction behavior among the factors which compose the

For example, consider a factorial design of 25. In this planning, five degrees of freedom correspond to the main effects, 10 degrees of freedom correspond to second order interactions and 16 degrees of freedom correspond to the highest order of interactions. In initial system or project studies, there is little interest in the highest level of interactions [35]. If interactions can be disregarded, fractional factorial design involving fewer executions for a complete set of 2k executions can be used in order to obtain information about the main

Thus, fractional factorial design provides a means by which to obtain estimates of main effects and, perhaps, second order interactions, with a fraction of the computational force

According to [18], the greatest application of fractional factorial designs is in screening experiments, where many factors are present in a system and the objective is to identify which factors indeed exercise a significant effect over the given response of interest. For the factors identified as significant through the use of fractional designs, the author recommends a more careful analysis with the use of other designs, such as full factorial

In fractional factorial design, a subset of 2k-p is constructed from a set of all of the possible

For this type of factorial design, the analyst must be attentive to its resolution. According to [35], the concept of resolution design is the form in which the fractional factorial designs are related in accordance with the associative standards which they produce. A design's

factorial design with resolution III, with half of the experiments used in full factorial design [18]. The designs with resolution of III, IV and V are particularly important; they are listed

 **Design resolution III:** These are the designs in which the main effect is associated with another main effect, but these main effects are associated with second order interactions

 **Design resolution IV:** These are the designs in which the main effect is associated with any other main effect or any other second order interaction, but second order

�����represents a

points for a 2k design, and a simulation is executed for only the chosen points [4].

resolution is represented by a Roman subscript number; for example, 2���

and second order interactions may be associated.

interactions are associated with each other.

system, this can be disregarded [35]. Instead, fractional factorial design can be used.


In the first step, the analyst should define the optimization objectives as well as verify if the simulation model is verified and validated. In doing so, the model will be apt to proceed to the following step, fractional factorial design. In this phase, the analyst should determine the model's input factors, their levels and select the response variables for analysis.

Once these initial steps have been completed, fractional factorial design can be applied. During execution of the experiments, the analyst should return to the simulation step and carry out the experiments in the simulation package. With the experiments done, the data should be analyzed using statistical means, determining the factors' significance level as well as its lower order interactions. At the end of this phase, the non-significant factors can be removed from analysis.

The third stage defined in the sensitivity analysis phase may be omitted by the analyst, depending on the degree of precision desired or for the cases in which the simulation models demand a lot of computational time in order to be processed. In this stage, a full factorial design is generated with experimental data, and only factors which show to be significant in the previous steps are tested. Depending on the necessity for executing more experiments in order to comply with full factorial design, the analyst will have to return to the simulation phase again to execute new experiments. In this stage, the residues should be analyzed in order to validate the results. Statistical analysis should be conducted once again in order to finalize this stage.

In the following stage, optimization simulation is utilized again. The significant factors after full factorial design application are utilized to configure the optimization tool, which is then executed. Many different configurations are tested for the input parameters until the optimizer converges on a solution. It is up to the analyst to evaluate the results and generate his or her conclusions and recommendations.

In order to demonstrate the utilization of this method in sensitivity analysis of discrete-event simulation models, two simulation models will be used as study objects in this chapter.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 75

The model in question was verified and validated, thus being ready for study. In order to verify the computational model and correct possible occurrences, some simulator resources such as counters, variables and signals were used, besides the conventional animation. Once the model was validated and verified, simulation was carried out. Statistical tests were utilized which compared the results obtained in the simulation model with data from the real system. The model was considered valid based on statistical tests which did not indicate

These conditions are indispensable for conducting sensitivity analysis. The utilization of a non-validated model would lead to erroneous conclusions and undesirable decision-making. For more information about validation and verification, readers are recommended to consult

[3]. Figure 4 presents an image of the model implemented in the simulation software.

**Figure 4.** Representation of the real system implemented in the simulation software

31 types of operations possible to be carried out depending on the type of product.

For the case in question, discrete variables were defined, with little variation in the lower and higher levels. This fact is justified due to the fact that the majority of simulation optimization problems work with such conditions; however, the experimentation can be conducted using other variable types and a greater variation between the upper and lower

For this study object, two levels were defined for each factor. For [18], when the experiment's objective is factor screening or process characterizations, it is common to utilize a small number of variables for each factor. According to the author, lower and higher levels are the most sufficient to obtain valid conclusions about the factors' significance. For example, even if the factor "Type 1 operators" exhibited the possibility of

The quality control station possesses the following characteristics:

 7 inspection posts; 3 operators;

19 types of products to be tested;

limits. Other types of applications can be seen in [43].

a statistical difference between real and simulated data.

**Figure 3.** Sequence of steps in order to conduct sensitivity analysis

It should be highlighted here that, in spite of the proposed modeling being applied to great variety of discrete-event simulation models, this approach could be unviable for models which demand a greater amount of computational time to be processed. In these cases, other types of experimental design which involve a smaller number of experiments could be used, such as Plackett-Burman Design; however, according to [18], such designs should be used with great care, as they possess limitations which should be evaluated carefully. An example of this shortcoming is the inability of certain designs to analyze main effects.

#### **7. Modeled systems**

The simulation models presented in this chapter were implemented in the software ProModel® and optimized in the package's optimization software, SimRunner®. However, it should be highlighted that the results presented could have been obtained using other commercial simulation packages. Likewise, a commercial statistical software package was utilized to analyze the data.

#### **7.1. Study object 1**

The first simulation model represents a quality control station from a telecommunications components factory. The cell is responsible for testing a diverse range of products before shipping them to final clients. This cell receives approximately 75% of the products from the six production lines in the company.

The model in question was verified and validated, thus being ready for study. In order to verify the computational model and correct possible occurrences, some simulator resources such as counters, variables and signals were used, besides the conventional animation. Once the model was validated and verified, simulation was carried out. Statistical tests were utilized which compared the results obtained in the simulation model with data from the real system. The model was considered valid based on statistical tests which did not indicate a statistical difference between real and simulated data.

These conditions are indispensable for conducting sensitivity analysis. The utilization of a non-validated model would lead to erroneous conclusions and undesirable decision-making. For more information about validation and verification, readers are recommended to consult [3]. Figure 4 presents an image of the model implemented in the simulation software.

**Figure 4.** Representation of the real system implemented in the simulation software

The quality control station possesses the following characteristics:


74 Discrete Event Simulations – Development and Applications

**Figure 3.** Sequence of steps in order to conduct sensitivity analysis

**7. Modeled systems** 

utilized to analyze the data.

six production lines in the company.

**7.1. Study object 1** 

It should be highlighted here that, in spite of the proposed modeling being applied to great variety of discrete-event simulation models, this approach could be unviable for models which demand a greater amount of computational time to be processed. In these cases, other types of experimental design which involve a smaller number of experiments could be used, such as Plackett-Burman Design; however, according to [18], such designs should be used with great care, as they possess limitations which should be evaluated carefully. An

example of this shortcoming is the inability of certain designs to analyze main effects.

The simulation models presented in this chapter were implemented in the software ProModel® and optimized in the package's optimization software, SimRunner®. However, it should be highlighted that the results presented could have been obtained using other commercial simulation packages. Likewise, a commercial statistical software package was

The first simulation model represents a quality control station from a telecommunications components factory. The cell is responsible for testing a diverse range of products before shipping them to final clients. This cell receives approximately 75% of the products from the


For the case in question, discrete variables were defined, with little variation in the lower and higher levels. This fact is justified due to the fact that the majority of simulation optimization problems work with such conditions; however, the experimentation can be conducted using other variable types and a greater variation between the upper and lower limits. Other types of applications can be seen in [43].

For this study object, two levels were defined for each factor. For [18], when the experiment's objective is factor screening or process characterizations, it is common to utilize a small number of variables for each factor. According to the author, lower and higher levels are the most sufficient to obtain valid conclusions about the factors' significance. For example, even if the factor "Type 1 operators" exhibited the possibility of


hiring 1 to 4 operators, for the purposes of the experimental matrix, only two levels would be considered: the lower level being 1 and the upper level being 4.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 77

The cell presents the following characteristics:

46 types of possible processes throughout the system.

**Table 3.** Experimental factors for the second study object

**8.1. Identification of significant factors** 

this situation, it is usually best to keep the number of factors levels low.

Variable Factor Lower level (-) Upper level (+) A Type 1 operators 1 2 B Type 2 operators 1 2 C Type 3 operators 1 2 D Type 1 machines 1 2 E Type 2 machines 1 2 F Type 3 machines 1 2 G Type 4 machines 1 2 H Type 5 machines 1 2 J Type 1 inspection posts 1 2 K Type 2 inspection posts 1 2 L Type 3 inspection posts 1 2 M Type 4 inspection posts 1 2

The optimum set of parameters is determined by three approaches similar to the first

According to [4], in simulation experimental designs provide a way to decide which specific configurations to simulate before the runs are performed so that the desired information can be obtained with the fewest simulation runs. For instance, considering the second application where there are 12 factors, if a full factorial experiment were chosen, 212 = 4096 runs would be necessary. Therefore, a screening experiment must be considered. Screening or characterization experiments are experiments in which many factors are considered and the objective is to identify those factors (if any) that have large effects [18]. Typically, screening experiment involves using fractional factorial designs and is performed in the early stages of the project when many factors are likely considered to have little or no effect on the response [18]. According to this author, in

For the first study object, ten experimental factors, each with two levels, were defined, as seen in Table 2. When considering full factorial design, a total number of 210 = 1024

8 different types of products;

 41 machines; 3 operators;

application.

**8. Experimentation** 

**8.2. Study object 1** 

**Table 2.** Experimental factors for the first study object

The optimum set of variables will be determined using three approaches. The first one performs several experiments to identify the main factors. After identifying the statistically significant simulation factors by using a two sample t hypothesis test (a usual procedure from any statistic package), the original fractional factorial design can -be converted to full factorial, eliminating the non-statistically significant terms. As these parameters are also important to the simulation arrangement, despite not being statistically significant, they can be kept constant in proper levels. Comparatively, a second approach can be established by using the main factors identified at the experiments DOE as input for the optimization via SimRunner. Finally, the third approach is performed using all ten factors in optimization via Simrunner.

### **7.2. Study object 2**

The second study object represents an automotive components production cell. The objective in this study object is to find the best combination of input variables which maximizes cell throughput. As with the previous case, the model was verified and validated, being ready for sensitivity analysis and optimization. Figure 5 shows an image of the model implemented in the simulation software.

**Figure 5.** Representation of the real system implemented in the simulation software

The cell presents the following characteristics:

41 machines;

76 Discrete Event Simulations – Development and Applications

**Table 2.** Experimental factors for the first study object

the model implemented in the simulation software.

**7.2. Study object 2** 

be considered: the lower level being 1 and the upper level being 4.

hiring 1 to 4 operators, for the purposes of the experimental matrix, only two levels would

Variable Factor Lower level (-) Upper level (+) A Type 1 operators 1 2 B Type 2 operators 1 2 C Type 3 operators 1 2 D Type 1 inspection posts 1 2 E Type 2 inspection posts 1 2 F Type 3 inspection posts 1 2 G Type 4 inspection posts 1 2 H Type 5 inspection posts 1 2 J Type 6 inspection posts 1 2 K Type 7 inspection posts 1 2

The optimum set of variables will be determined using three approaches. The first one performs several experiments to identify the main factors. After identifying the statistically significant simulation factors by using a two sample t hypothesis test (a usual procedure from any statistic package), the original fractional factorial design can -be converted to full factorial, eliminating the non-statistically significant terms. As these parameters are also important to the simulation arrangement, despite not being statistically significant, they can be kept constant in proper levels. Comparatively, a second approach can be established by using the main factors identified at the experiments DOE as input for the optimization via SimRunner. Finally, the third approach is performed using all ten factors in optimization via Simrunner.

The second study object represents an automotive components production cell. The objective in this study object is to find the best combination of input variables which maximizes cell throughput. As with the previous case, the model was verified and validated, being ready for sensitivity analysis and optimization. Figure 5 shows an image of

**Figure 5.** Representation of the real system implemented in the simulation software



**Table 3.** Experimental factors for the second study object

The optimum set of parameters is determined by three approaches similar to the first application.

## **8. Experimentation**

#### **8.1. Identification of significant factors**

According to [4], in simulation experimental designs provide a way to decide which specific configurations to simulate before the runs are performed so that the desired information can be obtained with the fewest simulation runs. For instance, considering the second application where there are 12 factors, if a full factorial experiment were chosen, 212 = 4096 runs would be necessary. Therefore, a screening experiment must be considered. Screening or characterization experiments are experiments in which many factors are considered and the objective is to identify those factors (if any) that have large effects [18]. Typically, screening experiment involves using fractional factorial designs and is performed in the early stages of the project when many factors are likely considered to have little or no effect on the response [18]. According to this author, in this situation, it is usually best to keep the number of factors levels low.

#### **8.2. Study object 1**

For the first study object, ten experimental factors, each with two levels, were defined, as seen in Table 2. When considering full factorial design, a total number of 210 = 1024 experiments would be necessary. In order to reduce the acceptable number of experiments, fractional factorial design is used.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 79

25 - - - + + - - - + - 93 26 + - - + + - - + - + 98 27 - + - + + - + - - + 99 28 + + - + + - + + + - 98 29 - - + + + - + + + + 95 30 + - + + + - + - - - 101 31 - + + + + - - + - - 98 32 + + + + + - - - + + 100 33 - - - - - + - - + + 99 34 + - - - - + - + - - 98 35 - + - - - + + - - - 97 36 + + - - - + + + + + 99 37 - - + - - + + + + - 98 38 + - + - - + + - - + 100 39 - + + - - + - + - + 98 40 + + + - - + - - + - 99 41 - - - + - + + + - + 98 42 + - - + - + + - + - 97 43 - + - + - + - + + - 100 44 + + - + - + - - - + 96 45 - - + + - + - - - - 99 46 + - + + - + - + + + 100 47 - + + + - + + - + + 98 48 + + + + - + + + - - 103 49 - - - - + + - - - - 99 50 + - - - + + - + + + 96 51 - + - - + + + - + + 101 52 + + - - + + + + - - 100 53 - - + - + + + + - + 99 54 + - + - + + + - + - 100 55 - + + - + + - + + - 98 56 + + + - + + - - - + 102 57 - - - + + + + + + - 96 58 + - - + + + + - - + 97 59 - + - + + + - + - + 95 60 + + - + + + - - + - 99 61 - - + + + + - - + + 97 62 + - + + + + - + - - 97 63 - + + + + + + - - - 99 64 + + + + + + + + + + 98

**Table 5.** The 2��

������ design matrix for principal fraction and results

Table 4 presents four factorial designs for 10 factors and their resolutions. As the objective of this analysis was to identify the model's sensitivity to certain factors, resolution IV was chosen. Resolution IV indicates that no main effect is associated with any other main effect or second order interaction, but there is interaction between certain second order interactions [18].


**Table 4.** Factorial designs for 10 factors and their resolutions




**Table 5.** The 2�� ������ design matrix for principal fraction and results

78 Discrete Event Simulations – Development and Applications

**Table 4.** Factorial designs for 10 factors and their resolutions

fractional factorial design is used.

interactions [18].

experiments would be necessary. In order to reduce the acceptable number of experiments,

Table 4 presents four factorial designs for 10 factors and their resolutions. As the objective of this analysis was to identify the model's sensitivity to certain factors, resolution IV was chosen. Resolution IV indicates that no main effect is associated with any other main effect or second order interaction, but there is interaction between certain second order

> **Fraction Resolution Design Executions** 1/8 V 2(10-3) 128 1/16 IV 2(10-4) 64 1/32 IV 2(10-5) 32 1/64 III 2(10-6) 16

**Experiment A B C D E F G H J K WIP**  1 - - - - - - + + + + 99 2 + - - - - - + - - - 95 3 - + - - - - - + - - 101 4 + + - - - - - - + + 104 5 - - + - - - - - + - 94 6 + - + - - - - + - + 98 7 - + + - - - + - - + 100 8 + + + - - - + + + - 99 9 - - - + - - - - - + 93 10 + - - + - - - + + - 95 11 - + - + - - + - + - 97 12 + + - + - - + + - + 98 13 - - + + - - + + - - 101 14 + - + + - - + - + + 101 15 - + + + - - - + + + 99 16 + + + + - - - - - - 98 17 - - - - + - + + - - 101 18 + - - - + - + - + + 95 19 - + - - + - - + + + 100 20 + + - - + - - - - - 98 21 - - + - + - - - - + 94 22 + - + - + - - + + - 95 23 - + + - + - + - + - 93 24 + + + - + - + + - + 99

Among the resolution IV designs presented in Table 4, the fractional factorial design 2�� ������ was chosen. This design, despite possessing a greater number of executions than the design 2�� ������, allows for the reduction of full factorial designs for six or less factors, without requiring new experiments; that is, the results of this design can be used in full factorial designs with six or less factors, without needing to conduct additional experiments. However, if preliminary studies show that there are a significant number of fewer variables (5 or less), the factorial design 2�� ������ could be chosen with no problems.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 81

1,267

������ fractional design with significance

With the aid of statistical software, it was possible to perform quantitative analysis of the

������ fractional design with a significance level of α = 5%

Lenth's PSE = 0,609375

By analyzing the figure, it can be seen that factor B (number of type 2 operators) and interaction CD (number of type 3 and type 1 inspection posts) are significant. According to [18], if the experimenter can reasonably assume that certain high-order interactions are negligible, information on the main effects and low-order interactions may be obtained. Otherwise, when there are several variables, the system or process is likely to be driven primarily by some of the main effects and low-order interactions. For this reason, it is reasonable to admit that factors A, E, F, G, H, J and K are not significant, although they are

0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6

**Effect**

Factors C and D may be considered significant, seeing as the interaction between the two factors are quite significant. In Figure 7, it can be seen that B and C exercise a positive effect over the WIP; shifting from the lower to upper level causes an increase in the WIP. Inversely, factor D exercises a negative effect; shifting from the lower to upper level causes WIP to fall. Analysis of interaction behavior in fractional factorial design is not recommended, since the effects possess the property of aliases. That is, according to [18], two or more factors are aliased when it is not possible to distinguish between the effects of overlapping factors. Only three main factors may be considered significant (B, C, D), and full factorial with these factors can be carried out with the data from the 64 experiments.

Factorial design's own structure helps explain why only factors B, C and D were chosen to compose the new factorial design. The main reason B is aliased with the triple interaction

still necessary for the simulation model and must be kept at the lower level (-).

stored data. Figure 6 presents the Pareto Chart for 2��

level of 5%.

**Figure 6.** Pareto chart for 2��

CG DG CF HJ BEF J BJ ACG CE AFG C BF G BC FK E EH F A FH EJ BCJ CJ HK GJ AGK AH AC B CD

**Term**

It is worth mentioning that factorials with a resolution less than IV should be omitted because, in these types of design, the main effects are associated with second level interactions, and second level interactions could possess interaction between each other, as well, thus making these designs undesirable.

Table 5 shows the design matrix for the principal fraction and the results obtained for the WIP. Wip represents the total number of pieces in quality control inspection; its result is shown by the variable introduced in the simulation model which subtracts the pieces which leave the system (inspected pieces) from the total number of entities which enter the system (pieces to be inspected). This value is then offered at the conclusion of the simulation in which the report is generated. In the table, the best results attained with the experimentation are shown. In Table 5, the symbols – and + indicate the lower and upper levels shown in Table 2, respectively.

As an example, the number of operator types 1, 2 and 3 and the number of inspection posts types 1, 2 and 3 (A B C D E F) were defined in the simulator as being equal to the lower level (1); for the inspection posts types 4, 5, 6 and 7 (G H J K), the upper level was defined. A replica utilizing this configuration was run in the simulation software, and work in process (WIP) statistics were stored for analysis. This process was repeated 63 other times until all of the experimental matrix's configurations were run.

The 2�� ������ fractional factorial design used in this research was not replicated. Therefore, it is not possible to assess the significance of the main and interaction effects using the conventional bilateral t-test or ANOVA. The standard analysis procedure for a nonreplicated two-level design is a normal plot of the estimated factor's effects. However, these designs are so widely used in practice that many formal analysis procedures have been proposed to overcome the subjectivity of normal probability [18]. [44], for instance, recommend the use of Lenth's method, a graphical approach based on a Pareto Chart for the error term. If the error term has one or more degrees of freedom, the line on the graph is drawn at t, where t is the (1 - α/2) quartile of a t-distribution with a number of degrees of freedom equal to the number of effects/3. The vertical line in the Pareto Chart is the margin of error, defined as ME = t x PSE. Lenth's pseudo standard error (PSE) is based on sparseness of the effects principle, which assumes the variation in the smallest effects is due to the random error. To calculate PSE the following steps are necessary: (a) calculates the absolute value of the effects; (b) calculates S, which is 1.5 x median of the step (c); calculates the median of the effects that are less than 2.5 x S and (d) calculates PSE, which is 1.5 x median calculated in step (c).

With the aid of statistical software, it was possible to perform quantitative analysis of the stored data. Figure 6 presents the Pareto Chart for 2�� ������ fractional design with significance level of 5%.

80 Discrete Event Simulations – Development and Applications

well, thus making these designs undesirable.

the experimental matrix's configurations were run.

(5 or less), the factorial design 2��

Table 2, respectively.

median calculated in step (c).

The 2��

2��

Among the resolution IV designs presented in Table 4, the fractional factorial design 2��

was chosen. This design, despite possessing a greater number of executions than the design

������, allows for the reduction of full factorial designs for six or less factors, without requiring new experiments; that is, the results of this design can be used in full factorial designs with six or less factors, without needing to conduct additional experiments. However, if preliminary studies show that there are a significant number of fewer variables

It is worth mentioning that factorials with a resolution less than IV should be omitted because, in these types of design, the main effects are associated with second level interactions, and second level interactions could possess interaction between each other, as

Table 5 shows the design matrix for the principal fraction and the results obtained for the WIP. Wip represents the total number of pieces in quality control inspection; its result is shown by the variable introduced in the simulation model which subtracts the pieces which leave the system (inspected pieces) from the total number of entities which enter the system (pieces to be inspected). This value is then offered at the conclusion of the simulation in which the report is generated. In the table, the best results attained with the experimentation are shown. In Table 5, the symbols – and + indicate the lower and upper levels shown in

As an example, the number of operator types 1, 2 and 3 and the number of inspection posts types 1, 2 and 3 (A B C D E F) were defined in the simulator as being equal to the lower level (1); for the inspection posts types 4, 5, 6 and 7 (G H J K), the upper level was defined. A replica utilizing this configuration was run in the simulation software, and work in process (WIP) statistics were stored for analysis. This process was repeated 63 other times until all of

������ fractional factorial design used in this research was not replicated. Therefore, it is not possible to assess the significance of the main and interaction effects using the conventional bilateral t-test or ANOVA. The standard analysis procedure for a nonreplicated two-level design is a normal plot of the estimated factor's effects. However, these designs are so widely used in practice that many formal analysis procedures have been proposed to overcome the subjectivity of normal probability [18]. [44], for instance, recommend the use of Lenth's method, a graphical approach based on a Pareto Chart for the error term. If the error term has one or more degrees of freedom, the line on the graph is drawn at t, where t is the (1 - α/2) quartile of a t-distribution with a number of degrees of freedom equal to the number of effects/3. The vertical line in the Pareto Chart is the margin of error, defined as ME = t x PSE. Lenth's pseudo standard error (PSE) is based on sparseness of the effects principle, which assumes the variation in the smallest effects is due to the random error. To calculate PSE the following steps are necessary: (a) calculates the absolute value of the effects; (b) calculates S, which is 1.5 x median of the step (c); calculates the median of the effects that are less than 2.5 x S and (d) calculates PSE, which is 1.5 x

������ could be chosen with no problems.

������

**Figure 6.** Pareto chart for 2�� ������ fractional design with a significance level of α = 5%

By analyzing the figure, it can be seen that factor B (number of type 2 operators) and interaction CD (number of type 3 and type 1 inspection posts) are significant. According to [18], if the experimenter can reasonably assume that certain high-order interactions are negligible, information on the main effects and low-order interactions may be obtained. Otherwise, when there are several variables, the system or process is likely to be driven primarily by some of the main effects and low-order interactions. For this reason, it is reasonable to admit that factors A, E, F, G, H, J and K are not significant, although they are still necessary for the simulation model and must be kept at the lower level (-).

Factors C and D may be considered significant, seeing as the interaction between the two factors are quite significant. In Figure 7, it can be seen that B and C exercise a positive effect over the WIP; shifting from the lower to upper level causes an increase in the WIP. Inversely, factor D exercises a negative effect; shifting from the lower to upper level causes WIP to fall. Analysis of interaction behavior in fractional factorial design is not recommended, since the effects possess the property of aliases. That is, according to [18], two or more factors are aliased when it is not possible to distinguish between the effects of overlapping factors. Only three main factors may be considered significant (B, C, D), and full factorial with these factors can be carried out with the data from the 64 experiments.

Factorial design's own structure helps explain why only factors B, C and D were chosen to compose the new factorial design. The main reason B is aliased with the triple interaction

AGH, which can be disregarded according to the chosen resolution for factorial design and by the sparsity of effects principle [18]. In turn, interaction CD is aliased with two triple interactions, AFH and BFG, which can be disregarded, just as in the last case. It can also be aliased with the double interaction in JK; however, as the main factors J and K are not significant, this interaction may also be discarded. The alias structure used in this analysis is available in many statistical packages.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 83

Mean 0 StDev 2,044 N 64 AD 0,433 P-Value 0,295


1 5 10 15 20 25 30 35 40 45 50 55 60

Once the residual validity was verified, the results could be analyzed using DOE. The analyses continued to be carried out via graphical analysis due to its ease of

Number of runs about median: 36 Expected number of runs: 32,9 Longest run about median: 6 Approx P-Value for Clustering: 0,785 Approx P-Value for Mixtures: 0,215

**Figure 9.** Verification of the residues independence

**Observation**

Number of runs up or down: 43 Expected number of runs: 42,3 Longest run up or down: 4 Approx P-Value for Trends: 0,579 Approx P-Value for Oscillation: 0,421

Normal

**RESIDUALS**

99,9

**Percent**

0,1

5,0

2,5

0,0

**RESIDUALS**

comprehension.



**Figure 8.** Verification of the residues normality

Fractional factorial design 2�� ������ was converted into a full factorial design 23 with replicas. Residual analysis is also possible, since now there are replicas, seeing that the experimentation changed from fractional factorial design 2�� ������to full factorial design 23 (8 experiments).

**Figure 7.** Main effects plot for WIP

According to [18], the residues need to be normal, random and not correlated in order to validate the experimental values obtained. Figure 8 shows the verification of the residues normality.

Evaluating the normality probability graph, one can see that the data are adjusted to a normal distribution, as is evidence by the way the points fall over the line in the graph as well as analysis of the P-value. One can see the data points follow the straight line, and the P-value for the normality test was less than 0.05, leading to the conclusion that the data are normally distributed. Figure 9 shows the verification of the residues independence. The standardized graphs versus the observed values do not present any random patterns of grouping or bias.

**Figure 8.** Verification of the residues normality

available in many statistical packages.


A


J K



According to [18], the residues need to be normal, random and not correlated in order to validate the experimental values obtained. Figure 8 shows the verification of the residues

Evaluating the normality probability graph, one can see that the data are adjusted to a normal distribution, as is evidence by the way the points fall over the line in the graph as well as analysis of the P-value. One can see the data points follow the straight line, and the P-value for the normality test was less than 0.05, leading to the conclusion that the data are normally distributed. Figure 9 shows the verification of the residues independence. The standardized graphs versus the observed values do not present any random patterns of

experimentation changed from fractional factorial design 2��

Fractional factorial design 2��

98,5 98,0 97,5

98,5 98,0 97,5

**Mean**

98,5 98,0 97,5

**Figure 7.** Main effects plot for WIP

normality.

grouping or bias.

experiments).

AGH, which can be disregarded according to the chosen resolution for factorial design and by the sparsity of effects principle [18]. In turn, interaction CD is aliased with two triple interactions, AFH and BFG, which can be disregarded, just as in the last case. It can also be aliased with the double interaction in JK; however, as the main factors J and K are not significant, this interaction may also be discarded. The alias structure used in this analysis is

Residual analysis is also possible, since now there are replicas, seeing that the

������ was converted into a full factorial design 23 with replicas.


E F G H

B C D


������to full factorial design 23 (8

**Figure 9.** Verification of the residues independence

Once the residual validity was verified, the results could be analyzed using DOE. The analyses continued to be carried out via graphical analysis due to its ease of comprehension.

Figures 10 and 11 present the analysis for the new design. By analyzing figure 10, it can be verified that factor B (number of type 2 operators) and the interaction CD (number of type 3 operators and number of type 1 inspection posts) remained significant. In this new design, no other main factor or interaction demonstrated a significance level greater than 5%.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 85

**C**

**B**

For the second study object, 12 experimental factors were defined, as presented in Table 3. Each factor possesses two levels. Unlike the previous case, for this study the objective is to maximize the manufacturing cell's throughput. Thus, the significance of the 12 factors will be analyzed. Considering full factorial design, 212 = 4096 experiments would be necessary. Similar to the previous case, fractional factorial design was used to reduce the number of

Table 6 presents four factorial designs for 12 factors and their resolutions. As the analysis objective in this case is to identify the model's performance sensitivity to the factors, a

> **Fraction Resolution Design Executions**  1/256 III 2(12-8) 16 1/128 IV 2(12-7) 32 1/64 IV 2(12-6) 64 1/32 IV 2(12-5) 128

Out of the resolution IV plans presented in Table 6, fractional factorial design 2��

chosen. Researchers opted for this design due to its location between the fractional factorial

������ , thus enabling a reduction to full factorial design for six or less


**Interaction Plot**  Data Means

99




������ was

98

97 99

98

97

**D**


D

B

**Main Effects Plot**  Data Means


resolution IV design was chosen.

������ and 2��

**Figure 11.** Factorial and interaction plots for 23 full factorial design

**Table 6.** Factorial designs for 10 factors and their resolutions


C

99,0

98,5

98,0

97,5

**Mean**

99,0

98,5

98,0

97,5

**8.3. Study object 2** 

experiments.

designs 2��

**Figure 10.** Pareto Chart for full factorial design with significance level α = 5%

Analysis of Figure 11 shows that factors B and C exercise a positive effect over the WIP; that is, they should be kept at the lower level in order to minimize the WIP count. Factor D should be kept at the upper level, since it exercises a negative effect on the WIP count. By observing Figure 11, it can be seen that the CD interaction has a strong effect on diminishing the WIP when the main effects C and D remain at their own respective lower and upper levels.

There are strong indications about an improved configuration for the input variables in order to minimize the WIP count. However, these suppositions will be tested using commercial optimization software. First, 10 input variables will be utilized; afterwards, only the three variables which showed to be significant according to the study in this section will be evaluated, and the seven other factors will be fixed in the lower level, seeing as they are not significant.

**Figure 11.** Factorial and interaction plots for 23 full factorial design

#### **8.3. Study object 2**

84 Discrete Event Simulations – Development and Applications

than 5%.

BCD

BD

D

C

**Term**

levels.

not significant.

BC

B

CD

Figures 10 and 11 present the analysis for the new design. By analyzing figure 10, it can be verified that factor B (number of type 2 operators) and the interaction CD (number of type 3 operators and number of type 1 inspection posts) remained significant. In this new design, no other main factor or interaction demonstrated a significance level greater

2,003

0,0 0,5 1,0 1,5 2,0 2,5 3,0

Analysis of Figure 11 shows that factors B and C exercise a positive effect over the WIP; that is, they should be kept at the lower level in order to minimize the WIP count. Factor D should be kept at the upper level, since it exercises a negative effect on the WIP count. By observing Figure 11, it can be seen that the CD interaction has a strong effect on diminishing the WIP when the main effects C and D remain at their own respective lower and upper

There are strong indications about an improved configuration for the input variables in order to minimize the WIP count. However, these suppositions will be tested using commercial optimization software. First, 10 input variables will be utilized; afterwards, only the three variables which showed to be significant according to the study in this section will be evaluated, and the seven other factors will be fixed in the lower level, seeing as they are

**Figure 10.** Pareto Chart for full factorial design with significance level α = 5%

**Standardized Effect**

For the second study object, 12 experimental factors were defined, as presented in Table 3. Each factor possesses two levels. Unlike the previous case, for this study the objective is to maximize the manufacturing cell's throughput. Thus, the significance of the 12 factors will be analyzed. Considering full factorial design, 212 = 4096 experiments would be necessary. Similar to the previous case, fractional factorial design was used to reduce the number of experiments.

Table 6 presents four factorial designs for 12 factors and their resolutions. As the analysis objective in this case is to identify the model's performance sensitivity to the factors, a resolution IV design was chosen.


**Table 6.** Factorial designs for 10 factors and their resolutions

Out of the resolution IV plans presented in Table 6, fractional factorial design 2�� ������ was chosen. Researchers opted for this design due to its location between the fractional factorial designs 2�� ������ and 2�� ������ , thus enabling a reduction to full factorial design for six or less factors, without having to carry out new experiments. It should be noted here that if more than six factors are shown to be significant, another factorial design could be done while taking advantage of the data already acquired from the fractional factorial design 2�� ������ and realizing only the non-tested experiments.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 87

������ was utilized without replicas.

33 - - - - - + + - + - - - 400400 34 + - - - - + + + + - + + 423800 35 - + - - - + + + - + + - 421200 36 + + - - - + + - - + - + 421200 37 - - + - - + + + - + - + 400400 38 + - + - - + + - - + + - 397800 39 - + + - - + + - + - + + 405600 40 + + + - - + + + + - - - 444600 41 - - - + - + - - - + - - 371800 42 + - - + - + - + - + + + 410800 43 - + - + - + - + + - + - 403000 44 + + - + - + - - + - - + 382200 45 - - + + - + - + + - - + 384800 46 + - + + - + - - + - + - 387400 47 - + + + - + - - - + + + 395200 48 + + + + - + - + - + - - 421200 49 - - - - + + - - - - + + 397800 50 + - - - + + - + - - - - 423800 51 - + - - + + - + + + - + 382200 52 + + - - + + - - + + + - 413400 53 - - + - + + - + + + + - 395200 54 + - + - + + - - + + - + 405600 55 - + + - + + - - - - - - 392600 56 + + + - + + - + - - + + 416000 57 - - - + + + + - + + + + 395200 58 + - - + + + + + + + - - 429000 59 - + - + + + + + - - - + 429000 60 + + - + + + + - - - + - 434200 61 - - + + + + + + - - + - 410800 62 + - + + + + + - - - - + 423800 63 - + + + + + + - + + - - 418600 64 + + + + + + + + + + + + 449800

**Table 7.** The 2��

Chart for 2��

������ design matrix for principal fraction and results

With the help of statistical software, the data were analyzed. Figure 12 shows the Pareto

first figure, it can be seen that factors G (number of type 4 machines), A (number of type 1 operators), H (number of type 5 machines), B (number of type 2 operators), E (number of type 2 machines) and the interaction BG (number of type two operators and number of type 4 machines) are significant according to the adopted significance level. It can be said that the factors C, D, F, J, K, L and M are not significant, although they are necessary for the simulation and their values may be fixed at the lower level

������, fractional design with a significance level of 5%. By analyzing the

Similar to the last case, fractional factorial design 2��

Table 7 presents the experimental design matrix for the principal fraction and throughput. Throughput represents the number of pieces produced by the manufacturing cell. For this case, researchers needed to create a variable to store the number of pieces produced and present this value at the end of the simulation. The greatest value produced is highlighted.



**Table 7.** The 2�� ������ design matrix for principal fraction and results

86 Discrete Event Simulations – Development and Applications

and realizing only the non-tested experiments.

factors, without having to carry out new experiments. It should be noted here that if more than six factors are shown to be significant, another factorial design could be done while taking advantage of the data already acquired from the fractional factorial design 2��

Table 7 presents the experimental design matrix for the principal fraction and throughput. Throughput represents the number of pieces produced by the manufacturing cell. For this case, researchers needed to create a variable to store the number of pieces produced and present this value at the end of the simulation. The greatest value produced is highlighted.

**Experiment A B C D E F G H J K L M Throughput**  1 - - - - - - - - + + + + 369200 2 + - - - - - - + + + - - 408200 3 - + - - - - - + - - - + 392600 4 + + - - - - - - - - + - 390000 5 - - + - - - - + - - + - 392600 6 + - + - - - - - - - - + 390000 7 - + + - - - - - + + - - 400400 8 + + + - - - - + + + + + 416000 9 - - - + - - + - - - + + 403000 10 + - - + - - + + - - - - 413400 11 - + - + - - + + + + - + 429000 12 + + - + - - + - + + + - 429000 13 - - + + - - + + + + + - 413400 14 + - + + - - + - + + - + 408200 15 - + + + - - + - - - - - 421200 16 + + + + - - + + - - + + 429000 17 - - - - + - + - - + - - 390000 18 + - - - + - + + - + + + 429000 19 - + - - + - + + + - + - 426400 20 + + - - + - + - + - - + 431600 21 - - + - + - + + + - - + 416000 22 + - + - + - + - + - + - 410800 23 - + + - + - + - - + + + 410800 24 + + + - + - + + - + - - 434200 25 - - - + + - - - + - - - 392600 26 + - - + + - - + + - + + 408200 27 - + - + + - - + - + + - 400400 28 + + - + + - - - - + - + 395200 29 - - + + + - - + - + - + 405600 30 + - + + + - - - - + + - 397800 31 - + + + + - - - + - + + 395200 32 + + + + + - - + + - - - 429000

������

Similar to the last case, fractional factorial design 2�� ������ was utilized without replicas. With the help of statistical software, the data were analyzed. Figure 12 shows the Pareto Chart for 2�� ������, fractional design with a significance level of 5%. By analyzing the first figure, it can be seen that factors G (number of type 4 machines), A (number of type 1 operators), H (number of type 5 machines), B (number of type 2 operators), E (number of type 2 machines) and the interaction BG (number of type two operators and number of type 4 machines) are significant according to the adopted significance level. It can be said that the factors C, D, F, J, K, L and M are not significant, although they are necessary for the simulation and their values may be fixed at the lower level

(assuming a value of 1), since they do not exercise a significant effect over the model's throughput.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 89


E F G H

J K L M

Before analyzing the new design's results, the validity of the results was tested, as was the case with the previous experiment. Once the validity of the new design's residues' was

With the new design, the main factors A, B, E, G and H, the interactions BG, AGH, ABEG and AH showed to be significant, according to the 5% significance level (Figure 14). All of the main factors presented positive effects on the throughput, according to Figure 15; that is, shifting from the lower (–) and upper (+) level, the production cell's throughput increases.

It can then be concluded that, although the simulation model possesses 12 input variables which can be arranged in order to maximize throughput, only five variables significantly contribute to increased production. In following, a comparison will be performed between the commercial optimization software; first, all 12 input variables will be optimized, and

It is worth mentioning here another optimization approach which is commonly utilized in simulation optimization. By using full factorial design with replicas, it is possible to generate a metamodel for the response variable under analysis. With a mathematical model on hand, traditional optimization tools such as Microsoft Excel's Solver may be utilized in place of

Another approach which is commonly employed in the literature is Response Surface Methodology. As was cited in the previous strategy, a mathematical model, generally nonlinear of second-order, is generated through the experimental data and then is optimized. The shortcoming of this strategy is that the model must possess a robust fit, which allows

verified, it was then possible to statistically analyze the results with DOE.

then only the five variables which are statistically significant will be optimized.

simulation optimization tools. An example of such a technique can be seen in [43].

B C D



������ was converted into a full factorial design 25 with replicas.


A



**Figure 13.** Main effects plot for throughput

420000

410000

400000

420000

410000

**Mean**

400000

420000

410000

400000

Fractional factorial design 2��

By analyzing Figure 13, it can be seen that the main factors A, B, E, G, and G exercise a positive effect over throughput; that is, by altering the lower level to the upper level, there is an increase in throughput.

Interaction behavior analysis in fractional factorial design is not recommended, seeing that aliasing between effects tends to emerge. Thus, full factorial design with the significant factors will be carried out using the data from the 64 experiments carried out. As was the previous case, the alias structure for factorial design 2�� ������enables the explanation for why A, B, E, G and H were chosen to make up the full factorial design.

The main factor A is aliased with the two triple interactions BCH and HLM. Factor B is aliased with two other triple interactions, ACH and CLM. Factor E is aliased with the interactions DFG and FJK. Factor G is aliased with DEF and DKJ. Factor H is aliased with the interactions ABC and ALM. Finally, the interaction BG is aliased with four triple interactions, ADL, CEK, CFJ and DHM. All these interactions can be disregarded according to the chosen level of resolution for factorial design and the principle of sparsity of effects principle [18].

It can be concluded that, although the simulation model possesses 12 input variables which may be arranged in order to maximize total throughput, only 5 variables significantly contribute to increased throughput. In following, an optimization of this simulation model will be executed, using a commercial software package, first optimizing all 12 input variables, and then with only the 5 input variables which are statistically significant.

**Figure 12.** Pareto Chart for 2�� ������ fractional design with significance level of 5%

**Figure 13.** Main effects plot for throughput

previous case, the alias structure for factorial design 2��

5069

A, B, E, G and H were chosen to make up the full factorial design.

throughput.

an increase in throughput.

CD KL BD DF FL BJ M D BK C DG CK GM DK AD AK CG BE AE AF AH EK JM BM BG E B H A G

**Figure 12.** Pareto Chart for 2��

**Term**

(assuming a value of 1), since they do not exercise a significant effect over the model's

By analyzing Figure 13, it can be seen that the main factors A, B, E, G, and G exercise a positive effect over throughput; that is, by altering the lower level to the upper level, there is

Interaction behavior analysis in fractional factorial design is not recommended, seeing that aliasing between effects tends to emerge. Thus, full factorial design with the significant factors will be carried out using the data from the 64 experiments carried out. As was the

The main factor A is aliased with the two triple interactions BCH and HLM. Factor B is aliased with two other triple interactions, ACH and CLM. Factor E is aliased with the interactions DFG and FJK. Factor G is aliased with DEF and DKJ. Factor H is aliased with the interactions ABC and ALM. Finally, the interaction BG is aliased with four triple interactions, ADL, CEK, CFJ and DHM. All these interactions can be disregarded according to the chosen level of resolution for factorial design and the principle of sparsity of effects principle [18].

It can be concluded that, although the simulation model possesses 12 input variables which may be arranged in order to maximize total throughput, only 5 variables significantly contribute to increased throughput. In following, an optimization of this simulation model will be executed, using a commercial software package, first optimizing all 12 input

0 5000 10000 15000 20000

**Effect**

������ fractional design with significance level of 5%

variables, and then with only the 5 input variables which are statistically significant.

������enables the explanation for why

Lenth's PSE = 2437,5

Fractional factorial design 2�� ������ was converted into a full factorial design 25 with replicas. Before analyzing the new design's results, the validity of the results was tested, as was the case with the previous experiment. Once the validity of the new design's residues' was verified, it was then possible to statistically analyze the results with DOE.

With the new design, the main factors A, B, E, G and H, the interactions BG, AGH, ABEG and AH showed to be significant, according to the 5% significance level (Figure 14). All of the main factors presented positive effects on the throughput, according to Figure 15; that is, shifting from the lower (–) and upper (+) level, the production cell's throughput increases.

It can then be concluded that, although the simulation model possesses 12 input variables which can be arranged in order to maximize throughput, only five variables significantly contribute to increased production. In following, a comparison will be performed between the commercial optimization software; first, all 12 input variables will be optimized, and then only the five variables which are statistically significant will be optimized.

It is worth mentioning here another optimization approach which is commonly utilized in simulation optimization. By using full factorial design with replicas, it is possible to generate a metamodel for the response variable under analysis. With a mathematical model on hand, traditional optimization tools such as Microsoft Excel's Solver may be utilized in place of simulation optimization tools. An example of such a technique can be seen in [43].

Another approach which is commonly employed in the literature is Response Surface Methodology. As was cited in the previous strategy, a mathematical model, generally nonlinear of second-order, is generated through the experimental data and then is optimized. The shortcoming of this strategy is that the model must possess a robust fit, which allows for a satisfactory representation of the response factor. If the model is not robust, experimental strategies should be employed in order to redefine the experimental region, which many times is not applicable for simulation optimization problems.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 91

Through the sensitivity analysis, each simulation model's significant variables may be identified. Considering only these results, the best combination of input variables can be inferred in order to optimize the simulation models; however, there is no way of

One way of confirming these results is through optimization. An example application of simulation optimization will be employed as a means of evaluating the efficiency of

The adopted procedure is to optimize the study objects in two different ways. In the first case, all input variables will be optimized and, in the second case, only the factors selected in the sensitivity analysis will be optimized. Finally, the results attained will be compared in order to verify if the design techniques were advantageous to the process. The time involved in the process will not be the basis for comparison; thus it is evident that the number of

The simulation software package SimRunner® from the ProModel Corporation will be utilized for the execution of experiments; however, there are other simulation optimization software packages that could have been chosen for this investigation. SimRunner® integrates resources to analyze and optimize simulation models through multivariable optimization. This type of optimization tests multiple factor combinations in search of the system input variable configuration which leads to the best objective function value [20].

SimRunner® is based in a genetic algorithm and possesses three optimization profiles: Aggressive, Moderate and Cautious. These profiles are directly related to the confidence in the solution and the time necessary to find this solution. The cautious profile was chosen for this study in order to consider the greatest possible number of solutions and in turn,

The optimization objective for the first simulation model was to find the best input variable combination in order to minimize the system's work in process count. As presented in Table 2, this model possesses 10 input variables, being varied at the lower level (1) and the upper level (2). In the first optimization stage, 10 input model variables were selected and the optimizer was

The optimizer converged with 296 experiments. The best result obtained was 92, which was attained during experiment 261, as seen in Figure 16. The values found for the factors are

The sensitivity analysis for the first study object identified three factors with significant effects which can be utilized for simulation (Table 9). The other variables were maintained at

**9. Simulation model optimization** 

guaranteeing this affirmation based in only a sensitivity analysis.

fractional factorial design in the execution of sensitivity analysis.

experiments necessary will be reduced for the model to arrive at a solution.

guarantee a more comprehensive search and present better responses [45].

**9.1. Optimization of the first study object** 

configured. The results found can be seen in Figure 16.

their original values, defined as the lower level (\*).

shown in Table 8.

**Figure 14.** Pareto Chart for 25 full factorial design with significance level α = 5%

**Figure 15.** Factorial and interaction plots for 25 full factorial design

## **9. Simulation model optimization**

90 Discrete Event Simulations – Development and Applications

2,04

BGH ABEH AB GH BEH ABG BEGH AEGH BH AEG EH EG AEH ABE ABEGH ABH BEG EGH ABGH BE AE AH ABEG AGH BG E B H A G

**Mean**

1 2

G H

**Main Effects Plot for prod1** Data Means

A

1 2

**Term**

for a satisfactory representation of the response factor. If the model is not robust, experimental strategies should be employed in order to redefine the experimental region,

0 2 4 6 8 10 12

**Standardized Effect**

**A**

**B**

**E**

**G**

**H**

1 2 1 2 1 2 1 2

**Interaction Plot for prod1** Data Means

1 2 A

1 2 A

1 2 A

> 1 2 G

B

which many times is not applicable for simulation optimization problems.

**Figure 14.** Pareto Chart for 25 full factorial design with significance level α = 5%

1 2 1 2

B E

1 2

**Figure 15.** Factorial and interaction plots for 25 full factorial design

Through the sensitivity analysis, each simulation model's significant variables may be identified. Considering only these results, the best combination of input variables can be inferred in order to optimize the simulation models; however, there is no way of guaranteeing this affirmation based in only a sensitivity analysis.

One way of confirming these results is through optimization. An example application of simulation optimization will be employed as a means of evaluating the efficiency of fractional factorial design in the execution of sensitivity analysis.

The adopted procedure is to optimize the study objects in two different ways. In the first case, all input variables will be optimized and, in the second case, only the factors selected in the sensitivity analysis will be optimized. Finally, the results attained will be compared in order to verify if the design techniques were advantageous to the process. The time involved in the process will not be the basis for comparison; thus it is evident that the number of experiments necessary will be reduced for the model to arrive at a solution.

The simulation software package SimRunner® from the ProModel Corporation will be utilized for the execution of experiments; however, there are other simulation optimization software packages that could have been chosen for this investigation. SimRunner® integrates resources to analyze and optimize simulation models through multivariable optimization. This type of optimization tests multiple factor combinations in search of the system input variable configuration which leads to the best objective function value [20].

SimRunner® is based in a genetic algorithm and possesses three optimization profiles: Aggressive, Moderate and Cautious. These profiles are directly related to the confidence in the solution and the time necessary to find this solution. The cautious profile was chosen for this study in order to consider the greatest possible number of solutions and in turn, guarantee a more comprehensive search and present better responses [45].

#### **9.1. Optimization of the first study object**

The optimization objective for the first simulation model was to find the best input variable combination in order to minimize the system's work in process count. As presented in Table 2, this model possesses 10 input variables, being varied at the lower level (1) and the upper level (2).

In the first optimization stage, 10 input model variables were selected and the optimizer was configured. The results found can be seen in Figure 16.

The optimizer converged with 296 experiments. The best result obtained was 92, which was attained during experiment 261, as seen in Figure 16. The values found for the factors are shown in Table 8.

The sensitivity analysis for the first study object identified three factors with significant effects which can be utilized for simulation (Table 9). The other variables were maintained at their original values, defined as the lower level (\*).

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 93

**Figure 17.** Performance measures plot for optimization using significant factors

**Table 10.** The best solution for optimization using significant factors

**9.2. Optimization of the second study object** 

varied from the lower level (1) to the upper level (2).

Figure 18.

**Factor Variable Value** A Type 1 operators 1\* B Type 2 operators 1 C Type 3 operators 1 D Type 1 inspection posts 2 E Type 2 inspection posts 1\* F Type 3 inspection posts 1\* G Type 4 inspection posts 1\* H Type 5 inspection posts 1\* J Type 6 inspection posts 1\* K Type 7 inspection posts 1\* Note: The parameters identified with (\*) were not used as input for optimization. They were kept at lower level (1).

The optimization objective for the second simulation model was to find the best combination of model input variables in order to maximize the manufacturing cell's throughput. As presented in Table 3, the model possesses 12 input variables which are

In the first optimization phase, the 12 model input variables were selected and the optimization software was set up for experimentation. The results found can be seen in

The optimizer converged with 173 experiments. The best result found was 452,400, which was obtained in experiment 10 (Figure 18). The obtained values are shown in Table 11.

**Figure 16.** Performance measures plot for optimization using all factors


**Table 8.** The best solution for optimization using all factors


**Table 9.** Significant factors for the first study object

The results found can be seen in Figure 17.

Simrunner® converged after 8 experiments. The best result was 93, which was obtained in the seventh experiment, as shown in Figure 17. The values are shown in Table 10.

**Figure 17.** Performance measures plot for optimization using significant factors

**Figure 16.** Performance measures plot for optimization using all factors

**Table 8.** The best solution for optimization using all factors

**Table 9.** Significant factors for the first study object

The results found can be seen in Figure 17.

**Factor Variable Value** A Type 1 operators 2 B Type 2 operators 2 C Type 3 operators 1 D Type 1 inspection posts 2 E Type 2 inspection posts 1 F Type 3 inspection posts 1 G Type 4 inspection posts 1 H Type 5 inspection posts 2 J Type 6 inspection posts 1 K Type 7 inspection posts 1

**Factor Variables Value Range** B Type 2 operators 1 - 2 C Type 3 operators 1 - 2 D Type 1 inspection posts 1 - 2

Simrunner® converged after 8 experiments. The best result was 93, which was obtained in

the seventh experiment, as shown in Figure 17. The values are shown in Table 10.


Note: The parameters identified with (\*) were not used as input for optimization. They were kept at lower level (1).

**Table 10.** The best solution for optimization using significant factors

#### **9.2. Optimization of the second study object**

The optimization objective for the second simulation model was to find the best combination of model input variables in order to maximize the manufacturing cell's throughput. As presented in Table 3, the model possesses 12 input variables which are varied from the lower level (1) to the upper level (2).

In the first optimization phase, the 12 model input variables were selected and the optimization software was set up for experimentation. The results found can be seen in Figure 18.

The optimizer converged with 173 experiments. The best result found was 452,400, which was obtained in experiment 10 (Figure 18). The obtained values are shown in Table 11.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 95

The results are shown in Figure 19.

seen in Table 13.

**10. Results analysis** 

**10.1 First study object** 

**Figure 19.** Performance measures plot for optimization using significant factors

**Table 13.** Best solution for optimization using significant factors

Simrunner® converged after 31 experiments. The best value found was 449,800, which was attained in the eighth experiment carried out using the optimizer. The factors' values can be

**Factor Variable Value** A Type 1 operators 2 B Type 2 operators 2 C Type 3 operators 1\* D Type 1 machines 1\* E Type 2 machines 2 F Type 3 machines 1\* G Type 4 machines 2 H Type 5 machines 2 J Type 1 inspection posts 1\* K Type 2 inspection posts 1\* L Type 3 inspection posts 1\* M Type 4 inspection posts 1\* Note: The parameters identified with (\*) were not used as input for optimization. They were kept at lower level (1).

Table 14 presents a comparison of the results attained utilizing the three methods for the first study object. As far as the number of experiments executed, the advantage of using

**Figure 18.** Performance measures plot for optimization using all factors


**Table 11.** Best solution for optimization using all factors

The sensitivity analysis for the second study object identified five factors with significant effects which will be used for simulation optimization inputs (Table 12). The other model input variables were kept at their lower level (\*).


**Table 12.** Significant factors for the second study object

The results are shown in Figure 19.

94 Discrete Event Simulations – Development and Applications

**Figure 18.** Performance measures plot for optimization using all factors

**Table 11.** Best solution for optimization using all factors

input variables were kept at their lower level (\*).

**Table 12.** Significant factors for the second study object

**Factor Variable Value** A Type 1 operators 2 B Type 2 operators 2 C Type 3 operators 1 D Type 1 machines 2 E Type 2 machines 2 F Type 3 machines 2 G Type 4 machines 2 H Type 5 machines 1 J Type 1 inspection posts 2 K Type 2 inspection posts 2 L Type 3 inspection posts 2 M Type 4 inspection posts 2

The sensitivity analysis for the second study object identified five factors with significant effects which will be used for simulation optimization inputs (Table 12). The other model

> **Factor Variables Value range** A Number of type 1 operators 1 - 2 B Number of type 2 operators 1 - 2 E Number of type 2 machines 1 - 2 G Number of type 4 machines 1 - 2 H Number of type 5 machines 1 - 2

**Figure 19.** Performance measures plot for optimization using significant factors

Simrunner® converged after 31 experiments. The best value found was 449,800, which was attained in the eighth experiment carried out using the optimizer. The factors' values can be seen in Table 13.


Note: The parameters identified with (\*) were not used as input for optimization. They were kept at lower level (1).

**Table 13.** Best solution for optimization using significant factors

#### **10. Results analysis**

#### **10.1 First study object**

Table 14 presents a comparison of the results attained utilizing the three methods for the first study object. As far as the number of experiments executed, the advantage of using

sensitivity analysis to identify the significant factors becomes obvious. The commercial optimizer carried out 296 experiments when all of the input variables were chosen; when only the significant factors were utilized, merely 8 experiments were executed. Summing up the 64 experiments utilized with fractional factorial design, the result (72) is still four times smaller than the number of experiments executed with the optimizer when all 10 input variables were utilized.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 97

**Optimization using significative factors**

**Design of experiments** 

design with optimization. This possibility was not explored in detail here, as it was not this chapter's objective; however, many authors [1, 4, 18, 19] present optimization techniques

Table 15 shows a comparison of the results obtained using the three methods for the second

A 2 2 2 B 2 2 2 C 1 1\* 2 D 2 1\* 2 E 2 2 2 F 2 1\* 2 G 2 2 2 H 1 2 2 J 2 1\* 2 K 2 1\* 2 L 2 1\* 2 M 2 1\* 2 Result (WIP) 452400 449800 449800

**all factors**

Confidence Interval (95%) (445182 – 459617) (440960 – 458639) - Number of runs 173 31 64 Note: The parameters identified with (\*) were not used as input for optimization. They were kept at lower level (1).

In relation to the number of experiments executed, once again, the sensitivity analysis showed itself to be efficient. Along with the reduction in the number of factors, the number of experiments fell from 173 to 31. With the addition of 64 experiments using fractional factorial design, the method was efficient with 95 experiments, a little more than half of the

In relation to the optimization result, once again, it was necessary to perform a more detailed analysis of the responses. In spite of the average difference between the solutions presented with the optimizer (using all 12 and then only 5 significant factors) being around 2,600 pieces, it was verified that the solutions were within the same confidence interval. Thus the post-sensitivity analysis optimization solution's quality again showed itself to be

As was true with the last case, the results in Table 7 were found with only the use of fractional factorial design and were selected using the best results during the

**Table 15.** Optimization results for the three procedures analyzing the second study object

efficient when comparing the optimization of all of the input variables.

experiments when all factors were considered.

experimentation process.

**Parameter Optimization using**

using only experimental design.

**10.2. Second study object** 

study object.


Note: The parameters identified with (\*) were not used as input for optimization. They were kept at the lower level (1).

**Table 14.** Optimization results for the three procedures of the first study object

In respect to the responses found, it should be highlighted that, due to the simulation model's stochastic character, the response presented by the optimizer should be analyzed with care while considering both the average value and confidence interval of each result found.

Analyzing only the average optimization result found, it can be noted that the result found by the optimizer, when all 10 decision variables were manipulated was greater, reaching a WIP result of 92. However, when the response's confidence interval is analyzed, it can be said that the optimization's responses, when considering only the significant factors and their respective confidence intervals, are equal. The advantage of the response found using the sensitivity analysis is that only two factors (D and E) had to remain at the upper level, while the rest were kept at the lower level in order to minimize WIP.

The results in Table 5 were found with only the use of fractional factorial design and were selected using the best results during the experimentation process. This approach, however, does not take into consideration any simulation optimization approach and should be viewed with caution. The result using only DOE shows the possibility of using experimental design with optimization. This possibility was not explored in detail here, as it was not this chapter's objective; however, many authors [1, 4, 18, 19] present optimization techniques using only experimental design.

#### **10.2. Second study object**

96 Discrete Event Simulations – Development and Applications

variables were utilized.

**Parameter** 

Confidence Interval

found.

sensitivity analysis to identify the significant factors becomes obvious. The commercial optimizer carried out 296 experiments when all of the input variables were chosen; when only the significant factors were utilized, merely 8 experiments were executed. Summing up the 64 experiments utilized with fractional factorial design, the result (72) is still four times smaller than the number of experiments executed with the optimizer when all 10 input

> B 2 1 1 C 1 1 1 D 2 2 2 E 1 2 2

Result (WIP) 92 93 93

(95%) (83 – 100) (86 – 99) - Number of runs 296 8 64 Note: The parameters identified with (\*) were not used as input for optimization. They were kept at the lower level (1).

In respect to the responses found, it should be highlighted that, due to the simulation model's stochastic character, the response presented by the optimizer should be analyzed with care while considering both the average value and confidence interval of each result

Analyzing only the average optimization result found, it can be noted that the result found by the optimizer, when all 10 decision variables were manipulated was greater, reaching a WIP result of 92. However, when the response's confidence interval is analyzed, it can be said that the optimization's responses, when considering only the significant factors and their respective confidence intervals, are equal. The advantage of the response found using the sensitivity analysis is that only two factors (D and E) had to remain at the upper level,

The results in Table 5 were found with only the use of fractional factorial design and were selected using the best results during the experimentation process. This approach, however, does not take into consideration any simulation optimization approach and should be viewed with caution. The result using only DOE shows the possibility of using experimental

**Optimization using significant factors** 

**Design of experiments** 

1

1

1

1

2

1

**Optimization using all factors** 

A 2 1\*

F 1 1\*

G 1 1\*

H 2 1\*

J 1 1\*

K 1 1\*

**Table 14.** Optimization results for the three procedures of the first study object

while the rest were kept at the lower level in order to minimize WIP.

Table 15 shows a comparison of the results obtained using the three methods for the second study object.


Note: The parameters identified with (\*) were not used as input for optimization. They were kept at lower level (1).

**Table 15.** Optimization results for the three procedures analyzing the second study object

In relation to the number of experiments executed, once again, the sensitivity analysis showed itself to be efficient. Along with the reduction in the number of factors, the number of experiments fell from 173 to 31. With the addition of 64 experiments using fractional factorial design, the method was efficient with 95 experiments, a little more than half of the experiments when all factors were considered.

In relation to the optimization result, once again, it was necessary to perform a more detailed analysis of the responses. In spite of the average difference between the solutions presented with the optimizer (using all 12 and then only 5 significant factors) being around 2,600 pieces, it was verified that the solutions were within the same confidence interval. Thus the post-sensitivity analysis optimization solution's quality again showed itself to be efficient when comparing the optimization of all of the input variables.

As was true with the last case, the results in Table 7 were found with only the use of fractional factorial design and were selected using the best results during the experimentation process.

## **11. Conclusion**

The objective of this chapter was to present how experimental design and analysis techniques can be employed in order to identify significant variables in discrete-event simulation models, thus aiding simulation optimization searches for optimal solutions.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 99

use of Kriging metamodeling for simulation has established itself in the scientific simulation community, which can be seen in [46, 47, 48, 49], thus demonstrating its promising research

The use of discrete-event simulation along with optimization is still scarce; nonetheless, in the last decade, important studies about this area of operational research have started to be realized, supporting the wider acceptance of this approach while also investigating the barriers to its continuous improvement. The use of sensitivity analysis in DOE enables a reduction in search space while increasing the optimization process's efficiency and speed.

Simulation optimization helps take simulation from being merely a means of scenario evaluation to a much greater solution generator. In doing so, sensitivity analysis plays a crucial role in this process, as it helps overcome the time and computational potential barriers presented by simulation models with large numbers of variables, thus making

In respect to future research possibilities, a potentially rich area in terms of investigation would be the examination of experimental designs which would reduce the number of experiments needed in order to identify the significant factors in pre-optimization phases. Another point which could be investigated more profoundly is the inclusion of qualitative techniques, such as brainstorming, cause and effect diagrams and Soft Systems Methodology, in order to select the factors to be utilized in experimentation. A field which is

little-explored is sensitivity analysis in the optimization of multiple-objective models.

José Arnaldo Barra Montevechi, Rafael de Carvalho Miranda and Jonathan Daniel Friend *Universidade Federal de Itajubá (UNIFEI), Instituto de Engenharia de Produção e Gestão (IEPG),* 

The authors extend their sincere gratitude to the Brazilian educational funding agencies of FAPEMIG, CNPq, the engineering support program CAPES, and the company PADTEC for

[1] Banks J, Carson II JS, Nelson BL, Nicol DM (2005) Discrete-event Simulation. New

[2] Bruzzone AG, Bocca E, Longo F, Massei M (2007) Training and recruitment in logistics node design by using web-based simulation. Int. J. Internet Manuf. Serv. 1(1): 32-50. [3] Sargent RG (2009) Verification and validation of simulation models. In: Winter

optimization an even greater tool for aiding decision-making.

their continued support during the development of this project.

Simulation Conference, Proceedings... Austin, TX, USA.

[4] Law AM (2007) Simulation modeling and analysis. New York: McGraw-Hill.

field.

**Author details** 

*Itajubá, MG, Brazil* 

**12. References** 

Jersey: Prentice-Hall.

**Acknowledgement** 

To develop this application, the main concepts of simulation optimization were presented during this chapter. In following, two applications were developed to verify how fractional factorial design can be used in sensitivity analysis for simulation models, identifying its advantages, disadvantages and effects on model optimization.

For optimization, the identification of the significant variables was extremely important, as it enabled the reduction of the search space and computational potential necessary to perform the search for an optimal solution. Optimization was performed for the two applications in two distinct forms, utilizing simulation optimization.

The first simulation optimization approach was to use all of the models' input variables. No previous studies were performed to determine if the models' variables exercised significant effects on overall system performance. The second approach relied on sensitivity analysis in order to identify the variables which influenced system performance. After identifying the significant variables, model optimization was utilized while using the reduced search space. The third approach involved using only the experimentation's results without using any simulation optimization procedure.

By analyzing the results, the advantages of using sensitivity analysis become evident, not only due to the reduction in necessary computational potential for the optimization process, but also for the greater level of detail and knowledge acquired about the process under study. By using this experimentation, it is possible to verify the process's variables which exercise the greatest effects on overall system performance, thus being able to determine the effect each variable has on the process as well as their interactions. These interactions would be very difficult to define and easy to disregard in simulation projects without the use of DOE.

It should not go without noting, however, that one must take extreme caution during execution of these experiments, as just one experiment realized under incorrect conditions or out of the correct matrix order can lead to erroneous results. In order to analyze the results obtained during experimentation, the user should have a solid understanding of DOE and statistics. Those erroneous conclusions in simulation could possibly lead to incorrect implementation, which could generate very tangible costs in the real world. Thus, it is recommended that one research further the concepts shown in this chapter.

One approach that is commonly used for simulation model optimization which was not explored in great detail in this chapter is the development of a mathematical metamodel, which represents a determined model output, according to optimization. This approach has a vast field of application and could have been applied in the problems presented. It is also recommended that the reader study this topic further as well [18, 43, 46]. In this sense, the use of Kriging metamodeling for simulation has established itself in the scientific simulation community, which can be seen in [46, 47, 48, 49], thus demonstrating its promising research field.

The use of discrete-event simulation along with optimization is still scarce; nonetheless, in the last decade, important studies about this area of operational research have started to be realized, supporting the wider acceptance of this approach while also investigating the barriers to its continuous improvement. The use of sensitivity analysis in DOE enables a reduction in search space while increasing the optimization process's efficiency and speed.

Simulation optimization helps take simulation from being merely a means of scenario evaluation to a much greater solution generator. In doing so, sensitivity analysis plays a crucial role in this process, as it helps overcome the time and computational potential barriers presented by simulation models with large numbers of variables, thus making optimization an even greater tool for aiding decision-making.

In respect to future research possibilities, a potentially rich area in terms of investigation would be the examination of experimental designs which would reduce the number of experiments needed in order to identify the significant factors in pre-optimization phases. Another point which could be investigated more profoundly is the inclusion of qualitative techniques, such as brainstorming, cause and effect diagrams and Soft Systems Methodology, in order to select the factors to be utilized in experimentation. A field which is little-explored is sensitivity analysis in the optimization of multiple-objective models.

## **Author details**

98 Discrete Event Simulations – Development and Applications

simulation optimization procedure.

DOE.

advantages, disadvantages and effects on model optimization.

applications in two distinct forms, utilizing simulation optimization.

The objective of this chapter was to present how experimental design and analysis techniques can be employed in order to identify significant variables in discrete-event simulation models, thus aiding simulation optimization searches for optimal solutions.

To develop this application, the main concepts of simulation optimization were presented during this chapter. In following, two applications were developed to verify how fractional factorial design can be used in sensitivity analysis for simulation models, identifying its

For optimization, the identification of the significant variables was extremely important, as it enabled the reduction of the search space and computational potential necessary to perform the search for an optimal solution. Optimization was performed for the two

The first simulation optimization approach was to use all of the models' input variables. No previous studies were performed to determine if the models' variables exercised significant effects on overall system performance. The second approach relied on sensitivity analysis in order to identify the variables which influenced system performance. After identifying the significant variables, model optimization was utilized while using the reduced search space. The third approach involved using only the experimentation's results without using any

By analyzing the results, the advantages of using sensitivity analysis become evident, not only due to the reduction in necessary computational potential for the optimization process, but also for the greater level of detail and knowledge acquired about the process under study. By using this experimentation, it is possible to verify the process's variables which exercise the greatest effects on overall system performance, thus being able to determine the effect each variable has on the process as well as their interactions. These interactions would be very difficult to define and easy to disregard in simulation projects without the use of

It should not go without noting, however, that one must take extreme caution during execution of these experiments, as just one experiment realized under incorrect conditions or out of the correct matrix order can lead to erroneous results. In order to analyze the results obtained during experimentation, the user should have a solid understanding of DOE and statistics. Those erroneous conclusions in simulation could possibly lead to incorrect implementation, which could generate very tangible costs in the real world. Thus,

One approach that is commonly used for simulation model optimization which was not explored in great detail in this chapter is the development of a mathematical metamodel, which represents a determined model output, according to optimization. This approach has a vast field of application and could have been applied in the problems presented. It is also recommended that the reader study this topic further as well [18, 43, 46]. In this sense, the

it is recommended that one research further the concepts shown in this chapter.

**11. Conclusion** 

José Arnaldo Barra Montevechi, Rafael de Carvalho Miranda and Jonathan Daniel Friend *Universidade Federal de Itajubá (UNIFEI), Instituto de Engenharia de Produção e Gestão (IEPG), Itajubá, MG, Brazil* 

## **Acknowledgement**

The authors extend their sincere gratitude to the Brazilian educational funding agencies of FAPEMIG, CNPq, the engineering support program CAPES, and the company PADTEC for their continued support during the development of this project.

## **12. References**

	- [5] Jahangirian M, Eldabi T, Naseer A, Stergioulas LK, Young T (2010) Simulation in manufacturing and business: A review. Eur J Oper Res. 203(1):1-13.

Sensitivity Analysis in Discrete-Event Simulation Using Design of Experiments 101

[25] Coello CAC, Lamont GB, Van Veldhuizen DA (2007) Evolutionary Algorithms for Solving Multi-Objective Problems (Genetic and Evolutionary Computation). New York:

[26] Holland JH (1992) Adaptation in Natural and Articial Systems. Cambridge: MIT Press. [27] Martí R, Laguna M, Glover F (2006) Principles of Scatter Search. Eur J Oper Res.169(2):

[28] Glover F, Laguna M, Martí R (2005) Principles of Tabu Search, In: Gonzalez T, editor, Approximation Algorithms and Metaheuristics. London: Chapman & Hall/CRC. [29] Ripley B (1996) Pattern Recognition and Neural Networks. Cambridge: University

[30] Aarts EHL, Korst J, Michiels W (2005) Simulated Annealing. In: Burke EK, Kendall G, editors. Introductory tutorials in optimisation, decision support and search

[31] April J, Glover F, Kelly JP, Laguna M (2003) Practical introduction to simulation optimization. In: Winter Simulation Conference, Proceedings... New Orleans, LA, USA. [32] Banks, J. Panel Session: The Future of Simulation. In: Winter Simulation Conference,

[33] Tyni T, Ylinen J (2006) Evolutionary bi-objective optimization in the elevator car routing

[34] Montgomery DC (2009) Introduction to Statistical Quality Control. New York: John

[35] Montgomery DC, Runger GC (2003) Applied Statistics and Probability for Engineers.

[36] Kelton WD (2003) Designing simulation experiments. In: Winter Simulation

[37] Kleijnen JPC, Sanchez SM, Lucas TW, Cioppa TM (2005) State-of-the-Art Review: A User's Guide to the Brave New World of Designing Simulation Experiments. J. Comput.

[38] Montevechi JAB, Pinho AF, Leal F, Marins FAZ (2007) Application of design of experiments on the simulation of a process in an automotive industry. In: Winter

[39] Sanchez SM, Moeeni F, Sanchez PJ (2006) So many factors, so little time…. Simulation

[40] Kleijnen JPC (2001) Experimental designs for sensitivity analysis of simulation models.

[41] Chung CA (2004) Simulation Modeling Handbook: a practical approach. Washington,

[42] Landsheer JA, Wittenboer GVD, Maassen GH (2006) Additive and multiplicative effects in a fixed 2 x 2 design using ANOVA can be difficult to differentiate: demonstration

[43] Montevechi JAB, Almeida Filho RG, Paiva AP, Costa RFS, Medeiros AL (2010) Sensitivity analysis in discrete-event simulation using fractional factorial designs. J.

methodologies. New York: Springer. pp. 187-211.

Conference, Proceedings... New Orleans, LA, USA.

In: Eurosim, Proceedings… Delft, Netherlands.

and mathematical reason. Soc Sci Res. 35: 279–294.

Simulation Conference, Proceedings... Washington, DC, USA.

experiments in the frequency domain. Int J Prod Econ. 103: 149–165.

Proceedings... Arlington, VA, USA, 2001.

problem. Eur J Oper Res. 169(3):960-977.

New York: John Wiley & Sons, Inc.

Wiley & Sons, Inc.

17(3): 263–289.

DC: CRC Press.

Simulat. 4: 128-142.

Springer.

359-372.

Press.


[25] Coello CAC, Lamont GB, Van Veldhuizen DA (2007) Evolutionary Algorithms for Solving Multi-Objective Problems (Genetic and Evolutionary Computation). New York: Springer.

100 Discrete Event Simulations – Development and Applications

215.

Hill.

McGraw-Hill.

New Jersey: Prentice-Hall.

John Wiley & Sons, Inc.

Eng. 54: 327‐339.

458.

Proceedings... Monterey, CA, USA.

Conference, Proceedings... San Diego, CA, USA.

using simulation. Utha: Promodel Corporation.

Simulation Conference, Proceedings... Orlando, FL, USA.

Simulation. New York: John Wiley & Sons. pp. 307-333.

Conference, Proceedings... San Diego, CA, USA.

Conference, Proceedings... Dallas, TX, USA.

York: John Wiley & Sons. p. 173-223.

[5] Jahangirian M, Eldabi T, Naseer A, Stergioulas LK, Young T (2010) Simulation in

[8] Harrel CR, Mott JRA, Bateman RE, Bowden RG, Gogg TJ (1996) System improvement

[9] Fu MC (2002) Optimization for Simulation: Theory vs. Practice. J. Comput. 14(3): 192-

[10] Fu MC, Andradóttir S, Carson JS, Glover F, Harrell CR, Ho, YC, Kelly JP, Robinson SM (2000) Integrating optimization and simulation: research and practice. In: Winter

[11] Harrel CR, Ghosh BK, Bowden R (2004) Simulation Using Promodel. New York:

[12] Law AM, Kelton WD (2000) Simulation modeling and analysis. New York: McGraw-

[13] Banks J, Carson II JS, Nelson BL, Nicol DM (2000) Discrete event system simulation.

[14] April J, Better M, Glover F, Kelly JP, Laguna M (2005) Enhancing business process management with simulation optimization. In: Winter Simulation Conference,

[15] Andradóttir S (1998) Simulation optimization. In: Banks J, editor, Handbook of

[16] Biles WE (1979) Experimental design in computer simulation. In: Winter Simulation

[17] Biles WE (1984). Experimental design in computer simulation. In: Winter Simulation

[18] Montgomery DC (2005) Design and Analysis of Experiments. New York: New York:

[19] Kleijnen JPC (1998). Experimental design for sensitivity analysis, optimization, and validation of simulation models. In: Banks J, editor. Handbook of Simulation. New

[20] Carson Y, Maria A (1997) Simulation optimization: methods and applications. In:

[21] Azadeh A, Tabatabaee M, Maghsoudi A (2009) Design of Intelligent Simulation Software with Capability of Optimization. Aust. J. Basic Appl. Sci. 3(4): 4478-4483. [22] Fu MC (1994) Optimization via simulation: A review. Ann Oper Res. 53:199-247.

[23] Rosen SL, Harmonosky CH, Traband MT (2007) Optimization of Systems with Multiple Performance Measures via Simulation: Survey and Recommendations. Comput. Ind.

[24] Bettonvil BWM, Castillo E, Kleijnen JPC (2009) Statistical testing of optimality conditions in multiresponse simulation-based optimization. Eur J Oper Res. 199(2): 448-

Winter Simulation Conference, Proceedings... Atlanta, GA, USA.

[6] Ryan J, Heavey C (2006) Process modeling for simulation. Comput Ind. 57(5): 437-450. [7] Law AM, Mccomas, MG (2002) Simulation-Based Optimization, In: Winter Simulation

manufacturing and business: A review. Eur J Oper Res. 203(1):1-13.

	- [44] Ye KQ, Hamada M (2001) A step-down Lenth method for analyzing unreplicated factorial designs. J Qual Technol. 33(2):140–153.

**Section 2** 

**Novel Integration of Discrete Event** 

**Simulation with Other Modeling Techniques** 

