[人工智能]python DEA: by-product NDDF DEA模型

点赞发Nature
关注中Science

by-product DEA模型取消了常用的DEA模型中的弱处置性（weak disposability）,可以用来模拟副产品与主产品相关性不强的情景。

如，在电力生产中，若安装了APCD，则空气污染物的产生与电力生产量就不具有强的相关性，应当采用by-product模型。

#!/usr/bin/env python
# coding: utf-8


import numpy as np
import pandas as pd
import pulp




class DEAProblem:
    """
    Help on class DEAProblem

    DEAProblem(inputs, outputs, bad_outs, weight_vector, directional_factor=None, returns='CRS',
                 in_weights=[0, None], out_weights=[0, None],badout_weights=[0, None])

    DEAProblem solves DEA model using directional distance function.

    Parameters:
    inputs: input data, DataFrame data
    outputs: output data, DataFrame data
    bad_outs: undesirable output data, DataFrame data
    weight_vector: weights for individual inputs and outputs. List data
    """

    def __init__(
        self,
        n_inputs,
        p_inputs,
        outputs,
        bad_outs,
        weight_vector,
        directional_factor=None,
        returns="CRS",
        in_weights=[0, None],
        out_weights=[0, None],
        badout_weights=[0, None],
    ):
        self.n_inputs = n_inputs
        self.p_inputs = p_inputs
        self.outputs = outputs
        self.bad_outs = bad_outs
        self.returns = returns
        self.weight_vector = (
            weight_vector  # weight vector in directional distance function
        )

        self.J, self.nI = self.n_inputs.shape  # no of DMUs, non-polluting inputs
        _, self.pI = self.p_inputs.shape  # no of polluting inputs
        _, self.R = self.outputs.shape  # no of outputs
        _, self.S = self.bad_outs.shape  # no of bad outputs
        self._ni = range(self.nI)  # iterate over non-polluting inputs
        self._pi = range(self.pI)  # iterable over polluting 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.gx1 = self.n_inputs
            self.gx2 = self.p_inputs
            self.gy = self.outputs
            self.gb = self.bad_outs
        else:
            self.gx1 = directional_factor[: self.nI]
            self.gx2 = directional_factor[: (self.nI + self.pI)]
            self.gy = directional_factor[
                (self.nI + self.pI) : ((self.nI + self.pI) + self.J)
            ]
            self.gy = directional_factor[((self.nI + self.pI) + 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 by-product technology model. Reference: doi.org/10.1111/deci.12421
        """
        # Set up pulp
        prob = pulp.LpProblem("".join(["DMU_", str(j0)]), pulp.LpMaximize)
        self.weights_1 = pulp.LpVariable.dicts(
            "Weight_non_polluting", (self._j), lowBound=self._in_weights[0]
        )  # define the weight for non-pulluting unit in the by-product technology model

        self.weights_2 = pulp.LpVariable.dicts(
            "Weight_polluting", (self._j), lowBound=self._in_weights[0]
        )  # define the weight for pulluting unit in the by-product technology model

        #         self.betax1 = pulp.LpVariable.dicts(
        #             "scalingFactor_nx", (self._ni), lowBound=0, upBound=1
        #         ) # scaling factor for non-polluting inputs
        #         self.betax2 = pulp.LpVariable.dicts(
        #             "scalingFactor_px", (self._pi), lowBound=0, upBound=1
        #         ) # scaling factor for polluting inputs

        self.beta = pulp.LpVariable("scalingFactor", lowBound=0)  # scaling factor

        #         self.betab = pulp.LpVariable.dicts(
        #             "scalingFactor_b", (self._s), lowBound=0, upBound=1
        #         ) # scaling factor for undesirable factor

        # Set up objective function
        prob += pulp.lpSum(self.beta)

        # Set up constraints
        for ni in self._ni:
            prob += (
                pulp.lpSum(
                    [
                        (self.weights_1[j0] * self.n_inputs.values[j0][ni])
                        for j0 in self._j
                    ]
                )
                <= self.n_inputs.values[j0][ni]
            )
        for pi in self._pi:
            prob += (
                pulp.lpSum(
                    [
                        (self.weights_1[j0] * self.p_inputs.values[j0][pi])
                        for j0 in self._j
                    ]
                )
                <= self.p_inputs.values[j0][pi]
            )
        # strong disposability for desirable output in non-pulluting process
        for r in self._r:
            prob += (
                pulp.lpSum(
                    [
                        (self.weights_1[j0] * self.outputs.values[j0][r])
                        for j0 in self._j
                    ]
                )
                >= self.outputs.values[j0][r] + self.beta * self.gy.values[j0][r]
            )

        for pi in self._pi:
            prob += (
                pulp.lpSum(
                    [
                        (self.weights_2[j0] * self.p_inputs.values[j0][pi])
                        for j0 in self._j
                    ]
                )
                >= self.p_inputs.values[j0][pi]
            )
        # strong disposability for undesirable output in polluting process
        for s in self._s:
            prob += (
                pulp.lpSum(
                    [
                        (self.weights_2[j0] * self.bad_outs.values[j0][s])
                        for j0 in self._j
                    ]
                )
                <= self.bad_outs.values[j0][s] - self.beta * self.gb.values[j0][s]
            )

        # Set returns to scale
        if self.returns == "VRS":
            prob += sum([self.weights_1[j] for j in self.weights_1]) == 1
            prob += sum([self.weights_2[j] for j in self.weights_2]) == 1

        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(pulp.PULP_CBC_CMD(msg=0))
            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)
#             for name, c in list(problem.constraints.items()):
#                 print(name, ":", c, "\t", c.pi, "\t\t", c.slack)
        return sol_status, sol_efficiency, sol_weights

使用的案例可以见我的Github

代码相关的文章工作还在进行中，代码有用的话，请考虑引用我的论文~

Cite my works

F. Yang, H. Lee. An innovative provincial CO2 emission quota allocation scheme for Chinese low-carbon transition. Technological forecast and social change. 2022. In press.

————————
我目前奥尔堡大学做博士后(我的google scholar, 我的Github)，关注能源转型过程中的环境、经济问题。

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

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