# This file contains functions for the Enhanced Russell Graph Slack Based Measure DEA model
"""
EnhanceRussellGraphDEAModel
An data structure representing an Enhanced Russell Graph Slack Based Measure DEA model.
"""
struct EnhancedRussellGraphDEAModel <: AbstractTechnicalDEAModel
n::Int64
m::Int64
s::Int64
rts::Symbol
dmunames::Union{Vector{String},Nothing}
eff::Vector
beta::Vector
tX::Matrix
tY::Matrix
slackX::Matrix
slackY::Matrix
lambda::SparseMatrixCSC{Float64, Int64}
Xtarget::Matrix
Ytarget::Matrix
end
"""
deaerg(X, Y)
Compute data envelopment analysis Enhanced Russell Graph Slack Based Measure for inputs `X` and outputs `Y`.
# Optional Arguments
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.
# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];
julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353] ;
julia> deaerg(X, Y, rts = :VRS)
Enhanced Russell Graph Slack Based Measure DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Orientation = Graph; Returns to Scale = VRS
───────────────────────────────────────
efficiency beta slackX1 slackY1
───────────────────────────────────────
1 1.0 1.0 0.0 0.0
2 1.0 1.0 0.0 0.0
3 1.0 1.0 0.0 0.0
4 1.0 1.0 0.0 0.0
5 0.4 0.6 2.0 2.0
6 0.47619 1.0 7.33333 0.0
7 0.857143 1.0 2.0 0.0
8 0.2 0.4706 5.412 2.647
───────────────────────────────────────
```
"""
function deaerg(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
rts::Symbol = :CRS, slack::Bool = true,
Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
names::Union{Vector{String},Nothing} = nothing,
optimizer::Union{DEAOptimizer,Nothing} = nothing)::EnhancedRussellGraphDEAModel
# Check parameters
nx, m = size(X, 1), size(X, 2)
ny, s = size(Y, 1), size(Y, 2)
if Xref === nothing Xref = X end
if Yref === nothing Yref = Y end
nrefx, mref = size(Xref, 1), size(Xref, 2)
nrefy, sref = size(Yref, 1), size(Yref, 2)
if nx != ny
error("number of observations is different in inputs and outputs")
end
if nrefx != nrefy
error("number of observations is different in inputs reference set and ouputs reference set")
end
if m != mref
error("number of inputs in evaluation set and reference set is different")
end
if s != sref
error("number of outputs in evaluation set and reference set is different")
end
# Default optimizer
if optimizer === nothing
optimizer = DEAOptimizer(GLPK.Optimizer)
end
# Compute efficiency for each DMU
n = nx
nref = nrefx
effi = zeros(n)
betai = zeros(n)
tXi = zeros(n, m)
tYi = zeros(n, s)
mui = spzeros(n, nref)
for i=1:n
# Value of inputs and outputs to evaluate
x0 = X[i,:]
y0 = Y[i,:]
# Create the optimization model
deamodel = newdeamodel(optimizer)
@variable(deamodel, beta >= 0)
@variable(deamodel, tX[1:m] >= 0)
@variable(deamodel, tY[1:s] >= 0)
@variable(deamodel, mu[1:nref] >= 0)
@objective(deamodel, Min, beta - 1 / m * sum(tX[t] / x0[t] for t in 1:m))
@constraint(deamodel, beta + 1 / s * sum(tY[t] / y0[t] for t in 1:s) == 1)
@constraint(deamodel, [j in 1:m], sum(Xref[t,j] * mu[t] for t in 1:nref) == beta * x0[j] - tX[j])
@constraint(deamodel, [j in 1:s], sum(Yref[t,j] * mu[t] for t in 1:nref) == beta * y0[j] + tY[j])
# Add return to scale constraints
if rts == :CRS
# No contraint to add for constant returns to scale
elseif rts == :VRS
@constraint(deamodel, sum(mu) == beta)
else
error("Invalid returns to scale $rts. Returns to scale should be :CRS or :VRS")
end
# Optimize and return results
JuMP.optimize!(deamodel)
effi[i] = JuMP.objective_value(deamodel)
betai[i] = JuMP.value(beta)
tXi[i,:] = JuMP.value.(tX)
tYi[i,:] = JuMP.value.(tY)
mui[i,:] = JuMP.value.(mu)
# Check termination status
if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
@warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
end
end
slackX = tXi ./ betai
slackY = tYi ./ betai
lambdaeff = mui ./ betai
# Get X and Y targets
Xtarget = X - slackX
Ytarget = Y + slackY
return EnhancedRussellGraphDEAModel(n, m, s, rts, names, effi, betai, tXi, tYi, slackX, slackY, lambdaeff, Xtarget, Ytarget)
end
function Base.show(io::IO, x::EnhancedRussellGraphDEAModel)
compact = get(io, :compact, false)
n = nobs(x)
m = ninputs(x)
s = noutputs(x)
dmunames = names(x)
eff = efficiency(x)
beta = x.beta
slackX = slacks(x, :X)
slackY = slacks(x, :Y)
hasslacks = ! isempty(slackX)
if !compact
print(io, "Enhanced Russell Graph Slack Based Measure DEA Model \n")
print(io, "DMUs = ", n)
print(io, "; Inputs = ", m)
print(io, "; Outputs = ", s)
print(io, "\n")
print(io, "Orientation = Graph")
print(io, "; Returns to Scale = ", string(x.rts))
print(io, "\n")
show(io, CoefTable(hcat(eff, beta, slackX, slackY), ["efficiency"; "beta"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
end
end
function efficiency(model::EnhancedRussellGraphDEAModel, type::Symbol)::Vector
if type == :beta return model.beta end
error(typeof(model), " has no efficiency type ", string(type))
end