Data envelopment analysis (DEA) is a mathematical programming methodology for evaluating and measuring the relative efficiencies of a set of decision making units (DMUs) that use multiple inputs to produce multiple outputs. Due to its solid underlying mathematical basis and wide applications to real-world problems, much effort has been devoted to the DEA methods since the pioneering work of [1, 2] summarized the major research in DEA over the last 30 years, which provided a good research framework.

In the conventional DEA models, all the data are assumed to have the form of specific numerical values which are "known exactly". However, this situation may not always be true. In reality, the data of real-world problems more often than not are uncertain—not known exactly at the time the problem is being solved. In applications one cannot ignore the possibility that even a small uncertainty in the data can make the nominal optimal solution to the problem completely meaningless from a practical viewpoint. [3] showed that a small perturbation on data of linear programming problem could lead heavily infeasible solution. As a linear programming based approach, DEA will never be able to escape from the impact of uncertainty, i.e., a small perturbation on output data of DMU could make a big change on the efficiencies, so the results of the ranking could be unreliable. Consequently, there exists a real need of a methodology capable of generating a robust solution, one that is immunized against the effect of data uncertainty.

In this paper, we consider the perturbation in output data and propose a new robust DEA model based on the adaptation of recently developed robust optimization approaches proposed by [4]; [5] proposed a robust DEA model with consideration of uncertainty on output parameters for the performance assessment of electricity distribution companies. However, the assumption of symmetric data uncertainty made in their paper could be too restrictive for many real-world applications; [6] considered the DEA with uncertain data, they proposed a second order cone model. In this paper, we relax the assumption of symmetric data uncertainty and construct the robust formulation for asymmetric data uncertainty, the model we present is a linear programming problem.

This paper is organized as follows. We first present the fundamentals of robust optimization in Section 2. In Section 3, we illustrate the classical DEA model and propose the robust DEA model when the output data are uncertain. In Section 4, we demonstrate some experimental results. Finally, in Section 5, we sum up our conclusions.

To present the robust optimization method, consider a linear programming problem

$ \begin{eqnarray} &\max &\quad c^{T}x, \nonumber\\ &\text{s.t.} &Ax \leq b, \\ && l\leq x\leq u, \nonumber \end{eqnarray} $

(2.1)

where $c, l, u\in \mathbb{R}^{n}$, $b\in \mathbb{R}^{m}$, $A=(a_{ij})$ is a $m\times n$ matrix, and $x\in \mathbb{R}^{n}$ is the vector of decision variables. We assume, without loss of generality, that only the elements of the matrix $A$ are subject to uncertainty. In fact, if $c$ and $b$ are also uncertain, we can rewrite the problem as

$ \begin{eqnarray*} &\max &\quad \tilde{c}^{T}\tilde{x}, \nonumber\\ &\text{s.t.} &\tilde{A} \tilde{x} \leq 0, \\ && \tilde{l}\leq \tilde{x}\leq \tilde{u} \nonumber \end{eqnarray*} $

with $\tilde{x}=(z, x, y)$, $\tilde{c}=(1, 0, 0)$, $\tilde{A}= \left( \begin{array}{ccc} 0&A&-b \\ 1&-c^{T}&0 \end{array} \right) $, $\tilde{l}=(z_{L}, l, 1)$, and $\tilde{u}=(z_{U}, u, 1)$, where $z_{L}$ and $z_{U}$ are finite constants as long as all the components of $c$ are bounded. In this new formulation, only the matrix $\tilde{A}$ contains uncertain data.

Consider a particular row $i$ of the matrix $A$ and let $J_{i}$ be the set of coefficients in row $i$ that are subject to uncertainty. Each entry $a_{ij}, j\in J_{i}$ is modeled as a symmetric and bounded random variable that takes values in $[\bar{a}_{ij}-\hat{a}_{ij}, \bar{a}_{ij}+\hat{a}_{ij}]$.

Given the data uncertainty structure for $A$, the traditional linear optimization methodology elect to solve the nominal formulation, where each random $a_{ij}$ is replaced by its mean value:

$ \begin{eqnarray} &\max &\quad c^{T}x, \nonumber\\ &\text{s.t.} &\sum\limits_{j}\bar{a}_{ij}x_{j}\leq b_{i}, \quad \forall i, \\ && l\leq x\leq u. \nonumber \end{eqnarray} $

(2.2)

