# A New Approach to Reducing Search Space and Increasing Efficiency in Simulation Optimization Problems via the Fuzzy-DEA-BCC

Authors: Rafael de Carvalho Miranda; José Arnaldo Barra Montevechi; Aneirson Francisco da Silva; Fernando Augusto Silva Marins
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The development of discrete-event simulation software was one of the most successful interfaces in operational research with computation. As a result, research has been focused on the development of new methods and algorithms with the purpose of increasing simulation optimization efficiency and reliability. This study aims to define optimum variation intervals for each decision variable through a proposed approach which combines the data envelopment analysis with the Fuzzy logic (Fuzzy-DEA-BCC), seeking to improve the decision-making units’ distinction in the face of uncertainty. In this study, Taguchi’s orthogonal arrays were used to generate the necessary quantity of DMUs, and the output variables were generated by the simulation. Two study objects were utilized as examples of mono- and multiobjective problems. Results confirmed the reliability and applicability of the proposed method, as it enabled a significant reduction in search space and computational demand when compared to conventional simulation optimization techniques.


---

## Body

## 1. Introduction

The development of discrete-event simulation (DES) software was one of the greatest successes in bringing the realms of operational research (OR) and computation together, according to Fu [1]. Hillier and Lieberman [2] highlight that simulation is an extremely versatile technique which enables the investigation of practically any type of stochastic system. This versatility has turned simulation into the most commonly used OR technique for stochastic systems.Nevertheless, the simulation optimization integration grew stronger from the 90s due to the development of commercial software packages which marketed integrated optimization routines, thus making it considerably easier to carry out a decision-making analysis [2–4].Thus, Azadeh et al. [5] affirm that simulation optimization is one of the most important OR tools that has come about in recent years. Previous methodologies demanded complex alterations which were frequently economically and temporally unviable, especially for problems with a large number of decision variables.According to Medaglia et al. [6], simulation optimization aims to find the best values for simulation model input parameters in the search of one or more desired outputs. This is generally a slow process that takes a large amount of time and the accomplishment of innumerous experiments.In spite of the advances in simulation model optimization software, a common criticism brought up is that, when dealing with more than one input variable, the optimization process becomes very slow [7–9]. Furthermore, according to Hillier and Lieberman [2], computational simulation packages may be considered as relatively slow and costly when applied in studies of stochastic and dynamic systems. In such systems where a random behavior is prevalent, there is a tendency towards elevated expenses and allocation of qualified labor and expertise, along with extra time for analysis and programming. All these factors lead to a considerably higher computational demand.Following this thought, Kleijnen et al. [10] recognize that simulation optimization problems are generally difficult to be resolved and present disadvantages such as the fact that the model outputs are implicit functions and exposed to noise. These authors highlight the difficulties involved in analyzing outputs for stochastic simulation models due to the existing variation in each replication.Taking all of this into account, this study aims to propose a procedure for the simulation model optimization which reduces search space and enables an increase in efficiency, measured through reduced computational demand by means of efficiency evaluation in the face of uncertainty.In order to do meet this paper’s premises, the procedure herein proposed makes use of the Taguchi orthogonal arrays [11] to represent the experimental region and to test each scenario and its replications using DES and Fuzzy logic combined with DEA. The method is called Fuzzy-DEA-BCC [12, 13] and incorporated the concept of superefficiency. A further description is offered later on in the paper.The use of the Taguchi orthogonal arrays [11] as means of representing the experimental region of a simulation optimization program is justified by the reduced number of experiments for such experimental designs, given that they are saturated fractional factorial arrays, thus allowing the analysis of k factors with n levels, testing all levels for each factor in a balanced manner [14].The arrays for this method were chosen in function of the number of inputs and outputs to follow the rule which determines that, in order to utilize classic DEA (CCR and BCC) models, the minimum number of decision-making units (DMUs) must be equal to or greater than three times the sum of the total number of input and output variables [12, 15, 16].According to Cook and Seiford [17], the DEA provides a method which allows an identification of the DMUs which serve as benchmark for the others under analysis, forming an efficiency boundary. In the specific case of the DES, Weng et al. [18] assert that the DEA allows the evaluation of relative efficiency of a group of entities using multiple inputs and outputs without being aware of their relationship.For this study, Fuzzy-DEA-BCC models were chosen due to the stochastic and nonlinear nature of the DES and the fact that the set of DMUs may generally present different characteristics (number of personnel, machines, throughput, profit, work in progress, leading time, etc.) and tend to have different yields or outputs on different scales. The latter justifies the use of the DEA-BCC model [18].Finally, in order to rank the scenarios generated by the experimental matrix and reduce search space, the concept of superefficiency was used, as proposed by Andersen and Petersen [19]. Xue and Harker [20] point out that superefficiency is capable of differentiating DMUs, thus permitting a classification in terms of efficiency.In order to reach the proposed objectives, this paper is divided into 5 sections. Section2 presents the theoretical base, focusing on the DEA, and Section 3 describes the optimization procedure, integrating the previously presented tools. Section 4 details the application of the new optimization procedure and discusses the results. Finally, Section 5 presents the study conclusions.


## 2. Data Envelopment Analysis

### 2.1. Deterministic Data Envelopment Analysis

Classic DEA models were introduced by Charnes et al. [15], which use constant returns of scale and were denominated as CCR models in homage to their creators [12]. These models were then extended by Banker et al. [12] to variable returns of scale, this time dubbed as BCC, again in reference to their originators.According to Cook and Seiford [17], DEA is a nonparametric methodology which comparatively measures the efficiency of each decision-making unit (DMU). Another fact that deserves highlighting is that, through its use, it is possible to avoid the problems created by incommensurability (different units of measurement) among the elements of the input and output matrices.In the original model from Charnes et al. [15], the input and output variable weights may be obtained from the Fractional Programming model solution, given by
(1)
E
j
=
∑
r
=
1
s
u
r
·
y
r
0
∑
i
=
1
m
v
i
·
x
i
0
subject to:
(2)
∑
r
=
1
s
u
r
y
r
j
∑
i
=
1
m
v
i
x
i
j
≤
1
,
j
=
1,2
,
…
,
n
,
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
v
i
≥
0
,
i
=
1,2
,
…
,
m
.With DMU0 to DMU being under evaluation, E
j is the relative efficiency of DMU0; y
r
0 and x
i
o are the input and output data for DMU0; j is the index for the DMU, j
=
1,2
,
…
n; r is the output index, with r
=
1,2
,
…
,
s; i is the input index, i
=
1,2
,
…
,
m; y
r
j is the value of rth output for jth DMU; x
i
j is the value of ith input for jth DMU; u
r is the weight associated with rth output; v
i is the weight associated with ith input.It is observed that ifE
j
=
1, DMU0 is going to be efficient when compared to the other units considered in the model, and, if 0
<
E
j
<
1, this DMU is deemed inefficient.This model (1)-(2) is not linear, thus having multiple solutions; nonetheless, it can be linearized, generating a DEA-CCD model, given by the following equation:
(3)
max
⁡
E
j
=
∑
r
=
1
s
u
r
·
y
r
0
subject
to:
∑
i
=
1
m
v
i
x
i
o
=
1
h
h
h
h
h
∑
r
=
1
s
u
r
y
r
j
-
∑
i
=
1
m
v
i
x
i
j
≤
0
,
j
=
1,2
,
…
,
n
,
h
h
h
h
h
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
h
h
h
h
h
v
i
≥
0
,
i
=
1,2
,
…
,
m
.DEA-BCC models, with variable returns of scale, can be expressed by the following equations:(4)
max
⁡
E
j
=
∑
r
=
1
s
u
r
·
y
r
0
+
c
0
(5)
subject
to:
∑
i
=
1
m
v
i
x
i
o
=
1
h
h
h
h
h
h
(6)
h
h
h
h
h
h
h
h
h
h
∑
r
=
1
s
u
r
y
r
j
-
∑
i
=
1
m
v
i
x
i
j
+
c
0
≤
0
,
h
h
h
h
h
h
h
h
h
h
h
j
=
1,2
,
…
,
n
.
(7)
h
h
h
h
h
h
h
h
h
u
r
≥
0
,
r
=
1,2
,
…
,
s
.
(8)
h
h
h
h
h
h
h
h
h
v
i
≥
0
,
i
=
1,2
,
…
,
m
.
(9)
c
0
free
.
hhhhA main difference between DEA-BCC and DEA-CCR models is the addition ofc
0 (free) variable which indicates variable returns of scale. Banker et al. [12] affirm that a DMU considered efficient in a BCC model will also be considered efficient in the CCR model; the inverse, however, is not necessarily true.According to Cooper et al. [16], in order to provide a suitable discrimination of the DMUs in traditional DEA models, the following must be verified: number of DMUs ≥ maximum {(product of the number of inputs and outputs), 3·(number of inputs + number of outputs)}.Classic DEA models consider DMUs withE
j
=
1 as efficient and DMUs with 0
≤
E
j
<
1 as inefficient. It is within the realm of possibility that multiple DMUs could be efficient; that is, this would not be a good discrimination of DMUs. In order to deal with this limitation, Andersen and Petersen [19] proposed the concept of superefficiency in order to help differentiate DMUs which present E
j
=
1.For the superefficiency evaluation to be employed in DEA-BCC models, constraint (6) must be removed in order for the DMU under analysis to attain greater scores than 1, being that the DMUs considered inefficient in traditional evaluation models will continue to be inefficient, but those which posed scores equal to one will be able to demonstrate scores above 1, thus enabling the elaboration of a ranking.


### 2.2. Fuzzy-DEA

