点赞发Nature
关注中Science

上一版非径向距离函数在这
普通非径向距离函数

现在很多研究在测度效率时要考虑非期望产出，如环境技术中CO2排放，银行业中的不良贷款等，这里我写了一个考虑非期望产出的非径向距离函数

import numpy as np
import pandas as pd
import pulp

class DEAProblem:
    def __init__(
        self,
        inputs,
        outputs,
        bad_outs,
        weight_vector,
        directional_factor=None,
        returns="CRS",
        disp="weak disposability",
        in_weights=[0, None],
        out_weights=[0, None],
        badout_weights=[0, None],
    ):
        self.inputs = inputs
        self.outputs = outputs
        self.bad_outs = bad_outs
        self.returns = returns
        self.weight_vector = (
            weight_vector  # weight vector in directional distance function
        )
        self.disp = disp

        self.J, self.I = self.inputs.shape  # no of DMUs, inputs
        _, self.R = self.outputs.shape  # no of outputs
        _, self.S = self.bad_outs.shape  # no of bad outputs
        self._i = range(self.I)  # inputs
        self._r = range(self.R)  # outputs
        self._s = range(self.S)  # bad_output
        self._j = range(self.J)  # DMUs
        if directional_factor == None:
            self.gx = self.inputs
            self.gy = self.outputs
            self.gb = self.bad_outs
        else:
            self.gx = directional_factor[: self.I]
            self.gy = directional_factor[self.I : (self.I + self.J)]
            self.gy = directional_factor[(self.I + self.J) :]

        self._in_weights = in_weights  # input weight restrictions
        self._out_weights = out_weights  # output weight restrictions
        self._badout_weights = badout_weights  # bad output weight restrictions

        # creates dictionary of pulp.LpProblem objects for the DMUs
        self.dmus = self._create_problems()

    def _create_problems(self):
        """
        Iterate over the DMU and create a dictionary of LP problems, one
        for each DMU.
        """

        dmu_dict = {}
        for j0 in self._j:
            dmu_dict[j0] = self._make_problem(j0)
        return dmu_dict

    def _make_problem(self, j0):
        """
        Create a pulp.LpProblem for a DMU.
        """

        # Set up pulp
        prob = pulp.LpProblem("".join(["DMU_", str(j0)]), pulp.LpMaximize)
        self.weights = pulp.LpVariable.dicts(
            "Weight", (self._j), lowBound=self._in_weights[0]
        )
        self.betax = pulp.LpVariable.dicts(
            "scalingFactor_x", (self._i), lowBound=0, upBound=1
        )

        self.betay = pulp.LpVariable.dicts("scalingFactor_y", (self._r), lowBound=0)

        self.betab = pulp.LpVariable.dicts(
            "scalingFactor_b", (self._s), lowBound=0, upBound=1
        )

        # Set returns to scale
        if self.returns == "VRS":
            prob += pulp.lpSum([self.weights[j] for j in self.weights]) == 1

        # Set up objective function
        prob += pulp.lpSum(
            [(self.weight_vector[i] * self.betax[i]) for i in self._i]
            + [(self.weight_vector[self.I + r] * self.betay[r]) for r in self._r]
            + [
                (self.weight_vector[self.I + self.R + s] * self.betab[s])
                for s in self._s
            ]
        )

        # Set up constraints
        for i in self._i:
            prob += (
                pulp.lpSum(
                    [(self.weights[j0] * self.inputs.values[j0][i]) for j0 in self._j]
                )
                <= self.inputs.values[j0][i] - self.betax[i] * self.gx.values[j0][i]
            )
        for r in self._r:
            prob += (
                pulp.lpSum(
                    [(self.weights[j0] * self.outputs.values[j0][r]) for j0 in self._j]
                )
                >= self.outputs.values[j0][r] + self.betay[r] * self.gy.values[j0][r]
            )

        if self.disp == "weak disposability":
            for s in self._s:  # weak disposability
                prob += (
                    pulp.lpSum(
                        [
                            (self.weights[j0] * self.bad_outs.values[j0][s])
                            for j0 in self._j
                        ]
                    )
                    == self.bad_outs.values[j0][s]
                    - self.betab[s] * self.gb.values[j0][s]
                )

        elif self.disp == "strong disposability":
            for s in self._s:  # strong disposability
                prob += (
                    pulp.lpSum(
                        [
                            (self.weights[j0] * self.bad_outs.values[j0][s])
                            for j0 in self._j
                        ]
                    )
                    >= self.bad_outs.values[j0][s]
                    - self.betab[s] * self.gb.values[j0][s]
                )

        return prob

    def solve(self):
        """
        Iterate over the dictionary of DMUs' problems, solve them, and collate
        the results into a pandas dataframe.
        """

        sol_status = {}
        sol_weights = {}
        sol_efficiency = {}

        for ind, problem in list(self.dmus.items()):
            problem.solve()
            sol_status[ind] = pulp.LpStatus[problem.status]
            sol_weights[ind] = {}
            for v in problem.variables():
                sol_weights[ind][v.name] = v.varValue
            sol_efficiency[ind] = pulp.value(problem.objective)
        return sol_status, sol_efficiency, sol_weights