Small data uncertainty is just ignored as if the given ("nominal") data were exact, and the resulting nominal solution is what recommended for use, in hope that small data uncertainties will not affect significantly the feasibility and optimality properties of this solution. However, this hope is not necessarily justified. The analysis of linear optimization problems from the NETLIB collection reported in Ben-Tal & Nemirovski [3] revealed that for 13 of 94 NETLIB problems, random 0.01percent perturbations of the uncertain data can make the nominal optimal solution severely infeasible: with a non-negligible probability, it violates some of the constraints by 50 percent and more. Consequently, this leads one to consider a solution that is guaranteed to satisfy $Ax\leq b$ for all realizations of the random $a_{ij}$'s, while maximizing the objective value. To obtain such a solution, Soyster [7] proposed the following robust formulation :

$ \begin{eqnarray} &\max &\quad c^{T}x, \nonumber\\[2pt] &\text{s.t.} &\sum\limits_{j}\bar{a}_{ij}x_{j}+\sum\limits_{j\in J_{i}}\hat{a}_{ij}y_{j} \leq b_{i}, \quad \forall i, \nonumber\\[2pt] &&-y_{j}\leq x_{j}\leq y_{j}, \quad \forall j, \\[2pt] && l\leq x\leq u, \nonumber\\[2pt] && y\geq 0. \nonumber \end{eqnarray} $

(2.3)

It can be proved that the optimal solution of (2.3) is a feasible solution of (2.1) for every possible realization of $A$. Although the Soyster's method provides the most robust solution, it is also the most conservative in practice in the sense that the robust solution has an objective function value much worse than the objective function value of the solution of the nominal linear optimization problem. To address this conservatism, [8] proposed the following robust problem:

$ \begin{eqnarray} &\max &\quad c^{T}x, \nonumber\\[2pt] &\text{s.t.} &\sum\limits_{j}\bar{a}_{ij}x_{j}+\sum\limits_{j\in J_{i}}\hat{a}_{ij}y_{j}+\Omega_{i}\sqrt{\sum\limits_{j\in J_{i}}\hat{a}_{ij}^{2}z_{ij}^{2}} \leq b_{i}, \quad \forall i, \nonumber\\[2pt] &&-y_{ij}\leq x_{j}-z_{ij}\leq y_{ij}, \quad \forall i, j\in J_{i}, \\[2pt] && l\leq x\leq u, \nonumber\\[2pt] && y\geq 0, \nonumber \end{eqnarray} $

(2.4)

where $\Omega_{i}$, which is dependent upon the user's risk preference, is a user defined parameter and adjusts the trade-off between robustness and optimality. The authors have shown that the probability that the $i$ constraint is violated at most $\text{exp}(-\Omega_{i}^2)/2$. This robust model is less conservative than Model (2.3) as every feasible solution of the latter problem is a feasible solution to the former problem. However, Model (2.4) is a second order cone problem, which is more demanding computationally than the linear Model (2.3).

To overcome this difficulty, [4] presented a new robust formulation. One attractive aspect of this method is that the new formulation is a linear programming problem. The authors introduced a budget parameter $\Gamma_{i}\in [0, |J_{i}|]$ for each row $i=1, \ldots, m$ of the matrix $A$, which is a user defined parameter that adjusts the robustness of the model and interpreted as the maximum number of uncertain parameters allowed to take their worst case value. The role of $\Gamma_{i}$ is in the following way:

If $\Gamma_{i}=0$, each $a_{ij}, j\in J_{i}$, is forced to take its mean value $\bar{a}_{ij}$. If $\Gamma_{i}=|J_{i}|$, each $a_{ij}, j\in J_{i}$, can take values from its range $[\bar{a}_{ij}-\hat{a}_{ij}, \bar{a}_{ij}+\hat{a}_{ij}]$. If $0< \Gamma_{i}<|J_{i}|$, $\lfloor\Gamma_{i}\rfloor$ elements among $a_{ij}, \forall j\in J_{i}$, can take values from their respective ranges; furthermore, if $\Gamma_{i}$ is not an integer, one other random element, says $a_{it_{i}}$, can take values from its reduced range $[\bar{a}_{it_{i}}-(\Gamma_{i}-\lfloor\Gamma_{i}\rfloor)\hat{a}_{it_{i}}, \bar{a}_{it_{i}}+ (\Gamma_{i}-\lfloor\Gamma_{i}\rfloor)\hat{a}_{it_{i}}]$; the remaining $|J_{i}|-\lceil\Gamma_{i} \rceil$ random elements are forced to take their respectively mean values. The role of the parameter $\Gamma_{i}$ is to adjust the robustness against the level of conservatism of the solution. As surmised in [4], nature could be restricted in its behavior in that only a subset of the random elements actually deviate from their respective mean values, in order to adversely affect the solution.