Hatami-Marbini et al. [21] undertook an expansive literature review about the Fuzzy Theory combined with DEA models. The authors’ motivations were based on the fact that, in general, the estimation of input and output DMU values in real problems is difficult, which could generate efficiency values with a low level of reliability; one approach to deal with the aspects of uncertainty involves the adoption of Fuzzy Theory concepts.Kao and Liu [13] assert that measuring a DMU’s efficiency is a difficult task which involves complex economic variables such as interest rates and tax rates and employment levels and demand. According to these authors, efficiency measurement becomes even more difficult when analyzing multiple inputs and outputs.In this context, Wen et al. [22] comment that DMUs can be classified in two categories: efficient and inefficient. However, the incorporation of uncertainty as an error in measurement of inputs and outputs can make the calculation of efficiency more reliable and robust.Fuzzy Theory has come to be used with the objective of modeling uncertainty in DEA models [23]. Fuzzy-DEA models are based on Fuzzy Linear Programming, with DEA-BCC models with Fuzzy coefficients having special emphasis in this paper (Fuzzy-DEA-BCC), proposed by Kao and Liu [13]. Consider the following:
(10)
E
~
j
=
Max
⁡
∑
r
=
1
s
u
r
y
~
r
0
+
c
0
subject to:
(11)
∑
i
=
1
m
x
~
i
0
·
v
i
=
1
(12)
∑
r
=
1
s
u
r
y
~
r
j
-
∑
i
=
1
m
v
i
x
~
i
j
+
c
0
≤
0
,
j
=
1,2
,
…
,
n
,
(13)
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
(14)
v
i
≥
0
,
i
=
1,2
,
…
,
m
.The value of the objective function (10) may have a greater value than 1 due to constraints (11)-(12) which involve Fuzzy parameters and are solved using probability [23]. With the incorporation of Fuzzy coefficients, DEA-BCC models cannot be resolved using traditional linear programming (LP) techniques. Hatami-Marbini et al. [21] list and describe the following main approaches which deal with Fuzzy-DEA:(i)
theα-level based approach;
(ii)
the tolerance approach;
(iii)
the Fuzzy ranking approach;
(iv)
the possibility approach.In this study, the approach based on theα-level was adopted and is described afterwards. This α-level application is the most common for Fuzzy-DEA models according to Hatami-Marbini et al. [21]. To apply this method, the idea is to convert the Fuzzy-DEA model into a pair of Parametric Programming Problems to find upper and lower boundaries for the DMU efficiency score membership functions [23].For the purposes of this research, triangular membership functions were utilized. According to Liang and Wang [24], such functions well represent human expertise in adequately judging the behavior of common variables in a range of practical situations. Along these same lines of thought, Aouni et al. [25] show multiple applications for Fuzzy triangular numbers which validate and justify the adoption of such a method in conjunction with goal programming (GP) models. Another justification arises from the fact that it is a linear function, thus easing the optimization process through traditional LP means [21].To illustrate the process of Fuzzy-DEA modeling, the example in Kao and Liu’s [13] paper was used, in which they dealt with a model with a variable return of scale. Figure 1 shows the DEA-CCR and DEA-BCC models with four DMUs, denominated as A, B, C, and D, with only a single input and a single output for each DMU and with the output of DMU B being Fuzzy.Figure 1
Production frontier when DMU B output is Fuzzy. Source: adapted from Kao and Liu [13].The input values for DMUs A, B, C, and D are 10, 20, 30, and 50, respectively. With these inputs, outputs 5, 9, and 15 from DMUs A, C, and D are also produced, respectively. These inputs and outputs are associated with points A = (10; 5), C = (30; 9), and D = (50; 15) in Figure1. For DMU B, the output associated with input value 20 is represented by the trapezoidal Fuzzy number (5; 6; 8; 9), which is also illustrated in Figure 1 by points (20; 5), (20; 6), (20; 7.5), (20; 8), and (20; 9).Based on Figure1 and according to the DEA-BCC model [13], when the output of DMU B is less than or equal to 7.5, the production boundary is defined by the line segments connecting point (10; 0) to A = (10; 5) and from A to D = (50; 15). The efficiency scores of DMUs A, C, and D are 1, 0.9, and 1, respectively; the efficiency of DMU B is between 5/7.5 = 0.67 and 7.5/7.5 = 1, depending on its output.If an increase in the output associated with DMU B went from 7.5 to 9, the production boundary would be represented by the line segments linking points (10; 0) to A from A to B = (20; 9) and from B to D = (50; 15). In this case, DMUs A, B, and D have an efficient output bound, while the efficiency of C will be, according to output orientation, between values of 9/10= 0.9 and 9/11 = 0.82.In this example, regardless of the output value of DMU B, DMUs A and D are always efficient. In other words, with the combination (generation of scenarios) of each Fuzzy output for DMU B, the effects of uncertainty on DMUs A and D are always the same and present 100% efficiency.Conducting an analysis under the assumption of constant return of scale (DEA-CCR model), as seen in Figure1, the production boundary will be the solid line which links the origin point (0; 0) with point A = (10; 6). In this case, only DMU A will always be efficient (with an efficiency score of 100%), regardless of the output value of DMU B. However, the efficiency of DMU B will vary from 1/2 = 0.5 to 9/10 = 0.9, and DMUs C and D will have an efficiency of 3/5 = 0.6; that is, there are not effects of the Fuzzy output of DMU B over the efficiencies of DMUs C and D.As a further example, the membership function linked to the output of DMU B may be illustrated, given that it is a trapezoidal function, shown by the following equation:(15)
μ
Y
~
B
(
y
)
=
{
y
-
5
,
5
≤
y
≤
6
,
1
,
6
≤
y
≤
7
,
9
-
y
,
8
≤
y
≤
9
.Theα-level base is defined as the interval [
5
+
α
,
9
-
α
], being (
E
B
)
α
L
e
(
E
B
)
α
U, respectively. The upper and lower bounds for the variation (uncertainty) of α-level, which are associated with a scenario generation, vary between pessimistic (α
=
0) and optimistic (α
=
1) for an efficiency analysis.As previously stated, theα-level based approach was adopted which, according to Hatami-Marbini et al. [21], is the most popular Fuzzy-DEA model with many referenced publications [13]. The value of α
∈
[
0
,
1
] allows the generation of scenarios, that is, different efficiency values respecting the variation range determined by the membership function. In such models, X
~
i
j and Y
~
r
j are, respectively, the Fuzzy parameters for ith input and rth output of the jth DMU. These values are approximately known and can be represented by Fuzzy sets through means of membership functions μ
X
~
i
j and μ
Y
~
r
j. Thus, Fuzzy-DEA models can be formed using S
(
X
~
i
j
) and S
(
Y
~
r
j
) which are the values of these parameters for a given S scenario. With the α-level which generates a set of scenarios for Y
~
r
j and X
~
i
j, the formulas, as defined by Kao and Liu [13], are as follows:
(16)
(
X
i
j
)
α
=
{
x
i
j
∈
S
(
X
~
i
j
)
∣
μ
X
~
i
j
(
x
i
j
)
≥
α
}
,
h
h
h
h
h
h
h
h
i
i
i
=
1,2
,
…
m
,
j
=
1,2
,
…
,
n
,
(
Y
r
j
)
α
=
{
y
r
j
∈
S
(
Y
~
r
j
)
∣
μ
Y
~
r
j
(
y
r
j
)
≥
α
}
,
h
h
h
h
h
h
h
h
i
r
=
1,2
,
…
,
s
,
j
=
1,2
,
…
,
n
.Fuzzy-DEA can be transformed into a family of DEA models with different levels of uncertainty:{
(
X
i
j
)
α
∣
0
≤
α
≤
1
} and {
(
Y
r
j
)
α
∣
0
≤
α
≤
1
}. The results for each scenario identify the uncertainty variation range in the model’s input and output data [13, 26–28], with the α-level being defined by the following equation, according to Kao and Liu [13]:
(17)
(
X
i
j
)
α
=
[
{
min
⁡
x
i
j
x
i
j
∈
S
(
X
~
i
j
)
∣
μ
X
~
i
j
(
x
i
j
)
≥
α
}
,
h
h
h
h
max
⁡
x
i
j
{
x
i
j
∈
S
(
X
~
i
j
)
∣
μ
X
~
i
j
(
x
i
j
)
≥
α
}
]
,
(
Y
r
j
)
α
=
[
{
min
⁡
y
r
j
y
r
j
∈
S
(
Y
~
r
j
)
∣
μ
Y
~
r
j
(
y
r
j
)
≥
α
}
,
h
h
h
h
max
⁡
y
r
j
{
y
r
j
∈
S
(
Y
~
r
j
)
∣
μ
Y
~
r
j
(
y
r
j
)
≥
α
}
]
.Kao and Liu [13], based on Yager [29], Zadeh [28], and Zimmermann [30], established that a membership function which defines the efficiency of DMU j can be expressed by
(18)
μ
E
~
j
(
z
)
=
sup
⁡
x
,
y
min
⁡
=
{
μ
X
~
i
j
(
x
i
j
)
,
μ
Y
~
r
j
(
y
r
j
)
,
∀
i
,
r
,
j
∣
z
=
E
j
(
x
,
y
)
}
,
with E
j
(
x
,
y
) obtained by (4)–(9). The approach for constructing the membership function μ
E
~
j proposed in this study adopted an α-level μ
E
~
j as being the value of efficiency in each scenario E
~
j obtained by (10)–(14). For further details, Kao and Liu [13], Hatami-Marbini et al. [21], and Kao and Lin [31] are recommended reading materials.Based on the models developed by Banker et al. [12], in accordance with (4)–(9), and Kao and Liu [13], in accordance with (10)–(14), the Fuzzy-DEA-BCC model was developed. Afterwards, the indices, parameters, auxiliary variables and decision variables, objective functions, and model constraints are proposed, considering DMU0 as the one under analysis.(i)
Indices:
(a)
j is the DMU index.
(b)
r is the output index.
(c)
i is the input index.
(ii)
Parameters:
(a)
y
~
r
0 and x
~
i
0 are, respectively, the upper and lower bounds in the definition intervals of the triangular membership function for the rth Fuzzy output and the ith Fuzzy input for DMU0, considering the average as the most probable value without uncertainty.
(b)
Y
~
r
0 and X
~
i
0 are, respectively, the upper and lower bounds in the definition intervals of the triangular membership function for the rth Fuzzy output and the ith Fuzzy input for DMU0, considering the average as the most probable value without uncertainty.
(c)
y
~
r
j is the lower bound for the definition interval of the triangular membership function of the rth Fuzzy output for the jth DMU, considering the average as the most probable value without uncertainty.
(d)
Y
~
r
j is the upper bound of the triangular membership function definition interval of the rth Fuzzy output for the jth DMU, considering the average as the most probable value without uncertainty.
(e)
x
~
i
j is the lower bound of the triangular membership function definition interval of the rth Fuzzy output for the jth DMU, considering the average as the most probable value without uncertainty.
(f)
X
~
i
j is the upper bound of the triangular membership function definition interval of the ith Fuzzy input for the jth DMU, considering the average as the most probable value without uncertainty.
(g)
α is the value chosen for the α-level based approach, with α
∈ variation [
0
,
1
].
(h)
Ψ
i
o is the α coefficient in constraints linked to the ith Fuzzy input for DMU0.
(i)
ρ
j
0 is the α coefficient in constraints linked to the jth Fuzzy output for DMU0.
(j)
P
r
j is the α coefficient for constraints linked to the rth Fuzzy output of the jth DMU.
(k)
Ψ
i
j is the α coefficient for constraints linked to the ith Fuzzy input of the jth DMU.
(iii)
Decision variables:
(a)
u
r is the weight associated with the rth output.
(b)
v
i is the weight associated with the ith input.The Fuzzy-DEA-BCC modelα
L for an efficiency analysis of the DMUs in a pessimistic scenario is as follows:
(19)
max
⁡
E
j
=
∑
r
=
1
s
u
r
(
Y
~
r
0
-
ρ
r
0
·
α
)
+
c
0
subject
to:
∑
i
=
1
m
v
i
(
x
~
i
o
+
ψ
i
0
·
α
)
=
1
h
h
h
h
h
∑
r
=
1
s
u
r
(
Y
~
r
j
-
ρ
r
j
·
α
)
h
h
h
h
h
h
h
-
∑
i
=
1
m
v
i
(
x
~
i
j
+
ψ
i
j
·
α
)
+
c
0
≤
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
j
=
1,2
,
…
,
n
,
h
h
h
h
h
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
h
h
h
h
h
v
i
≥
0
,
i
=
1,2
,
…
,
m
,
h
h
h
h
h
c
0
free
.The Fuzzy-DEA-BCC modelα
U for an efficiency analysis of the DMUs in an optimistic scenario is as follows:
(20)
max
⁡
E
j
=
∑
r
=
1
s
u
r
(
y
~
r
0
+
ρ
r
0
·
α
)
+
c
0
subject
to:
∑
i
=
1
m
v
i
(
X
~
i
o
-
ψ
i
0
·
α
)
=
1
h
h
h
h
h
h
∑
r
=
1
s
u
r
(
y
~
r
j
+
ρ
r
j
·
α
)
-
∑
i
=
1
m
v
i
(
X
~
i
j
-
ψ
i
j
·
α
)
+
c
0
≤
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
j
=
1,2
,
…
,
n
,
h
h
h
h
h
h
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
h
h
h
h
h
h
v
i
≥
0
,
i
=
1,2
,
…
,
m
,
h
h
h
h
h
h
c
0
free
.Figure2 geometrically contemplates the position of the DEA models parameters via the triangular membership function. x
~
i
0 or x
~
i
j and y
~
r
0 or y
~
r
j correspond to the lower variation bound associated with inputs and outputs of the DMUs. X
~
i
0 or X
~
i
j and Y
~
r
0 or Y
~
r
j correspond to the upper variation bound of the DMU’s inputs and outputs. The mean value of the triangular membership function is associated with the values of these parameters in a scenario without uncertainty.Figure 2
Variation limits for the incorporation of uncertainty.


## 2.1. Deterministic Data Envelopment Analysis

Classic DEA models were introduced by Charnes et al. [15], which use constant returns of scale and were denominated as CCR models in homage to their creators [12]. These models were then extended by Banker et al. [12] to variable returns of scale, this time dubbed as BCC, again in reference to their originators.According to Cook and Seiford [17], DEA is a nonparametric methodology which comparatively measures the efficiency of each decision-making unit (DMU). Another fact that deserves highlighting is that, through its use, it is possible to avoid the problems created by incommensurability (different units of measurement) among the elements of the input and output matrices.In the original model from Charnes et al. [15], the input and output variable weights may be obtained from the Fractional Programming model solution, given by
(1)
E
j
=
∑
r
=
1
s
u
r
·
y
r
0
∑
i
=
1
m
v
i
·
x
i
0
subject to:
(2)
∑
r
=
1
s
u
r
y
r
j
∑
i
=
1
m
v
i
x
i
j
≤
1
,
j
=
1,2
,
…
,
n
,
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
v
i
≥
0
,
i
=
1,2
,
…
,
m
.With DMU0 to DMU being under evaluation, E
j is the relative efficiency of DMU0; y
r
0 and x
i
o are the input and output data for DMU0; j is the index for the DMU, j
=
1,2
,
…
n; r is the output index, with r
=
1,2
,
…
,
s; i is the input index, i
=
1,2
,
…
,
m; y
r
j is the value of rth output for jth DMU; x
i
j is the value of ith input for jth DMU; u
r is the weight associated with rth output; v
i is the weight associated with ith input.It is observed that ifE
j
=
1, DMU0 is going to be efficient when compared to the other units considered in the model, and, if 0
<
E
j
<
1, this DMU is deemed inefficient.This model (1)-(2) is not linear, thus having multiple solutions; nonetheless, it can be linearized, generating a DEA-CCD model, given by the following equation:
(3)
max
⁡
E
j
=
∑
r
=
1
s
u
r
·
y
r
0
subject
to:
∑
i
=
1
m
v
i
x
i
o
=
1
h
h
h
h
h
∑
r
=
1
s
u
r
y
r
j
-
∑
i
=
1
m
v
i
x
i
j
≤
0
,
j
=
1,2
,
…
,
n
,
h
h
h
h
h
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
h
h
h
h
h
v
i
≥
0
,
i
=
1,2
,
…
,
m
.DEA-BCC models, with variable returns of scale, can be expressed by the following equations:(4)
max
⁡
E
j
=
∑
r
=
1
s
u
r
·
y
r
0
+
c
0
(5)
subject
to:
∑
i
=
1
m
v
i
x
i
o
=
1
h
h
h
h
h
h
(6)
h
h
h
h
h
h
h
h
h
h
∑
r
=
1
s
u
r
y
r
j
-
∑
i
=
1
m
v
i
x
i
j
+
c
0
≤
0
,
h
h
h
h
h
h
h
h
h
h
h
j
=
1,2
,
…
,
n
.
(7)
h
h
h
h
h
h
h
h
h
u
r
≥
0
,
r
=
1,2
,
…
,
s
.
(8)
h
h
h
h
h
h
h
h
h
v
i
≥
0
,
i
=
1,2
,
…
,
m
.
(9)
c
0
free
.
hhhhA main difference between DEA-BCC and DEA-CCR models is the addition ofc
0 (free) variable which indicates variable returns of scale. Banker et al. [12] affirm that a DMU considered efficient in a BCC model will also be considered efficient in the CCR model; the inverse, however, is not necessarily true.According to Cooper et al. [16], in order to provide a suitable discrimination of the DMUs in traditional DEA models, the following must be verified: number of DMUs ≥ maximum {(product of the number of inputs and outputs), 3·(number of inputs + number of outputs)}.Classic DEA models consider DMUs withE
j
=
1 as efficient and DMUs with 0
≤
E
j
<
1 as inefficient. It is within the realm of possibility that multiple DMUs could be efficient; that is, this would not be a good discrimination of DMUs. In order to deal with this limitation, Andersen and Petersen [19] proposed the concept of superefficiency in order to help differentiate DMUs which present E
j
=
1.For the superefficiency evaluation to be employed in DEA-BCC models, constraint (6) must be removed in order for the DMU under analysis to attain greater scores than 1, being that the DMUs considered inefficient in traditional evaluation models will continue to be inefficient, but those which posed scores equal to one will be able to demonstrate scores above 1, thus enabling the elaboration of a ranking.


## 2.2. Fuzzy-DEA

Hatami-Marbini et al. [21] undertook an expansive literature review about the Fuzzy Theory combined with DEA models. The authors’ motivations were based on the fact that, in general, the estimation of input and output DMU values in real problems is difficult, which could generate efficiency values with a low level of reliability; one approach to deal with the aspects of uncertainty involves the adoption of Fuzzy Theory concepts.Kao and Liu [13] assert that measuring a DMU’s efficiency is a difficult task which involves complex economic variables such as interest rates and tax rates and employment levels and demand. According to these authors, efficiency measurement becomes even more difficult when analyzing multiple inputs and outputs.In this context, Wen et al. [22] comment that DMUs can be classified in two categories: efficient and inefficient. However, the incorporation of uncertainty as an error in measurement of inputs and outputs can make the calculation of efficiency more reliable and robust.Fuzzy Theory has come to be used with the objective of modeling uncertainty in DEA models [23]. Fuzzy-DEA models are based on Fuzzy Linear Programming, with DEA-BCC models with Fuzzy coefficients having special emphasis in this paper (Fuzzy-DEA-BCC), proposed by Kao and Liu [13]. Consider the following:
(10)
E
~
j
=
Max
⁡
∑
r
=
1
s
u
r
y
~
r
0
+
c
0
subject to:
(11)
∑
i
=
1
m
x
~
i
0
·
v
i
=
1
(12)
∑
r
=
1
s
u
r
y
~
r
j
-
∑
i
=
1
m
v
i
x
~
i
j
+
c
0
≤
0
,
j
=
1,2
,
…
,
n
,
(13)
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
(14)
v
i
≥
0
,
i
=
1,2
,
…
,
m
.The value of the objective function (10) may have a greater value than 1 due to constraints (11)-(12) which involve Fuzzy parameters and are solved using probability [23]. With the incorporation of Fuzzy coefficients, DEA-BCC models cannot be resolved using traditional linear programming (LP) techniques. Hatami-Marbini et al. [21] list and describe the following main approaches which deal with Fuzzy-DEA:(i)
theα-level based approach;
(ii)
the tolerance approach;
(iii)
the Fuzzy ranking approach;
(iv)
the possibility approach.In this study, the approach based on theα-level was adopted and is described afterwards. This α-level application is the most common for Fuzzy-DEA models according to Hatami-Marbini et al. [21]. To apply this method, the idea is to convert the Fuzzy-DEA model into a pair of Parametric Programming Problems to find upper and lower boundaries for the DMU efficiency score membership functions [23].For the purposes of this research, triangular membership functions were utilized. According to Liang and Wang [24], such functions well represent human expertise in adequately judging the behavior of common variables in a range of practical situations. Along these same lines of thought, Aouni et al. [25] show multiple applications for Fuzzy triangular numbers which validate and justify the adoption of such a method in conjunction with goal programming (GP) models. Another justification arises from the fact that it is a linear function, thus easing the optimization process through traditional LP means [21].To illustrate the process of Fuzzy-DEA modeling, the example in Kao and Liu’s [13] paper was used, in which they dealt with a model with a variable return of scale. Figure 1 shows the DEA-CCR and DEA-BCC models with four DMUs, denominated as A, B, C, and D, with only a single input and a single output for each DMU and with the output of DMU B being Fuzzy.Figure 1
Production frontier when DMU B output is Fuzzy. Source: adapted from Kao and Liu [13].The input values for DMUs A, B, C, and D are 10, 20, 30, and 50, respectively. With these inputs, outputs 5, 9, and 15 from DMUs A, C, and D are also produced, respectively. These inputs and outputs are associated with points A = (10; 5), C = (30; 9), and D = (50; 15) in Figure1. For DMU B, the output associated with input value 20 is represented by the trapezoidal Fuzzy number (5; 6; 8; 9), which is also illustrated in Figure 1 by points (20; 5), (20; 6), (20; 7.5), (20; 8), and (20; 9).Based on Figure1 and according to the DEA-BCC model [13], when the output of DMU B is less than or equal to 7.5, the production boundary is defined by the line segments connecting point (10; 0) to A = (10; 5) and from A to D = (50; 15). The efficiency scores of DMUs A, C, and D are 1, 0.9, and 1, respectively; the efficiency of DMU B is between 5/7.5 = 0.67 and 7.5/7.5 = 1, depending on its output.If an increase in the output associated with DMU B went from 7.5 to 9, the production boundary would be represented by the line segments linking points (10; 0) to A from A to B = (20; 9) and from B to D = (50; 15). In this case, DMUs A, B, and D have an efficient output bound, while the efficiency of C will be, according to output orientation, between values of 9/10= 0.9 and 9/11 = 0.82.In this example, regardless of the output value of DMU B, DMUs A and D are always efficient. In other words, with the combination (generation of scenarios) of each Fuzzy output for DMU B, the effects of uncertainty on DMUs A and D are always the same and present 100% efficiency.Conducting an analysis under the assumption of constant return of scale (DEA-CCR model), as seen in Figure1, the production boundary will be the solid line which links the origin point (0; 0) with point A = (10; 6). In this case, only DMU A will always be efficient (with an efficiency score of 100%), regardless of the output value of DMU B. However, the efficiency of DMU B will vary from 1/2 = 0.5 to 9/10 = 0.9, and DMUs C and D will have an efficiency of 3/5 = 0.6; that is, there are not effects of the Fuzzy output of DMU B over the efficiencies of DMUs C and D.As a further example, the membership function linked to the output of DMU B may be illustrated, given that it is a trapezoidal function, shown by the following equation:(15)
μ
Y
~
B
(
y
)
=
{
y
-
5
,
5
≤
y
≤
6
,
1
,
6
≤
y
≤
7
,
9
-
y
,
8
≤
y
≤
9
.Theα-level base is defined as the interval [
5
+
α
,
9
-
α
], being (
E
B
)
α
L
e
(
E
B
)
α
U, respectively. The upper and lower bounds for the variation (uncertainty) of α-level, which are associated with a scenario generation, vary between pessimistic (α
=
0) and optimistic (α
=
1) for an efficiency analysis.As previously stated, theα-level based approach was adopted which, according to Hatami-Marbini et al. [21], is the most popular Fuzzy-DEA model with many referenced publications [13]. The value of α
∈
[
0
,
1
] allows the generation of scenarios, that is, different efficiency values respecting the variation range determined by the membership function. In such models, X
~
i
j and Y
~
r
j are, respectively, the Fuzzy parameters for ith input and rth output of the jth DMU. These values are approximately known and can be represented by Fuzzy sets through means of membership functions μ
X
~
i
j and μ
Y
~
r
j. Thus, Fuzzy-DEA models can be formed using S
(
X
~
i
j
) and S
(
Y
~
r
j
) which are the values of these parameters for a given S scenario. With the α-level which generates a set of scenarios for Y
~
r
j and X
~
i
j, the formulas, as defined by Kao and Liu [13], are as follows:
(16)
(
X
i
j
)
α
=
{
x
i
j
∈
S
(
X
~
i
j
)
∣
μ
X
~
i
j
(
x
i
j
)
≥
α
}
,
h
h
h
h
h
h
h
h
i
i
i
=
1,2
,
…
m
,
j
=
1,2
,
…
,
n
,
(
Y
r
j
)
α
=
{
y
r
j
∈
S
(
Y
~
r
j
)
∣
μ
Y
~
r
j
(
y
r
j
)
≥
α
}
,
h
h
h
h
h
h
h
h
i
r
=
1,2
,
…
,
s
,
j
=
1,2
,
…
,
n
.Fuzzy-DEA can be transformed into a family of DEA models with different levels of uncertainty:{
(
X
i
j
)
α
∣
0
≤
α
≤
1
} and {
(
Y
r
j
)
α
∣
0
≤
α
≤
1
}. The results for each scenario identify the uncertainty variation range in the model’s input and output data [13, 26–28], with the α-level being defined by the following equation, according to Kao and Liu [13]:
(17)
(
X
i
j
)
α
=
[
{
min
⁡
x
i
j
x
i
j
∈
S
(
X
~
i
j
)
∣
μ
X
~
i
j
(
x
i
j
)
≥
α
}
,
h
h
h
h
max
⁡
x
i
j
{
x
i
j
∈
S
(
X
~
i
j
)
∣
μ
X
~
i
j
(
x
i
j
)
≥
α
}
]
,
(
Y
r
j
)
α
=
[
{
min
⁡
y
r
j
y
r
j
∈
S
(
Y
~
r
j
)
∣
μ
Y
~
r
j
(
y
r
j
)
≥
α
}
,
h
h
h
h
max
⁡
y
r
j
{
y
r
j
∈
S
(
Y
~
r
j
)
∣
μ
Y
~
r
j
(
y
r
j
)
≥
α
}
]
.Kao and Liu [13], based on Yager [29], Zadeh [28], and Zimmermann [30], established that a membership function which defines the efficiency of DMU j can be expressed by
(18)
μ
E
~
j
(
z
)
=
sup
⁡
x
,
y
min
⁡
=
{
μ
X
~
i
j
(
x
i
j
)
,
μ
Y
~
r
j
(
y
r
j
)
,
∀
i
,
r
,
j
∣
z
=
E
j
(
x
,
y
)
}
,
with E
j
(
x
,
y
) obtained by (4)–(9). The approach for constructing the membership function μ
E
~
j proposed in this study adopted an α-level μ
E
~
j as being the value of efficiency in each scenario E
~
j obtained by (10)–(14). For further details, Kao and Liu [13], Hatami-Marbini et al. [21], and Kao and Lin [31] are recommended reading materials.Based on the models developed by Banker et al. [12], in accordance with (4)–(9), and Kao and Liu [13], in accordance with (10)–(14), the Fuzzy-DEA-BCC model was developed. Afterwards, the indices, parameters, auxiliary variables and decision variables, objective functions, and model constraints are proposed, considering DMU0 as the one under analysis.(i)
Indices:
(a)
j is the DMU index.
(b)
r is the output index.
(c)
i is the input index.
(ii)
Parameters:
(a)
y
~
r
0 and x
~
i
0 are, respectively, the upper and lower bounds in the definition intervals of the triangular membership function for the rth Fuzzy output and the ith Fuzzy input for DMU0, considering the average as the most probable value without uncertainty.
(b)
Y
~
r
0 and X
~
i
0 are, respectively, the upper and lower bounds in the definition intervals of the triangular membership function for the rth Fuzzy output and the ith Fuzzy input for DMU0, considering the average as the most probable value without uncertainty.
(c)
y
~
r
j is the lower bound for the definition interval of the triangular membership function of the rth Fuzzy output for the jth DMU, considering the average as the most probable value without uncertainty.
(d)
Y
~
r
j is the upper bound of the triangular membership function definition interval of the rth Fuzzy output for the jth DMU, considering the average as the most probable value without uncertainty.
(e)
x
~
i
j is the lower bound of the triangular membership function definition interval of the rth Fuzzy output for the jth DMU, considering the average as the most probable value without uncertainty.
(f)
X
~
i
j is the upper bound of the triangular membership function definition interval of the ith Fuzzy input for the jth DMU, considering the average as the most probable value without uncertainty.
(g)
α is the value chosen for the α-level based approach, with α
∈ variation [
0
,
1
].
(h)
Ψ
i
o is the α coefficient in constraints linked to the ith Fuzzy input for DMU0.
(i)
ρ
j
0 is the α coefficient in constraints linked to the jth Fuzzy output for DMU0.
(j)
P
r
j is the α coefficient for constraints linked to the rth Fuzzy output of the jth DMU.
(k)
Ψ
i
j is the α coefficient for constraints linked to the ith Fuzzy input of the jth DMU.
(iii)
Decision variables:
(a)
u
r is the weight associated with the rth output.
(b)
v
i is the weight associated with the ith input.The Fuzzy-DEA-BCC modelα
L for an efficiency analysis of the DMUs in a pessimistic scenario is as follows:
(19)
max
⁡
E
j
=
∑
r
=
1
s
u
r
(
Y
~
r
0
-
ρ
r
0
·
α
)
+
c
0
subject
to:
∑
i
=
1
m
v
i
(
x
~
i
o
+
ψ
i
0
·
α
)
=
1
h
h
h
h
h
∑
r
=
1
s
u
r
(
Y
~
r
j
-
ρ
r
j
·
α
)
h
h
h
h
h
h
h
-
∑
i
=
1
m
v
i
(
x
~
i
j
+
ψ
i
j
·
α
)
+
c
0
≤
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
j
=
1,2
,
…
,
n
,
h
h
h
h
h
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
h
h
h
h
h
v
i
≥
0
,
i
=
1,2
,
…
,
m
,
h
h
h
h
h
c
0
free
.The Fuzzy-DEA-BCC modelα
U for an efficiency analysis of the DMUs in an optimistic scenario is as follows:
(20)
max
⁡
E
j
=
∑
r
=
1
s
u
r
(
y
~
r
0
+
ρ
r
0
·
α
)
+
c
0
subject
to:
∑
i
=
1
m
v
i
(
X
~
i
o
-
ψ
i
0
·
α
)
=
1
h
h
h
h
h
h
∑
r
=
1
s
u
r
(
y
~
r
j
+
ρ
r
j
·
α
)
-
∑
i
=
1
m
v
i
(
X
~
i
j
-
ψ
i
j
·
α
)
+
c
0
≤
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
j
=
1,2
,
…
,
n
,
h
h
h
h
h
h
u
r
≥
0
,
r
=
1,2
,
…
,
s
,
h
h
h
h
h
h
v
i
≥
0
,
i
=
1,2
,
…
,
m
,
h
h
h
h
h
h
c
0
free
.Figure2 geometrically contemplates the position of the DEA models parameters via the triangular membership function. x
~
i
0 or x
~
i
j and y
~
r
0 or y
~
r
j correspond to the lower variation bound associated with inputs and outputs of the DMUs. X
~
i
0 or X
~
i
j and Y
~
r
0 or Y
~
r
j correspond to the upper variation bound of the DMU’s inputs and outputs. The mean value of the triangular membership function is associated with the values of these parameters in a scenario without uncertainty.Figure 2
Variation limits for the incorporation of uncertainty.


## 3. Problem Description and Modeling

### 3.1. Proposal of Integrating Fuzzy-DEA with Simulation Optimization Problems

This paper’s proposal of integrating Fuzzy-DEA-BCC models in simulation optimization problems is based on the following four techniques:(i)
discrete-event simulation to represent the real system to be optimized and conduct scenario simulation and data collection;
(ii)
the Taguchi orthogonal arrays to generate experimental matrices and define simulation runs to be executed in order to represent the search space;
(iii)
Fuzzy-DEA-BCC to analyze the efficiency of each generated scenario, taking into account the uncertainty present in each of them and classifying them in terms of the most efficient DMUs;
(iv)
an optimization procedure via simulation which can perform searches for optimal solutions.The use of optimization assumes that the simulation model is constructed, verified, and validated, thus assuring that the model adequately simulates the reality of the phenomenon under study. It is also suggested that the response variables are either discrete or integers. The steps for their use are presented in Figure3.Figure 3
Procedure proposed for an optimization via simulation utilizing Fuzzy-DEA-BCC and superefficiency.The application phases for the proposed procedure are described below.Step  1. Define the simulation model decision variables (x
1
,
x
2
,
x
3
,
…
,
x
n) and the variation ranges for each variable (lower level ≤ x
i ≤ upper level, com 1
≤
i
≤
n).Step  2. Determine the output variables (one or more) to be optimized (y
1
,
y
2
,
y
3
,
…
,
y
m).Step  3. Select the Taguchi orthogonal arrays in function of the number of decision variables and their variation limits. This selection must obey the fundamental rule established for the minimum number of DMUs to be analyzed through Fuzzy-DEA-BCC analysis [16]. After the arrays selection, generate an experimental matrix which represents the most diversified solution region as possible and explore all levels of each decision variable, if possible.Step  4. Execute the experiments in a discrete-event simulator and store the maximum, minimum, and mean values for each output variable to be optimized for analysis.Step 5. Fuzzy analysis is as follows. Based on the experiments carried out in the previous step, insert the minimum, maximum, and mean values for each experiment in the triangular membership function. The choice of this membership function is based on comments by Liang and Wang [24] who justify the use of Fuzzy triangular functions, given that they suitably mirror human judgments. Thus, according to these authors, the advantage of using triangular Fuzzy number lies in the fact that not only human expertise can be suitably represented, but also the model additionally accounts for uncertainty in the involved data and parameters. The optimization models dealing with uncertainty which contain present and future information cannot be perfectly accounted for and should be considered as uncertain [32]. This study utilized the α-level based approach given that, according to Wang and Liang [32], it is the most popular Fuzzy-DEA model due to the great number of publications using Fuzzy-DEA in current scientific literature [21]. For more information on these approaches, please see Hatami-Marbini et al. [21].Step  6. Determination of efficiency for each scenario by means of Fuzzy-DEA-BCC for the simulated results is as follows. Upon the definition of maximum, minimum, and mean values linked to each membership function associated with each analyzed condition, the Fuzzy-DEA-BCC model was applied using the α-level based approach, varying the α value (0,0.1,0.2
,
…
,
1), carrying out 11 scenarios with different superefficiency values for both pessimistic and optimistic scenarios.Step  7. Ranking the most efficient DMUs based on the concept of superefficiency is as follows. Given that there are 11 pessimistic scenarios and 11 optimistic scenarios, it was necessary to lump them together into a global scenario. The adopted method was to extract the geometric average of each scenario, thus arriving at a global scenario. Thus, the geometric averages of all linked scenarios were obtained and ranked from greatest to lowest. The first and second positions in the rank were chosen to reduce the range for each decision variable and, in turn, carry out the optimization simulation.Step  8. Reducing the variation range of each decision variable is as follows. The first and second positions of the rank were chosen to carry out the range reduction of each variable, thus taking out those variables with equal values for both DMUs. This was the value adopted for that variable.Step  9. Optimize the simulation model with a new variation range for each decision variable.Step  10. Analyze the results of the optimization and make decisions based on the results found.In order to exemplify the application of the proposed procedure, real situations involving the optimization of two simulation models are presented in the following sections. The utilized models were previously verified and validated, thus proving to be apt for optimization of the simulation.


## 3.1. Proposal of Integrating Fuzzy-DEA with Simulation Optimization Problems

This paper’s proposal of integrating Fuzzy-DEA-BCC models in simulation optimization problems is based on the following four techniques:(i)
discrete-event simulation to represent the real system to be optimized and conduct scenario simulation and data collection;
(ii)
the Taguchi orthogonal arrays to generate experimental matrices and define simulation runs to be executed in order to represent the search space;
(iii)
Fuzzy-DEA-BCC to analyze the efficiency of each generated scenario, taking into account the uncertainty present in each of them and classifying them in terms of the most efficient DMUs;
(iv)
an optimization procedure via simulation which can perform searches for optimal solutions.The use of optimization assumes that the simulation model is constructed, verified, and validated, thus assuring that the model adequately simulates the reality of the phenomenon under study. It is also suggested that the response variables are either discrete or integers. The steps for their use are presented in Figure3.Figure 3
Procedure proposed for an optimization via simulation utilizing Fuzzy-DEA-BCC and superefficiency.The application phases for the proposed procedure are described below.Step  1. Define the simulation model decision variables (x
1
,
x
2
,
x
3
,
…
,
x
n) and the variation ranges for each variable (lower level ≤ x
i ≤ upper level, com 1
≤
i
≤
n).Step  2. Determine the output variables (one or more) to be optimized (y
1
,
y
2
,
y
3
,
…
,
y
m).Step  3. Select the Taguchi orthogonal arrays in function of the number of decision variables and their variation limits. This selection must obey the fundamental rule established for the minimum number of DMUs to be analyzed through Fuzzy-DEA-BCC analysis [16]. After the arrays selection, generate an experimental matrix which represents the most diversified solution region as possible and explore all levels of each decision variable, if possible.Step  4. Execute the experiments in a discrete-event simulator and store the maximum, minimum, and mean values for each output variable to be optimized for analysis.Step 5. Fuzzy analysis is as follows. Based on the experiments carried out in the previous step, insert the minimum, maximum, and mean values for each experiment in the triangular membership function. The choice of this membership function is based on comments by Liang and Wang [24] who justify the use of Fuzzy triangular functions, given that they suitably mirror human judgments. Thus, according to these authors, the advantage of using triangular Fuzzy number lies in the fact that not only human expertise can be suitably represented, but also the model additionally accounts for uncertainty in the involved data and parameters. The optimization models dealing with uncertainty which contain present and future information cannot be perfectly accounted for and should be considered as uncertain [32]. This study utilized the α-level based approach given that, according to Wang and Liang [32], it is the most popular Fuzzy-DEA model due to the great number of publications using Fuzzy-DEA in current scientific literature [21]. For more information on these approaches, please see Hatami-Marbini et al. [21].Step  6. Determination of efficiency for each scenario by means of Fuzzy-DEA-BCC for the simulated results is as follows. Upon the definition of maximum, minimum, and mean values linked to each membership function associated with each analyzed condition, the Fuzzy-DEA-BCC model was applied using the α-level based approach, varying the α value (0,0.1,0.2
,
…
,
1), carrying out 11 scenarios with different superefficiency values for both pessimistic and optimistic scenarios.Step  7. Ranking the most efficient DMUs based on the concept of superefficiency is as follows. Given that there are 11 pessimistic scenarios and 11 optimistic scenarios, it was necessary to lump them together into a global scenario. The adopted method was to extract the geometric average of each scenario, thus arriving at a global scenario. Thus, the geometric averages of all linked scenarios were obtained and ranked from greatest to lowest. The first and second positions in the rank were chosen to reduce the range for each decision variable and, in turn, carry out the optimization simulation.Step  8. Reducing the variation range of each decision variable is as follows. The first and second positions of the rank were chosen to carry out the range reduction of each variable, thus taking out those variables with equal values for both DMUs. This was the value adopted for that variable.Step  9. Optimize the simulation model with a new variation range for each decision variable.Step  10. Analyze the results of the optimization and make decisions based on the results found.In order to exemplify the application of the proposed procedure, real situations involving the optimization of two simulation models are presented in the following sections. The utilized models were previously verified and validated, thus proving to be apt for optimization of the simulation.


## 4. Application and Data Analysis

### 4.1. Monoobjective Case

The modeled situation corresponds to a quality control cell in a fiber-optic transponder company. The cell is responsible for a series of tests which will ensure the approval or failure of the equipment produced at the plant. All verification and validation phases were appropriately undertaken, thus providing a consistent model.For this study object, the decision variables were defined as the number of operators responsible for carrying out quality control tests, denoted by(
x
1
,
x
2
,
x
3
), and the number of pieces of equipment for test types 1, 2, and 3, denoted by (
x
4
,
x
5
,
x
6
). All variables were defined as integers, with a lower bound equal to one and an upper one equal to five. Table 1 presents this information.Table 1
Decision variables, types, and limits for the study object.
Variable
Type
Lower bound
Upper bound
x
1
Number of operators, type 1
Integer
1
5
x
2
Number of operators, type 2
Integer
1
5
x
3
Number of operators, type 3
Integer
1
5
x
4
Number of pieces of equipment for test, type 1
Integer
1
5
x
5
Number of pieces of equipment for test, type 2
Integer
1
5
x
6
Number of pieces of equipment for test, type 3
Integer
1
5The optimization objective was to find the combination of decision variables which would maximize the total number of inspected products in the quality control center, denoted by (y
1). For the problem in question, considering the number of decision variables (
6
) and their maximum variation (1–5), there are a total of 15.625 = (56) possible scenarios for the search space.Considering the quantity of decision variables, the variation of levels of each variable, and the rule for the minimum number of DMUs proposed by Cooper et al. [16], the orthogonal array L25 was chosen. Seeing that there are six input variables and one output variable, there would be a minimum of 21 DMUs (runs) by following the classic rule, justifying the use of array L25. With a defined array, the experimental matrix was generated and presented in Table 2.Table 2
Experimental matrix and results.
DMU
Input variables
Output variables
x
1
x
2
x
3
x
4
x
5
x
6
y
1 minimum
y
1 average
y
1 maximum
1
1
1
1
1
1
1
3,922
3,977
4,032
2
1
2
2
2
2
2
3,936
3,979
4,022
3
1
3
3
3
3
3
3,919
3,978
4,037
4
1
4
4
4
4
4
3,932
3,985
4,038
5
1
5
5
5
5
5
3,939
3,982
4,025
6
2
1
2
3
4
5
3,939
3,982
4,025
7
2
2
3
4
5
1
3,925
3,978
4,031
8
2
3
4
5
1
2
3,930
3,982
4,034
9
2
4
5
1
2
3
3,938
3,979
4,020
10
2
5
1
2
3
4
3,920
3,983
4,046
11
3
1
3
5
2
4
3,918
3,982
4,046
12
3
2
4
1
3
5
3,910
3,979
4,048
13
3
3
5
2
4
1
3,924
3,983
4,042
14
3
4
1
3
5
2
3,926
3,981
4,036
15
3
5
2
4
1
3
3,937
3,983
4,029
16
4
1
4
2
5
3
3,938
3,980
4,022
17
4
2
5
3
1
4
3,930
3,978
4,026
18
4
3
1
4
2
5
3,932
3,982
4,032
19
4
4
2
5
3
1
3,922
3,981
4,040
20
4
5
3
1
4
2
3,921
3,982
4,043
21
5
1
5
4
3
2
3,937
3,980
4,023
22
5
2
1
5
4
3
3,932
3,982
4,032
23
5
3
2
1
5
4
3,930
3,982
4,034
24
5
4
3
2
1
5
3,933
3,982
4,031
25
5
5
4
3
2
1
3,933
3,982
4,031In the following, 25 scenarios were simulated from array L25 with 30 replications of a month’s time of operation in quality control. The simulations were carried out on a computer with an Intel processor (Core 2 Duo) 1.58 GHZ, a 2 GB RAM, and a 64-bit Microsoft operational system platform.Data for each output variable were stored for the superefficiency calculations. Nearly 26 minutes were spent to process the 25 scenarios, considering all replications. Maximum, minimum, and average values were stored for 30 replications for each scenario with the goal of using them in the triangular membership function. Results for output variabley
1 are shown in Table 2.For the calculation of superefficiency related to each DMU with a Fuzzy-DEA-BCC model, the softwareThe General Algebraic Modeling (GAMS) [33] was used, version 22.8.1, using the solver CPLEX, version 11.0, adapted for the specific calculation.With these results, the superefficiency value of the DMU can be related to each scenario for the pessimistic and optimistic scenarios by means of anα variation, summing up to 11 scenarios. These values are presented in Tables 3 and 4.Table 3
Superefficiency under uncertainty for a pessimistic scenario.
Scenario
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11
α
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
DMU 1
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
DMU 2
1.750
1.742
1.733
1.715
1.691
1.625
1.535
1.407
1.253
1.079
1.000
DMU 3
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 4
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.207
6.000
6.714
2.833
DMU 5
2.000
2.000
2.000
2.000
2.000
2.000
2.000
1.000
1.000
1.000
1.000
DMU 6
5.000
5.000
5.000
5.000
5.000
5.000
5.000
4.000
2.929
2.462
1.446
DMU 7
1.014
1.012
1.010
1.007
1.004
1.000
1.000
1.000
1.000
1.000
1.000
DMU 8
1.023
1.033
1.044
1.057
1.074
1.094
1.121
1.156
1.202
1.268
1.316
DMU 9
2.222
2.178
2.133
2.089
2.044
1.976
1.883
1.676
1.528
1.330
1.052
DMU 10
1.000
1.000
1.000
1.000
1.000
1.000
1.047
1.161
1.306
1.480
2.000
DMU 11
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.064
1.294
1.571
DMU 12
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 13
1.000
1.000
1.000
1.000
1.000
1.000
1.067
1.156
1.280
1.475
2.000
DMU 14
1.000
1.000
1.000
1.000
1.000
1.021
1.049
1.085
1.128
1.181
1.200
DMU 15
1.933
2.037
2.208
2.625
3.111
3.500
3.750
4.167
3.182
2.500
2.000
DMU 16
1.333
1.329
1.323
1.316
1.307
1.296
1.282
1.263
1.233
1.161
1.074
DMU 17
1.222
1.207
1.190
1.168
1.143
1.111
1.071
1.018
1.000
1.000
1.000
DMU 18
1.308
1.310
1.313
1.315
1.319
1.323
1.327
1.332
1.328
1.296
1.250
DMU 19
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.063
1.139
1.224
1.318
DMU 20
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.011
1.167
1.288
1.357
DMU 21
1.318
1.320
1.322
1.325
1.329
1.333
1.351
1.367
1.367
1.346
1.211
DMU 22
1.263
1.271
1.280
1.290
1.303
1.319
1.341
1.352
1.367
1.389
1.400
DMU 23
1.222
1.246
1.275
1.311
1.356
1.415
1.494
1.611
1.676
1.711
1.400
DMU 24
1.231
1.238
1.245
1.255
1.267
1.281
1.300
1.325
1.360
1.430
1.500
DMU 25
1.683
1.698
1.714
1.732
1.758
1.800
1.882
2.026
1.769
1.442
1.333Table 4
Superefficiency under uncertainty for an optimistic scenario.
Scenario
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11
α
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
DMU 1
2.441
2.424
2.403
2.377
2.345
2.304
2.249
2.171
2.054
2.000
2.000
DMU 2
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 3
1.081
1.040
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 4
1.111
1.186
1.268
1.359
1.459
1.571
1.697
1.839
2.000
2.333
2.833
DMU 5
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 6
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.137
1.446
DMU 7
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 8
1.143
1.176
1.210
1.243
1.284
1.314
1.351
1.399
1.397
1.356
1.316
DMU 9
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.052
DMU 10
3.000
3.067
3.200
3.600
3.125
3.182
3.261
3.378
4.000
3.091
2.000
DMU 11
1.875
1.897
1.924
1.958
2.000
3.500
3.875
3.833
3.000
2.575
1.571
DMU 12
2.667
30.500
10.500
3.833
2.207
2.000
1.765
1.595
1.382
1.164
1.000
DMU 13
1.667
1.714
1.769
1.833
1.909
2.000
2.222
2.500
2.857
2.727
2.000
DMU 14
1.000
1.000
1.000
1.009
1.029
1.053
1.080
1.111
1.149
1.194
1.200
DMU 15
1.000
1.000
1.000
1.000
1.000
1.035
1.128
1.221
1.331
1.476
2.000
DMU 16
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.074
DMU 17
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 18
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.022
1.137
1.250
DMU 19
1.500
1.500
1.500
1.500
1.500
1.500
1.474
1.447
1.419
1.389
1.318
DMU 20
1.389
1.399
1.411
1.425
1.441
1.460
1.482
1.529
1.611
1.493
1.357
DMU 21
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.211
DMU 22
1.000
1.000
1.000
1.000
1.000
1.000
1.079
1.174
1.278
1.354
1.400
DMU 23
1.000
1.000
1.000
1.000
1.000
1.000
1.008
1.062
1.140
1.257
1.400
DMU 24
1.000
1.000
1.000
1.008
1.068
1.145
1.246
1.359
1.407
1.454
1.500
DMU 25
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.069
1.159
1.247
1.333For the calculation of superefficiency for each scenario, the average geometric value was extracted for both pessimistic and optimistic scenarios. The results found are presented in Table5. With these data, it was possible to calculate the average of each DMU and then rank them in function of their superefficiency values (see Table 5).Table 5
Global matrix for geometric averages for pessimistic and optimistic scenarios.
DMU
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11
Average superefficiency
Ranking
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1
2.210
2.202
2.192
2.180
2.166
2.147
2.121
2.084
2.027
2.000
2.000
2.121
1
2
1.323
1.320
1.317
1.310
1.300
1.275
1.239
1.186
1.119
1.039
1.000
1.221
16
3
1.040
1.020
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.005
24
4
1.054
1.089
1.126
1.166
1.208
1.254
1.303
1.489
3.464
3.958
2.833
1.813
5
5
1.414
1.414
1.414
1.414
1.414
1.414
1.414
1.000
1.000
1.000
1.000
1.264
11
6
2.236
2.236
2.236
2.236
2.236
2.236
2.236
2.000
1.711
1.673
1.446
2.044
2
7
1.007
1.006
1.005
1.004
1.002
1.000
1.000
1.000
1.000
1.000
1.000
1.002
25
8
1.081
1.102
1.124
1.147
1.174
1.199
1.231
1.272
1.296
1.311
1.316
1.205
18
9
1.491
1.476
1.461
1.445
1.430
1.406
1.372
1.295
1.236
1.153
1.052
1.347
10
10
1.732
1.751
1.789
1.897
1.768
1.784
1.847
1.980
2.286
2.139
2.000
1.907
4
11
1.369
1.377
1.387
1.399
1.414
1.871
1.969
1.958
1.786
1.825
1.571
1.630
7
12
1.633
5.523
3.240
1.958
1.486
1.414
1.328
1.263
1.176
1.079
1.000
1.918
3
13
1.291
1.309
1.330
1.354
1.382
1.414
1.540
1.700
1.912
2.006
2.000
1.567
8
14
1.000
1.000
1.000
1.005
1.015
1.037
1.064
1.098
1.138
1.188
1.200
1.068
22
15
1.390
1.427
1.486
1.620
1.764
1.903
2.057
2.255
2.058
1.921
2.000
1.807
6
16
1.155
1.153
1.150
1.147
1.143
1.139
1.132
1.124
1.110
1.078
1.074
1.128
21
17
1.106
1.099
1.091
1.081
1.069
1.054
1.035
1.009
1.000
1.000
1.000
1.049
23
18
1.144
1.145
1.146
1.147
1.148
1.150
1.152
1.154
1.165
1.214
1.250
1.165
19
19
1.225
1.225
1.225
1.225
1.225
1.225
1.214
1.240
1.271
1.304
1.318
1.245
14
20
1.179
1.183
1.188
1.194
1.200
1.208
1.218
1.244
1.371
1.387
1.357
1.248
13
21
1.148
1.149
1.150
1.151
1.153
1.155
1.162
1.169
1.169
1.160
1.211
1.162
20
22
1.124
1.127
1.131
1.136
1.142
1.149
1.203
1.260
1.321
1.372
1.400
1.215
17
23
1.106
1.116
1.129
1.145
1.164
1.189
1.227
1.308
1.382
1.466
1.400
1.239
15
24
1.109
1.112
1.116
1.125
1.163
1.211
1.273
1.342
1.383
1.442
1.500
1.252
12
25
1.297
1.303
1.309
1.316
1.326
1.342
1.372
1.472
1.432
1.341
1.333
1.349
9Through a superefficiency analysis, it was possible to rank the DMUs in order of efficiency. For the problem in question, DMU 1 is the most efficient, followed by DMU 6. Both are highlighted in Table5.With the identification of the two most efficient DMUs and based on the experimental matrix in Table2, a new interval can be identified for each decision variable in which a better set of solutions is expected. The new intervals for each decision variable are presented in Table 6. Variable x
2 stands out, as its value was already defined, being equal to x
2
=
1, reducing the number of decision variables to 5.Table 6
Decision variables, types, and new limits for the study object.
Variable
Type
Lower bound
Upper bound
x
1
Number of operators, type 1
Integer
1
2
x
3
Number of operators, type 3
Integer
1
2
x
4
Number of pieces of equipment for test, type 1
Integer
1
3
x
5
Number of pieces of equipment for test, type 2
Integer
1
4
x
6
Number of pieces of equipment for test, type 3
Integer
1
5With the reduction of the variation interval for each decision variable, the search space for the best solution was reduced from 15,625 to 240—a reduction of 98.4%.To confirm the efficiency of search space reduction, the optimizer SimRummer [34] was utilized for the simulation model optimization. SimRunner stands as popular optimization software [35] in operation and being sold in conjunction with the ProModel simulator. According to Kim et al. [36] and Banks et al. [4], this optimization software finds an optimal solution using as its search method a metaheuristic called Genetic Algorithm that as pointed by Ólafsson [37] mimics the process of natural selection.The optimizer was set to the same conditions and objectives; however, two optimizations were conducted for the same problem. One used a reduction in the variation interval (Table6), and the other one used the original problem (Table 1) conditions. The results found can be seen in Table 7.Table 7
Optimization results.
Decision variables
Solution
New range
Original range
x
1
2
5
x
2
1
3
x
3
1
2
x
4
2
2
x
5
1
1
x
6
3
5
Output variable
Responses
New range
Original range
y
1
3,986
3,987
Confidence interval fory
1 (95%)
(3,974–3,998)
(3,977–3,997)The responses presented by the optimizer for the decision variables were only equal for variablesx
4 and x
5. For the other variables, the optimizer arrived at different values. As for the solution of y
1, with search space reduction, the optimizer reached a value statistically equal to both optimized cases.The optimizer for the optimization problem with its range being reduced carried out 88 experiments before converging, which is equal to 36.67% of the experimental area (240 scenarios), taking a little more than 2.25 hours. For the problem with original range settings, the optimizer carried out 183 experiments, equal to slightly more than 1.17% of the total experimental area (15,625 scenarios), taking 4.8 hours.


### 4.2. Multiobjective Cases

The second study object represents a production cell in a Brazilian telecommunications company which produces fiber-optical equipment. This model, which is similar to the previous case, had already passed through the phases of verification and validation prior to being used for this paper in the simulation optimization technique. For the presented study object, the decision variables were defined as the number of pieces of inspection equipment, types 1 and 2, denoted by(
x
1
,
x
2
), and the number of employees, types 1, 2, and 3, denoted by (
x
3
,
x
4
,
x
5
), who carry out the activities presented in the model. The variables were defined as integers, with a lower bound of 1 and an upper bound of 5. Table 8 presents this information.Table 8
Decision variables, types, and limits for the second study object.
Variable
Type
Lower bound
Upper bound
x
1
Number of pieces of inspection equipment, type 1
Integer
1
5
x
2
Number of pieces of inspection equipment, type 2
Integer
1
5
x
3
Number of employees, type 1
Integer
1
5
x
4
Number of employees, type 2
Integer
1
5
x
5
Number of employees, type 3
Integer
1
5The optimization objective was to find the combination of decision variables which would maximize throughput (y
1) and production cell profit (y
2). It can be seen for the problem in question that, considering the number of variables (
5
) and their maximum variation (1–5), there are a total of 3,125 = (55) possible scenarios for the search space for the best configuration.In order to meet the fundamental rule for the number of DMUs, an L25 orthogonal array was used. An experimental matrix was generated which is presented in Table9. The scenarios for the experimental matrix were simulated in ProModel.Table 9
Experimental matrix and results.
DMU
Input variables
Output variables
x
1
x
2
x
3
x
4
x
5
y
1 minimum
y
1 average
y
1 maximum
y
2 minimum
y
2 average
y
2 maximum
1
1
1
1
1
1
643
863
1,083
381,096
682,082
983,067
2
1
2
2
2
2
1,295
1,395
1,495
826,452
1,097,340
1,368,228
3
1
3
3
3
3
1,188
1,381
1,574
852,107
1,080,373
1,308,639
4
1
4
4
4
4
1,235
1,399
1,563
726,952
1,091,904
1,456,857
5
1
5
5
5
5
1,148
1,403
1,658
690,053
1,088,349
1,486,646
6
2
1
2
3
4
1,064
1,286
1,508
748,800
1,007,407
1,266,015
7
2
2
3
4
5
1,219
1,383
1,547
842,802
1,067,833
1,292,864
8
2
3
4
5
1
1,175
1,391
1,607
783,439
1,081,793
1,380,148
9
2
4
5
1
2
573
868
1,163
465,318
681,706
898,094
10
2
5
1
2
3
895
1,140
1,385
515,279
896,302
1,277,326
11
3
1
3
5
2
1,153
1,404
1,655
862,630
1,095,382
1,328,135
12
3
2
4
1
3
736
847
958
425,035
660,137
895,239
13
3
3
5
2
4
1,129
1,383
1,637
804,861
1,077,000
1,349,140
14
3
4
1
3
5
894
1,167
1,440
532,523
911,184
1,289,846
15
3
5
2
4
1
1,275
1,397
1,519
840,255
1,096,732
1,353,210
16
4
1
4
2
5
1,178
1,396
1,614
704,903
1,084,938
1,464,973
17
4
2
5
3
1
1,178
1,409
1,640
710,034
1,101,006
1,491,979
18
4
3
1
4
2
1,013
1,123
1,233
516,328
871,771
1,227,214
19
4
4
2
5
3
1,215
1,385
1,555
727,103
1,071,027
1,414,953
20
4
5
3
1
4
586
828
1,070
426,589
638,185
849,782
21
5
1
5
4
3
1,264
1,406
1,548
720,561
1,093,357
1,466,152
22
5
2
1
5
4
974
1,157
1,340
582,466
888,310
1,194,154
23
5
3
2
1
5
561
829
1,097
406,997
635,928
864,860
24
5
4
3
2
1
1,109
1,301
1,493
738,318
1,014,712
1,291,106
25
5
5
4
3
2
1,214
1,379
1,544
834,734
1,068,090
1,301,447Thirty replications were simulated for each scenario, which corresponded to a month of operation in the production cell, and the data for each output variable were stored for the superefficiency calculation. Simulation of the 25 scenarios and their 30 replications took slightly more than 16.25 minutes. The results (minimum, maximum, and mean values) for the outputsy
1 and y
2 for each DMU are shown in Table 9.For the calculation of superefficiency related to each DMU for the Fuzzy-DEA-BCC model, the same procedure was employed. With the output variable results, each scenario can be related to each DMU superefficiency value for pessimistic and optimistic scenarios. In the following, the scenarios are lumped into a geometric average. The final result is presented in Table10. Thus, the average can be calculated for each DMU, making it possible to rank them in terms of superefficiency.Table 10
Global matrix for geometric average for pessimistic and optimistic scenarios.
DMU
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11
Average superefficiency
Ranking
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
4
2
3.803
3.661
3.512
3.360
3.203
3.043
3.006
2.915
2.758
2.700
2.792
3.159
1
3
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
18
4
1.019
1.022
1.024
1.028
1.036
1.055
1.168
1.085
1.051
1.021
1.000
1.046
16
5
1.008
1.010
1.012
1.016
1.067
1.285
1.432
1.547
1.591
1.674
2.200
1.349
8
6
1.281
1.278
1.269
1.258
1.248
1.238
1.228
1.219
1.209
1.200
1.191
1.238
11
7
0.805
0.808
0.811
0.811
0.810
0.811
0.814
0.815
0.813
0.812
0.812
0.811
24
8
1.270
1.271
1.271
1.271
1.271
1.271
1.271
1.271
1.271
1.271
1.271
1.271
10
9
1.016
1.015
1.014
1.013
1.012
1.010
1.009
1.009
1.009
1.009
1.009
1.011
17
10
1.030
1.046
1.066
1.087
1.108
1.130
1.153
1.176
1.189
1.188
1.188
1.124
15
11
1.952
2.000
2.069
2.178
2.353
2.691
7.030
1.683
1.665
1.647
1.431
2.427
3
12
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
19
13
0.977
0.982
0.982
0.983
0.984
0.984
0.985
0.986
0.987
0.987
0.988
0.984
22
14
1.652
1.614
1.574
1.535
1.495
1.455
1.414
1.316
1.213
1.126
1.086
1.407
7
15
1.539
1.555
1.566
1.582
1.614
1.649
2.709
1.953
1.854
1.655
1.578
1.750
5
16
1.448
1.461
1.481
1.500
1.496
1.493
1.494
1.493
1.491
1.491
1.490
1.485
6
17
1.694
1.704
1.719
1.739
3.999
2.856
2.873
3.015
3.138
2.808
2.009
2.505
2
18
1.106
1.114
1.123
1.132
1.142
1.152
1.164
1.169
1.173
1.176
1.176
1.148
14
19
0.847
0.858
0.868
0.877
0.887
0.900
0.902
0.919
0.961
0.971
0.975
0.906
23
20
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
20
21
1.313
1.316
1.318
1.318
1.315
1.316
1.318
1.322
1.337
1.352
1.400
1.330
9
22
1.123
1.126
1.138
1.150
1.164
1.179
1.195
1.213
1.232
1.240
1.240
1.182
13
23
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
21
24
1.235
1.231
1.227
1.224
1.220
1.216
1.213
1.210
1.207
1.205
1.205
1.218
12
25
0.750
0.753
0.756
0.755
0.757
0.757
0.757
0.758
0.758
0.758
0.758
0.756
25Based on the superefficiency value rankings (Table3) in this study, DMU 2 proved to be the most efficient followed by DMU 17. Both are highlighted in Table 10.Once these two DMUs are identified as the most efficient, a new interval for each decision variable can be redefined. In doing so, the optimization search space is reduced. Table9 shows that variable x
2 presented the same value for the most efficient DMUs, thus causing its value to be defined as x
2
=
2, reducing the number of variables from five to four. The new intervals for the other decision variables are presented in Table 11.Table 11
Decision variables, types, and new limits for the second study object.
Variable
Type
Lower bound
Upper bound
x
1
Number of pieces of inspection equipment, type 1
Integer
1
4
x
3
Number of employees, type 1
Integer
2
5
x
4
Number of employees, type 2
Integer
2
3
x
5
Number of employees, type 3
Integer
2
1By reducing the variation interval for each decision variable, the search space for the optimal solution was reduced from 3,125 to 64—a reduction of nearly 98%.To test the efficiency and robustness of the new search space, SimRunner was set to carry out the simulation model optimization, aiming to maximize total cell production(
y
1
) and total profit (
y
2
), for both the original variation range (Table 8) and the reduced variation range (Table 11). Results are seen in Table 12.Table 12
Optimization results.
Decision variables
Solution
New range
Original range
x
1
1
1
x
2
2
5
x
3
2
4
x
4
3
2
x
5
1
2
Output variable
Responses
New range
Original range
y
1
1,402
1,403
Confidence Interval fory
1 (95%)
(1,395–1,409)
(1,397–1,410)
y
2
1,103,794
1,103,595
Confidence Interval fory
2 (95%)
(1,096,972–1,110,617)
(1,095,465–1,111,727)The responses presented by the optimizer for the decision variables were only equal for variablex
1
=
1. For the other decision variables, the optimizer arrived at different values. In the case of the original variation range, the responses indicated the need to hire more employees and purchase more equipment, except for variable x
4. As for the solutions of y
1 and y
2, the results were statistically equal for a confidence level of 95%.Continuing the analysis, the optimizer took 38.25 minutes to execute 51 experiments before converging in the case of the reduced range problem which is equal to approximately 80% of the reduced experimental area. In the case of the original variation range, the optimizer took 2.1 hours to carry out 168 experiments before converging, which is equal to about 5.3% of the total experimental area.As it can be seen in Figure4 for both models tested with the proposed optimization procedure, the output variable responses were statistically identical, while the time needed to reach these results dropped considerably for reduced ranges.Figure 4
Results comparison with optimization of both simulation models.


## 4.1. Monoobjective Case

The modeled situation corresponds to a quality control cell in a fiber-optic transponder company. The cell is responsible for a series of tests which will ensure the approval or failure of the equipment produced at the plant. All verification and validation phases were appropriately undertaken, thus providing a consistent model.For this study object, the decision variables were defined as the number of operators responsible for carrying out quality control tests, denoted by(
x
1
,
x
2
,
x
3
), and the number of pieces of equipment for test types 1, 2, and 3, denoted by (
x
4
,
x
5
,
x
6
). All variables were defined as integers, with a lower bound equal to one and an upper one equal to five. Table 1 presents this information.Table 1
Decision variables, types, and limits for the study object.
Variable
Type
Lower bound
Upper bound
x
1
Number of operators, type 1
Integer
1
5
x
2
Number of operators, type 2
Integer
1
5
x
3
Number of operators, type 3
Integer
1
5
x
4
Number of pieces of equipment for test, type 1
Integer
1
5
x
5
Number of pieces of equipment for test, type 2
Integer
1
5
x
6
Number of pieces of equipment for test, type 3
Integer
1
5The optimization objective was to find the combination of decision variables which would maximize the total number of inspected products in the quality control center, denoted by (y
1). For the problem in question, considering the number of decision variables (
6
) and their maximum variation (1–5), there are a total of 15.625 = (56) possible scenarios for the search space.Considering the quantity of decision variables, the variation of levels of each variable, and the rule for the minimum number of DMUs proposed by Cooper et al. [16], the orthogonal array L25 was chosen. Seeing that there are six input variables and one output variable, there would be a minimum of 21 DMUs (runs) by following the classic rule, justifying the use of array L25. With a defined array, the experimental matrix was generated and presented in Table 2.Table 2
Experimental matrix and results.
DMU
Input variables
Output variables
x
1
x
2
x
3
x
4
x
5
x
6
y
1 minimum
y
1 average
y
1 maximum
1
1
1
1
1
1
1
3,922
3,977
4,032
2
1
2
2
2
2
2
3,936
3,979
4,022
3
1
3
3
3
3
3
3,919
3,978
4,037
4
1
4
4
4
4
4
3,932
3,985
4,038
5
1
5
5
5
5
5
3,939
3,982
4,025
6
2
1
2
3
4
5
3,939
3,982
4,025
7
2
2
3
4
5
1
3,925
3,978
4,031
8
2
3
4
5
1
2
3,930
3,982
4,034
9
2
4
5
1
2
3
3,938
3,979
4,020
10
2
5
1
2
3
4
3,920
3,983
4,046
11
3
1
3
5
2
4
3,918
3,982
4,046
12
3
2
4
1
3
5
3,910
3,979
4,048
13
3
3
5
2
4
1
3,924
3,983
4,042
14
3
4
1
3
5
2
3,926
3,981
4,036
15
3
5
2
4
1
3
3,937
3,983
4,029
16
4
1
4
2
5
3
3,938
3,980
4,022
17
4
2
5
3
1
4
3,930
3,978
4,026
18
4
3
1
4
2
5
3,932
3,982
4,032
19
4
4
2
5
3
1
3,922
3,981
4,040
20
4
5
3
1
4
2
3,921
3,982
4,043
21
5
1
5
4
3
2
3,937
3,980
4,023
22
5
2
1
5
4
3
3,932
3,982
4,032
23
5
3
2
1
5
4
3,930
3,982
4,034
24
5
4
3
2
1
5
3,933
3,982
4,031
25
5
5
4
3
2
1
3,933
3,982
4,031In the following, 25 scenarios were simulated from array L25 with 30 replications of a month’s time of operation in quality control. The simulations were carried out on a computer with an Intel processor (Core 2 Duo) 1.58 GHZ, a 2 GB RAM, and a 64-bit Microsoft operational system platform.Data for each output variable were stored for the superefficiency calculations. Nearly 26 minutes were spent to process the 25 scenarios, considering all replications. Maximum, minimum, and average values were stored for 30 replications for each scenario with the goal of using them in the triangular membership function. Results for output variabley
1 are shown in Table 2.For the calculation of superefficiency related to each DMU with a Fuzzy-DEA-BCC model, the softwareThe General Algebraic Modeling (GAMS) [33] was used, version 22.8.1, using the solver CPLEX, version 11.0, adapted for the specific calculation.With these results, the superefficiency value of the DMU can be related to each scenario for the pessimistic and optimistic scenarios by means of anα variation, summing up to 11 scenarios. These values are presented in Tables 3 and 4.Table 3
Superefficiency under uncertainty for a pessimistic scenario.
Scenario
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11
α
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
DMU 1
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
DMU 2
1.750
1.742
1.733
1.715
1.691
1.625
1.535
1.407
1.253
1.079
1.000
DMU 3
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 4
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.207
6.000
6.714
2.833
DMU 5
2.000
2.000
2.000
2.000
2.000
2.000
2.000
1.000
1.000
1.000
1.000
DMU 6
5.000
5.000
5.000
5.000
5.000
5.000
5.000
4.000
2.929
2.462
1.446
DMU 7
1.014
1.012
1.010
1.007
1.004
1.000
1.000
1.000
1.000
1.000
1.000
DMU 8
1.023
1.033
1.044
1.057
1.074
1.094
1.121
1.156
1.202
1.268
1.316
DMU 9
2.222
2.178
2.133
2.089
2.044
1.976
1.883
1.676
1.528
1.330
1.052
DMU 10
1.000
1.000
1.000
1.000
1.000
1.000
1.047
1.161
1.306
1.480
2.000
DMU 11
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.064
1.294
1.571
DMU 12
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 13
1.000
1.000
1.000
1.000
1.000
1.000
1.067
1.156
1.280
1.475
2.000
DMU 14
1.000
1.000
1.000
1.000
1.000
1.021
1.049
1.085
1.128
1.181
1.200
DMU 15
1.933
2.037
2.208
2.625
3.111
3.500
3.750
4.167
3.182
2.500
2.000
DMU 16
1.333
1.329
1.323
1.316
1.307
1.296
1.282
1.263
1.233
1.161
1.074
DMU 17
1.222
1.207
1.190
1.168
1.143
1.111
1.071
1.018
1.000
1.000
1.000
DMU 18
1.308
1.310
1.313
1.315
1.319
1.323
1.327
1.332
1.328
1.296
1.250
DMU 19
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.063
1.139
1.224
1.318
DMU 20
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.011
1.167
1.288
1.357
DMU 21
1.318
1.320
1.322
1.325
1.329
1.333
1.351
1.367
1.367
1.346
1.211
DMU 22
1.263
1.271
1.280
1.290
1.303
1.319
1.341
1.352
1.367
1.389
1.400
DMU 23
1.222
1.246
1.275
1.311
1.356
1.415
1.494
1.611
1.676
1.711
1.400
DMU 24
1.231
1.238
1.245
1.255
1.267
1.281
1.300
1.325
1.360
1.430
1.500
DMU 25
1.683
1.698
1.714
1.732
1.758
1.800
1.882
2.026
1.769
1.442
1.333Table 4
Superefficiency under uncertainty for an optimistic scenario.
Scenario
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11
α
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
DMU 1
2.441
2.424
2.403
2.377
2.345
2.304
2.249
2.171
2.054
2.000
2.000
DMU 2
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 3
1.081
1.040
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 4
1.111
1.186
1.268
1.359
1.459
1.571
1.697
1.839
2.000
2.333
2.833
DMU 5
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 6
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.137
1.446
DMU 7
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 8
1.143
1.176
1.210
1.243
1.284
1.314
1.351
1.399
1.397
1.356
1.316
DMU 9
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.052
DMU 10
3.000
3.067
3.200
3.600
3.125
3.182
3.261
3.378
4.000
3.091
2.000
DMU 11
1.875
1.897
1.924
1.958
2.000
3.500
3.875
3.833
3.000
2.575
1.571
DMU 12
2.667
30.500
10.500
3.833
2.207
2.000
1.765
1.595
1.382
1.164
1.000
DMU 13
1.667
1.714
1.769
1.833
1.909
2.000
2.222
2.500
2.857
2.727
2.000
DMU 14
1.000
1.000
1.000
1.009
1.029
1.053
1.080
1.111
1.149
1.194
1.200
DMU 15
1.000
1.000
1.000
1.000
1.000
1.035
1.128
1.221
1.331
1.476
2.000
DMU 16
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.074
DMU 17
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
DMU 18
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.022
1.137
1.250
DMU 19
1.500
1.500
1.500
1.500
1.500
1.500
1.474
1.447
1.419
1.389
1.318
DMU 20
1.389
1.399
1.411
1.425
1.441
1.460
1.482
1.529
1.611
1.493
1.357
DMU 21
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.211
DMU 22
1.000
1.000
1.000
1.000
1.000
1.000
1.079
1.174
1.278
1.354
1.400
DMU 23
1.000
1.000
1.000
1.000
1.000
1.000
1.008
1.062
1.140
1.257
1.400
DMU 24
1.000
1.000
1.000
1.008
1.068
1.145
1.246
1.359
1.407
1.454
1.500
DMU 25
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.069
1.159
1.247
1.333For the calculation of superefficiency for each scenario, the average geometric value was extracted for both pessimistic and optimistic scenarios. The results found are presented in Table5. With these data, it was possible to calculate the average of each DMU and then rank them in function of their superefficiency values (see Table 5).Table 5
Global matrix for geometric averages for pessimistic and optimistic scenarios.
DMU
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11
Average superefficiency
Ranking
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1
2.210
2.202
2.192
2.180
2.166
2.147
2.121
2.084
2.027
2.000
2.000
2.121
1
2
1.323
1.320
1.317
1.310
1.300
1.275
1.239
1.186
1.119
1.039
1.000
1.221
16
3
1.040
1.020
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.005
24
4
1.054
1.089
1.126
1.166
1.208
1.254
1.303
1.489
3.464
3.958
2.833
1.813
5
5
1.414
1.414
1.414
1.414
1.414
1.414
1.414
1.000
1.000
1.000
1.000
1.264
11
6
2.236
2.236
2.236
2.236
2.236
2.236
2.236
2.000
1.711
1.673
1.446
2.044
2
7
1.007
1.006
1.005
1.004
1.002
1.000
1.000
1.000
1.000
1.000
1.000
1.002
25
8
1.081
1.102
1.124
1.147
1.174
1.199
1.231
1.272
1.296
1.311
1.316
1.205
18
9
1.491
1.476
1.461
1.445
1.430
1.406
1.372
1.295
1.236
1.153
1.052
1.347
10
10
1.732
1.751
1.789
1.897
1.768
1.784
1.847
1.980
2.286
2.139
2.000
1.907
4
11
1.369
1.377
1.387
1.399
1.414
1.871
1.969
1.958
1.786
1.825
1.571
1.630
7
12
1.633
5.523
3.240
1.958
1.486
1.414
1.328
1.263
1.176
1.079
1.000
1.918
3
13
1.291
1.309
1.330
1.354
1.382
1.414
1.540
1.700
1.912
2.006
2.000
1.567
8
14
1.000
1.000
1.000
1.005
1.015
1.037
1.064
1.098
1.138
1.188
1.200
1.068
22
15
1.390
1.427
1.486
1.620
1.764
1.903
2.057
2.255
2.058
1.921
2.000
1.807
6
16
1.155
1.153
1.150
1.147
1.143
1.139
1.132
1.124
1.110
1.078
1.074
1.128
21
17
1.106
1.099
1.091
1.081
1.069
1.054
1.035
1.009
1.000
1.000
1.000
1.049
23
18
1.144
1.145
1.146
1.147
1.148
1.150
1.152
1.154
1.165
1.214
1.250
1.165
19
19
1.225
1.225
1.225
1.225
1.225
1.225
1.214
1.240
1.271
1.304
1.318
1.245
14
20
1.179
1.183
1.188
1.194
1.200
1.208
1.218
1.244
1.371
1.387
1.357
1.248
13
21
1.148
1.149
1.150
1.151
1.153
1.155
1.162
1.169
1.169
1.160
1.211
1.162
20
22
1.124
1.127
1.131
1.136
1.142
1.149
1.203
1.260
1.321
1.372
1.400
1.215
17
23
1.106
1.116
1.129
1.145
1.164
1.189
1.227
1.308
1.382
1.466
1.400
1.239
15
24
1.109
1.112
1.116
1.125
1.163
1.211
1.273
1.342
1.383
1.442
1.500
1.252
12
25
1.297
1.303
1.309
1.316
1.326
1.342
1.372
1.472
1.432
1.341
1.333
1.349
9Through a superefficiency analysis, it was possible to rank the DMUs in order of efficiency. For the problem in question, DMU 1 is the most efficient, followed by DMU 6. Both are highlighted in Table5.With the identification of the two most efficient DMUs and based on the experimental matrix in Table2, a new interval can be identified for each decision variable in which a better set of solutions is expected. The new intervals for each decision variable are presented in Table 6. Variable x
2 stands out, as its value was already defined, being equal to x
2
=
1, reducing the number of decision variables to 5.Table 6
Decision variables, types, and new limits for the study object.
Variable
Type
Lower bound
Upper bound
x
1
Number of operators, type 1
Integer
1
2
x
3
Number of operators, type 3
Integer
1
2
x
4
Number of pieces of equipment for test, type 1
Integer
1
3
x
5
Number of pieces of equipment for test, type 2
Integer
1
4
x
6
Number of pieces of equipment for test, type 3
Integer
1
5With the reduction of the variation interval for each decision variable, the search space for the best solution was reduced from 15,625 to 240—a reduction of 98.4%.To confirm the efficiency of search space reduction, the optimizer SimRummer [34] was utilized for the simulation model optimization. SimRunner stands as popular optimization software [35] in operation and being sold in conjunction with the ProModel simulator. According to Kim et al. [36] and Banks et al. [4], this optimization software finds an optimal solution using as its search method a metaheuristic called Genetic Algorithm that as pointed by Ólafsson [37] mimics the process of natural selection.The optimizer was set to the same conditions and objectives; however, two optimizations were conducted for the same problem. One used a reduction in the variation interval (Table6), and the other one used the original problem (Table 1) conditions. The results found can be seen in Table 7.Table 7
Optimization results.
Decision variables
Solution
New range
Original range
x
1
2
5
x
2
1
3
x
3
1
2
x
4
2
2
x
5
1
1
x
6
3
5
Output variable
Responses
New range
Original range
y
1
3,986
3,987
Confidence interval fory
1 (95%)
(3,974–3,998)
(3,977–3,997)The responses presented by the optimizer for the decision variables were only equal for variablesx
4 and x
5. For the other variables, the optimizer arrived at different values. As for the solution of y
1, with search space reduction, the optimizer reached a value statistically equal to both optimized cases.The optimizer for the optimization problem with its range being reduced carried out 88 experiments before converging, which is equal to 36.67% of the experimental area (240 scenarios), taking a little more than 2.25 hours. For the problem with original range settings, the optimizer carried out 183 experiments, equal to slightly more than 1.17% of the total experimental area (15,625 scenarios), taking 4.8 hours.


## 4.2. Multiobjective Cases

The second study object represents a production cell in a Brazilian telecommunications company which produces fiber-optical equipment. This model, which is similar to the previous case, had already passed through the phases of verification and validation prior to being used for this paper in the simulation optimization technique. For the presented study object, the decision variables were defined as the number of pieces of inspection equipment, types 1 and 2, denoted by(
x
1
,
x
2
), and the number of employees, types 1, 2, and 3, denoted by (
x
3
,
x
4
,
x
5
), who carry out the activities presented in the model. The variables were defined as integers, with a lower bound of 1 and an upper bound of 5. Table 8 presents this information.Table 8
Decision variables, types, and limits for the second study object.
Variable
Type
Lower bound
Upper bound
x
1
Number of pieces of inspection equipment, type 1
Integer
1
5
x
2
Number of pieces of inspection equipment, type 2
Integer
1
5
x
3
Number of employees, type 1
Integer
1
5
x
4
Number of employees, type 2
Integer
1
5
x
5
Number of employees, type 3
Integer
1
5The optimization objective was to find the combination of decision variables which would maximize throughput (y
1) and production cell profit (y
2). It can be seen for the problem in question that, considering the number of variables (
5
) and their maximum variation (1–5), there are a total of 3,125 = (55) possible scenarios for the search space for the best configuration.In order to meet the fundamental rule for the number of DMUs, an L25 orthogonal array was used. An experimental matrix was generated which is presented in Table9. The scenarios for the experimental matrix were simulated in ProModel.Table 9
Experimental matrix and results.
DMU
Input variables
Output variables
x
1
x
2
x
3
x
4
x
5
y
1 minimum
y
1 average
y
1 maximum
y
2 minimum
y
2 average
y
2 maximum
1
1
1
1
1
1
643
863
1,083
381,096
682,082
983,067
2
1
2
2
2
2
1,295
1,395
1,495
826,452
1,097,340
1,368,228
3
1
3
3
3
3
1,188
1,381
1,574
852,107
1,080,373
1,308,639
4
1
4
4
4
4
1,235
1,399
1,563
726,952
1,091,904
1,456,857
5
1
5
5
5
5
1,148
1,403
1,658
690,053
1,088,349
1,486,646
6
2
1
2
3
4
1,064
1,286
1,508
748,800
1,007,407
1,266,015
7
2
2
3
4
5
1,219
1,383
1,547
842,802
1,067,833
1,292,864
8
2
3
4
5
1
1,175
1,391
1,607
783,439
1,081,793
1,380,148
9
2
4
5
1
2
573
868
1,163
465,318
681,706
898,094
10
2
5
1
2
3
895
1,140
1,385
515,279
896,302
1,277,326
11
3
1
3
5
2
1,153
1,404
1,655
862,630
1,095,382
1,328,135
12
3
2
4
1
3
736
847
958
425,035
660,137
895,239
13
3
3
5
2
4
1,129
1,383
1,637
804,861
1,077,000
1,349,140
14
3
4
1
3
5
894
1,167
1,440
532,523
911,184
1,289,846
15
3
5
2
4
1
1,275
1,397
1,519
840,255
1,096,732
1,353,210
16
4
1
4
2
5
1,178
1,396
1,614
704,903
1,084,938
1,464,973
17
4
2
5
3
1
1,178
1,409
1,640
710,034
1,101,006
1,491,979
18
4
3
1
4
2
1,013
1,123
1,233
516,328
871,771
1,227,214
19
4
4
2
5
3
1,215
1,385
1,555
727,103
1,071,027
1,414,953
20
4
5
3
1
4
586
828
1,070
426,589
638,185
849,782
21
5
1
5
4
3
1,264
1,406
1,548
720,561
1,093,357
1,466,152
22
5
2
1
5
4
974
1,157
1,340
582,466
888,310
1,194,154
23
5
3
2
1
5
561
829
1,097
406,997
635,928
864,860
24
5
4
3
2
1
1,109
1,301
1,493
738,318
1,014,712
1,291,106
25
5
5
4
3
2
1,214
1,379
1,544
834,734
1,068,090
1,301,447Thirty replications were simulated for each scenario, which corresponded to a month of operation in the production cell, and the data for each output variable were stored for the superefficiency calculation. Simulation of the 25 scenarios and their 30 replications took slightly more than 16.25 minutes. The results (minimum, maximum, and mean values) for the outputsy
1 and y
2 for each DMU are shown in Table 9.For the calculation of superefficiency related to each DMU for the Fuzzy-DEA-BCC model, the same procedure was employed. With the output variable results, each scenario can be related to each DMU superefficiency value for pessimistic and optimistic scenarios. In the following, the scenarios are lumped into a geometric average. The final result is presented in Table10. Thus, the average can be calculated for each DMU, making it possible to rank them in terms of superefficiency.Table 10
Global matrix for geometric average for pessimistic and optimistic scenarios.
DMU
S1
S2
S3
S4
S5
S6
S7
S8
S9
S10
S11
Average superefficiency
Ranking
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
2.000
4
2
3.803
3.661
3.512
3.360
3.203
3.043
3.006
2.915
2.758
2.700
2.792
3.159
1
3
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
18
4
1.019
1.022
1.024
1.028
1.036
1.055
1.168
1.085
1.051
1.021
1.000
1.046
16
5
1.008
1.010
1.012
1.016
1.067
1.285
1.432
1.547
1.591
1.674
2.200
1.349
8
6
1.281
1.278
1.269
1.258
1.248
1.238
1.228
1.219
1.209
1.200
1.191
1.238
11
7
0.805
0.808
0.811
0.811
0.810
0.811
0.814
0.815
0.813
0.812
0.812
0.811
24
8
1.270
1.271
1.271
1.271
1.271
1.271
1.271
1.271
1.271
1.271
1.271
1.271
10
9
1.016
1.015
1.014
1.013
1.012
1.010
1.009
1.009
1.009
1.009
1.009
1.011
17
10
1.030
1.046
1.066
1.087
1.108
1.130
1.153
1.176
1.189
1.188
1.188
1.124
15
11
1.952
2.000
2.069
2.178
2.353
2.691
7.030
1.683
1.665
1.647
1.431
2.427
3
12
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
19
13
0.977
0.982
0.982
0.983
0.984
0.984
0.985
0.986
0.987
0.987
0.988
0.984
22
14
1.652
1.614
1.574
1.535
1.495
1.455
1.414
1.316
1.213
1.126
1.086
1.407
7
15
1.539
1.555
1.566
1.582
1.614
1.649
2.709
1.953
1.854
1.655
1.578
1.750
5
16
1.448
1.461
1.481
1.500
1.496
1.493
1.494
1.493
1.491
1.491
1.490
1.485
6
17
1.694
1.704
1.719
1.739
3.999
2.856
2.873
3.015
3.138
2.808
2.009
2.505
2
18
1.106
1.114
1.123
1.132
1.142
1.152
1.164
1.169
1.173
1.176
1.176
1.148
14
19
0.847
0.858
0.868
0.877
0.887
0.900
0.902
0.919
0.961
0.971
0.975
0.906
23
20
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
20
21
1.313
1.316
1.318
1.318
1.315
1.316
1.318
1.322
1.337
1.352
1.400
1.330
9
22
1.123
1.126
1.138
1.150
1.164
1.179
1.195
1.213
1.232
1.240
1.240
1.182
13
23
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
21
24
1.235
1.231
1.227
1.224
1.220
1.216
1.213
1.210
1.207
1.205
1.205
1.218
12
25
0.750
0.753
0.756
0.755
0.757
0.757
0.757
0.758
0.758
0.758
0.758
0.756
25Based on the superefficiency value rankings (Table3) in this study, DMU 2 proved to be the most efficient followed by DMU 17. Both are highlighted in Table 10.Once these two DMUs are identified as the most efficient, a new interval for each decision variable can be redefined. In doing so, the optimization search space is reduced. Table9 shows that variable x
2 presented the same value for the most efficient DMUs, thus causing its value to be defined as x
2
=
2, reducing the number of variables from five to four. The new intervals for the other decision variables are presented in Table 11.Table 11
Decision variables, types, and new limits for the second study object.
Variable
Type
Lower bound
Upper bound
x
1
Number of pieces of inspection equipment, type 1
Integer
1
4
x
3
Number of employees, type 1
Integer
2
5
x
4
Number of employees, type 2
Integer
2
3
x
5
Number of employees, type 3
Integer
2
1By reducing the variation interval for each decision variable, the search space for the optimal solution was reduced from 3,125 to 64—a reduction of nearly 98%.To test the efficiency and robustness of the new search space, SimRunner was set to carry out the simulation model optimization, aiming to maximize total cell production(
y
1
) and total profit (
y
2
), for both the original variation range (Table 8) and the reduced variation range (Table 11). Results are seen in Table 12.Table 12
Optimization results.
Decision variables
Solution
New range
Original range
x
1
1
1
x
2
2
5
x
3
2
4
x
4
3
2
x
5
1
2
Output variable
Responses
New range
Original range
y
1
1,402
1,403
Confidence Interval fory
1 (95%)
(1,395–1,409)
(1,397–1,410)
y
2
1,103,794
1,103,595
Confidence Interval fory
2 (95%)
(1,096,972–1,110,617)
(1,095,465–1,111,727)The responses presented by the optimizer for the decision variables were only equal for variablex
1
=
1. For the other decision variables, the optimizer arrived at different values. In the case of the original variation range, the responses indicated the need to hire more employees and purchase more equipment, except for variable x
4. As for the solutions of y
1 and y
2, the results were statistically equal for a confidence level of 95%.Continuing the analysis, the optimizer took 38.25 minutes to execute 51 experiments before converging in the case of the reduced range problem which is equal to approximately 80% of the reduced experimental area. In the case of the original variation range, the optimizer took 2.1 hours to carry out 168 experiments before converging, which is equal to about 5.3% of the total experimental area.As it can be seen in Figure4 for both models tested with the proposed optimization procedure, the output variable responses were statistically identical, while the time needed to reach these results dropped considerably for reduced ranges.Figure 4
Results comparison with optimization of both simulation models.


## 5. Conclusions and Recommendations for Future Research

This paper has proposed a method for reducing search space for simulation optimization problems which has provided search space reductions of roughly 98% and significant reductions in convergence time without compromising the response quality.These results were reached through the proposal of using the Fuzzy Theory along with a DEA-BCC model which enabled scenario analysis in the face of uncertainty, which is a common reality for discrete-event simulation models, considering that they most commonly deal with stochastic, dynamic, and interrelated environments.With an output variable analysis taking uncertainty into account, the use of Fuzzy-DEA-BCC allowed an efficiency analysis of pessimistic and optimistic scenarios. This, in turn, enabled the ranking of the most efficient DMUs.Upon determining the two most efficient DMUs, a new range for each decision variable was established, permitting the optimization software to concentrate on the region of the greatest efficiency, according to the analysis conducted with Fuzzy-DEA-BCC.As a means of validating this paper’s proposal, a widely used commercial optimizer was used to check whether the method was indeed able to limit the search region to the area containing the best solutions or not. In order to do so, the optimizer ran the simulation model under both variations: original and reduced with Fuzzy-DEA-BCC. For both study objects, the method was able to provide a response of equal quality from the reduced search space when both variation ranges were compared. There was a significant reduction in convergence time with reduced search space.Finally, it is worth mentioning that the Taguchi arrays proved to be practical and reliable, as they represented the search space, exploring the maximum possible diversity of levels present in each decision variable, which would have been difficult using classical experimental techniques that use two or three levels.The possibilities for future research include(i)
utilizing GPDEA-BCC [38] which improves DEA discrimination even without a number of DMUs to meet Cooper, Seiford, and Tone’s rule [16];
(ii)
conducting tests with continuous decision variables;
(iii)
applying other arrays or strategic experiments which substitute the necessity for orthogonal arrays;
(iv)
testing the method proposed in this paper, using other optimizers, such as multiple comparison procedure, optimal computing budget allocation, and nested partitions that seek to increase the efficiency of the optimization process.



---

*Source: 450367-2014-05-19.xml*