简单说，这是一个考虑弱处置性的NDDF模型，input就是决策单元的输入，output是决策单元的输出，bad_outs是决策单元的非期望产出， weight_vector是计算距离函数的权重向量。

你们可以测试一下，有什么问题、有什么思路都可以联系我交流，我还在持续更新中。

希望以后能做成一个比较完整的开源模块给大家用。

举个例子

X = pd.DataFrame(
    np.array(
        [
            [20, 300],
            [30, 200],
            [40, 100],
            [20, 200],
            [10, 400],
            [11, 222],
            [12, 321],
            [14, 231],
        ]
    )
)
y = pd.DataFrame(np.array([[20], [30], [40], [30], [50], [21], [32], [42]]))
b = pd.DataFrame(np.array([[10], [20], [10], [10], [10], [12], [-2], [-1]]))
weight = [1 / 6, 1 / 6, 1 / 3, 1 / 3]
names = pd.DataFrame(
    ["Bratislava", "Zilina", "Kosice", "Presov", "Poprad", "ala", "ba", "ca"],
    columns=["DMU"],
)

solve = DEAProblem(X, y, b, weight, disp="weak disposability").solve()

status = pd.DataFrame.from_dict(solve[0], orient="index", columns=["status"])
efficiency = pd.DataFrame.from_dict(solve[1], orient="index", columns=["efficiency"])
weights = pd.DataFrame.from_dict(solve[2], orient="index")
results = pd.concat([names, status, efficiency, weights], axis=1)

print(results.round(decimals=4))

 DMU   status  efficiency  Weight_0  Weight_1  Weight_2  Weight_3  
0  Bratislava  Optimal      0.9182       0.0       0.0    0.0728       0.0   
1      Zilina  Optimal      0.5078       0.0       0.0    0.0830       0.0   
2      Kosice  Optimal      0.0000       0.0       0.0    1.0000       0.0   
3      Presov  Optimal      0.4663       0.0       0.0    0.0830       0.0   
4      Poprad  Optimal      0.0000       0.0       0.0    0.0000       0.0   
5         ala  Optimal      0.5646       0.0       0.0    0.0000       0.0   
6          ba  Optimal      0.4458       0.0       0.0    0.0000       0.0   
7          ca  Optimal      0.0000       0.0       0.0    0.0000       0.0   

Cite the work

Yang, F.; Choi, Y.: Lee, H. Life-cycle data envelopment analysis to measure efficiency and cost-effectiveness of environmental regulation in China’s transport sector. Ecological indicators 2021
Choi, Y.; Yang, F.; Lee, H. On the Unbalanced Atmospheric Environmental Performance of Major Cities in China. Sustainability 2020, 12, 5391. https://doi.org/10.3390/su12135391

————————update 2020/12/4——————————
更新了VRS的计算方法
from

        if self.returns == "VRS":
            prob += pulp.lpSum([weight for weight in self.weights]) == 1

to

        if self.returns == "VRS":
            prob += pulp.lpSum([self.weights[j] for j in self.weights]) == 1

————————
我是仁荷大学的经管博士生(我的google scholar, 我的Github)，关注能源转型过程中的环境、经济问题。

专注于分享利用python科研的技巧，欢迎一起交流、学习、合作。

关于我的博客内容、其他相关的研究问题，有问题可以下方👇评论、或私信我~