Given $\Gamma_{i}$ for all $i$, [4] sought a solution that maximizes the objective value under the restriction that it must remain feasible to (2.1) as long as up to $\Gamma_{i}$ elements out of $a_{ij}, \forall j\in J_{i}$, are allowed to change. Let $S_{i}$ be a subset of $ J_{i}$, such that $|S_{i}|=\lfloor\Gamma_{i}\rfloor$, and let $t_{i}\in J_{i}\backslash S_{i}$. To obtain such a solution, [4] constructed the following robust formulation:

$ \begin{eqnarray} &\max &\quad c^{T}x, \nonumber\\[2pt] &\text{s.t.} &\sum\limits_{j}\bar{a}_{ij}x_{j}+\max\limits_{S_{i}\cup{t_{i}}\subseteq J_{i}} \left\{ \sum\limits_{j\in S_{i}}\hat{a}_{ij}y_{j}+ (\Gamma_{i}-\lfloor\Gamma_{i}\rfloor)\hat{a}_{it_{i}}y_{t_{i}} \right\} \leq b_{i}, \quad \forall i, \nonumber\\[2pt] &&-y_{j}\leq x_{j}\leq y_{j}, \quad \forall j, \\[2pt] && l\leq x\leq u, \nonumber\\[2pt] && y\geq 0. \nonumber \end{eqnarray} $

(2.5)

[4] proved that the probability that the $i$ constraint violation is bounded above by $\text{exp}(-\Gamma_{i}^2)/2|J_{i}|$. Thus, as $\Gamma_{i}$ increases, more protection is given and the solution is more robust. The authors showed that the above nonlinear robust formulation can be recast as an equivalent linear programming formulation:

$ \begin{eqnarray} &\max &\quad c^{T}x, \nonumber\\[2pt] &\text{s.t.} &\sum\limits_{j}\bar{a}_{ij}x_{j}+\Gamma_{i}z_{i}+ \sum\limits_{j\in J_{i}}p_{ij} \leq b_{i}, \quad \forall i, \nonumber\\[2pt] && z_{i}+p_{ij}\geq \hat{a}_{ij}y_{j}, \quad \forall i, \, \forall j\in J_{i}, \nonumber\\[2pt] && p_{ij}\geq 0, \quad \forall i, \, \forall j\in J_{i}, \nonumber\\[2pt] && z_{i}\geq 0, \quad \forall i, \nonumber\\[2pt] &&-y_{j}\leq x_{j}\leq y_{j}, \quad \forall j, \\[2pt] && l\leq x\leq u, \nonumber\\[2pt] && y\geq 0. \nonumber \end{eqnarray} $

(2.6)

Assume that we deal with a set of $n$ DMUs converting $m$ inputs into $s$ outputs, with input-output vectors $(x_{j}, y_{j})$; $j=1, \cdots, n$, in which $x_{j}=(x_{1j}, \cdots, x_{mj})^T$ and $y_{j}=(y_{1j}, \cdots, y_{sj})^T$. Define $X=[x_{1}, x_{2}, \ldots, x_{n}]$ and $Y=[y_{1}, y_{2}, \cdots, y_{n}]$ as $m\times n$ and $s\times n$ matrices of inputs and outputs, respectively. The original fractional DEA model which is an input-oriented CCR model is presented as :

$ \begin{eqnarray} &\max &\quad \frac{v^T y_{j_{0}}}{u^T x_{j_{0}}}, \nonumber\\[2pt] &\text{s.t.} &\frac{v^Ty_{j}}{u^Tx_{j}} \leq 1, \, j=1, \cdots, n, \\[2pt] && v, u \geq 0, \nonumber \end{eqnarray} $

(3.1)

which evaluates the relative efficiencies of $n$ DMUs by maximizing the ratio of weighted summation of outputs to weighted summation of inputs. $u$ and $v$ are weight vectors associated with inputs and outputs, respectively. In addition, $x_{j_{0}}$ and $y_{j_{0}}$ are the input and output for the DMU under evaluation. This model is a nonlinear programming problem, and it is equivalent to the following linear programming problem which is more computational convenient:

$ \begin{eqnarray} &\max &\quad \mu^T y_{j_{0}}, \nonumber\\[2pt] &\text{s.t.} &\mu^T y_{j}-\nu^T x_{j} \leq 0, \, j=1, \ldots, n, \\[2pt] && \nu^T x_{j_{0}}=1, \nonumber\\[2pt] && \mu, \nu\geq 0. \nonumber \end{eqnarray} $

(3.2)

We assume that only the output data cannot be exactly obtained due to the existence of uncertainty, so as to avoid the appearance of uncertainty in the input related equality constraint. They are only known to lie within the upper and lower bounds represented by the range $[y_{rj}^L, y_{rj}^U]$, where $y_{rj}^L>0$. Assume the mean of the random variable $y_{rj}$ is $\bar{y}_{rj}$. In order to avoid the appearance of uncertainty in objective function, we express the objective function as $\max z$, and add the constraint $z-\mu^Ty_{j_{0}}\leq0$ into the constraints. In order to generate an uncertainty-immune solution, we adopt the robust optimization technique proposed by [4]. Let $J_{j}$ be the set of coefficients in column $j$ of matrix $Y$ that are subject to uncertainty. Choose $\Gamma_{j}\in [0, |J_{j}|]$. Let

$ \begin{eqnarray*}\Re(\Gamma_{j})&=&\{y_j=(y_{1j}, \cdots, y_{sj})|y_{ij}\in [\bar{y}_{ij}-\beta_{ij}(\bar{y}_{ij}-y_{ij}^L), \bar{y}_{ij}+\beta_{ij}(y_{ij}^U-\bar{y}_{ij})], \\ && \forall i;0\leq \beta_{ij}\leq 1, \forall i; \sum\limits_{i\in J_{j}}\beta_{ij}\leq \Gamma_{j}, \text{ at most one }\beta_{ij} \text{is fractional}\}.\end{eqnarray*} $

The robust DEA model with output uncertainty is expressed as following:

$ \begin{eqnarray} &\max &\quad z, \nonumber\\[2pt] &\text{s.t.}&\min\limits_{y_{j_{0}}\in \Re(\Gamma_{j_{0}}) }\mu^Ty_{j_{0}}\geq z, \\[2pt] &&\max\limits_{y_{j}\in \Re(\Gamma_{j})}\mu^Ty_{j}-\nu^Tx_{j} \leq 0, \, j=1, \ldots, n, \nonumber\\[2pt] && \nu^Tx_{j_{0}}=1, \nonumber\\[2pt] && \mu, \nu\geq 0. \nonumber \end{eqnarray} $

(3.3)

Let $S_{j}$ be a subset of $J_{j}$ such that $|S_{j}|=\lfloor\Gamma_{j}\rfloor$, and let $t\in J_{j}\backslash S_{j}$. It can be deduced that

$\begin{eqnarray*} \min\limits_{y_{j_{0}}\in \Re(\Gamma_{j_{0}}) }\mu^Ty_{j_{0}}=\mu^T\bar{y}_{j_{0}}-\max\limits_{S_{j_{0}}\cup \{t\}\subseteq J_{j_{0}}}\{ \sum\limits_{i\in S_{j_{0}}}\mu_{i} y_{ij_{0}}^L+(\Gamma_{j_{0}}-\lfloor\Gamma_{j_{0}}\rfloor)\mu_{t}y_{tj_{0}}^L \}, \end{eqnarray*} $

and

$\begin{eqnarray*} \max\limits_{y_{j}\in \Re(\Gamma_{j}) }\mu^Ty_{j}=\mu^T\bar{y}_{j}+\max\limits_{S_{j}\cup \{t\}\subseteq J_j}\{ \sum\limits_{i\in S_{j}}\mu_{i} y_{ij}^U+(\Gamma_{j}-\lfloor\Gamma_{j}\rfloor)\mu_{t}y_{tj}^U \}.\end{eqnarray*} $

Let

$\begin{eqnarray*} \psi_{j}(y, \Gamma_{j})=\max\limits_{S_{j}\cup \{t\}\subseteq J_j}\{ \sum\limits_{i\in S_{j}}\mu_{i} y_{ij}^U+(\Gamma_{j}-\lfloor\Gamma_{j}\rfloor)\mu_{t}y_{tj}^U \}, \end{eqnarray*} $

the function $\psi_{j}(y, \Gamma_{j})$ equals to the objective function of the following linear programming problem:

$ \begin{eqnarray} &\max &\quad \sum\limits_{i\in J_{j}}\mu_{i} y_{ij}^U \omega_{ij}, \nonumber\\[2pt] &\text{s.t.}&\sum\limits_{i\in J_{j}}\omega_{ij} \leq \Gamma_{j}, \\[2pt] && 0\leq \omega_{ij}\leq 1, \, \forall i\in J_{j}. \nonumber \end{eqnarray} $

(3.4)

Its dual problem is

$ \begin{eqnarray} &\min &\quad \Gamma_{j} z_{j}+ \sum\limits_{i\in J_{j}}p_{ij}, \nonumber\\[2pt] &\text{s.t.}&z_{j}+p_{ij}\geq \mu_{i} y_{ij}^U, \forall i\in J_{j}, \\[2pt] && p_{ij}\geq 0, \, \forall i\in J_{j}, \nonumber \\[2pt] && z_{j}\geq 0. \nonumber \end{eqnarray} $

(3.5)

By strong duality of linear programming, since problem (3.4) is feasible and bounded for all $\Gamma_{j}\in [0, |J_{j}|]$, then the dual problem (3.5) is also feasible and bounded and their objective values coincide. So we have that $\psi_{j}(y, \Gamma_{j})$ is equal to the objective function value of problem (3.5). Substituting to problem (3.3) we obtain that the robust DEA model is equivalent to the following linear optimization problem:

$ \begin{eqnarray} &\max &\quad z, \nonumber\\[2pt] &\text{s.t.}&\mu^T \bar{y}_{j_{0}}-(\Gamma_{j_{0}}z_{j_{0}}+\sum\limits_{i\in J_{j_{0}}}p_{ij_{0}})\geq z, \nonumber\\[2pt] && \mu^T \bar{y}_{j}+\Gamma_{j}z_{j}+\sum\limits_{i\in J_{j}}p_{ij}-\nu^T x_{j}\leq 0, \\[2pt] && \nu^T x_{j_{0}}=1, \nonumber\\[2pt] && z_{j}+p_{ij}\geq \mu_{i}y_{ij}^U, \quad \forall i\in J_{j}, \, \forall j, \nonumber\\[2pt] && z_{j_{0}}+p_{ij_{0}}\geq \mu_{i}y_{ij_{0}}^L, \quad \forall i\in J_{j_{0}}, \nonumber\\[2pt] && p_{ij}\geq 0, \quad \forall i\in J_{j}, \, \forall j, \nonumber\\[2pt] && z_{j}\geq 0, \quad \forall j, \nonumber\\[2pt] && \mu\geq 0, \, \nu\geq 0.\nonumber \end{eqnarray} $

(3.6)

This linear programming model can be solved by the simplex method or interior point method.

In this section, a simple numerical example computed by the Matlab software is considered to clearly illustrate the proposed robust DEA approach. Now we look at a DEA example which includes five DMUs, each DMU has three inputs and two outputs. Assume that the input data $x_{ij}$ is deterministic and the output data $y_{ij}$ is uncertain, which lies within the range $[\bar{y}_{ij}-d^L_{ij}, \bar{y}_{ij}+d^U_{ij}]$, where $\bar{y}_{ij}$ is the nominal output data. Table 1 shows the input data and the nominal output data.

For simplicity, we assume that $d^L_{ij}$ and $d^U_{ij}$ are 0.5 and 1 respectively, $J_j$ is equal to 2, so $\Gamma_j$ can vary between 0 and 2. Denote $\Gamma=(\Gamma_1, \Gamma_2, \cdots, \Gamma_5)$, we increase $\Gamma$ from 0 to 1.75, the efficiency scores of the DMUs are summarized in Figure 1. It is evident that the efficiency scores of the DMUs are decreasing as $\Gamma$ increasing, we can call this the price of robustness. As shown, compared with the others, DMU 2 is the most efficient decision making unit.

In this paper, we analyze the impact of output data uncertainty on the result of DEA and propose a robust model based on the newly developed robust optimization approaches. The conventional robust DEA model assumed the symmetric data uncertainty: the ranges of the uncertain elements are symmetrically bounded around their means. We relax the assumption of symmetric data uncertainty to make our robust DEA model more suitable for applications. We construct the robust DEA model for asymmetric data uncertainty and obtain an equivalent linear programming formulation. We implement the proposed model in a numerical example and the results indicate that considering the output data uncertainties when applying DEA approach is very important, and using robust DEA approach could be more reliable for efficiency evaluation and ranking in multiple criteria decision making problems.