带有非期望产出的SBM模型(slacks-based measure)可用于处理多个投入和产出变量的效率测度问题，可用于工业经济绿色发展效率的分析。在绿色发展分析中，需要同时考虑经济收益（期望产出）和负面环境效应（非期望产出）的问题，一方面要提高经济收益，另一方面需要减少污染排放，是投入、期望产出和非期望产出三方权衡的问题。¹

这个帖子记录了和非期望产出SBM模型苦苦斗争一个星期的学习成果。

数学原理

DEA模型/CCR模型

DEA

SBM模型是数据包络分析(Data Envelopment Analysis, DEA)的一种。在DEA模型中，假定存在$n$个可比决策单元(decision making units, DMU)，记为$DM{U}_j(j=1,2,\dots ,n)$；

每个$DMU$有$m$种投入，记为$X_i(i=1,2,\dots ,m)$，每种投入的权重为$v_i$；

有$q$种产出，记为$Y_r(r=1,2,\dots ,q)$，每种产出的权重为$u_r$.

对第$k$个可比决策单元而言，其投入产出比（技术效率）为

$$H_k=\frac{\sum _{r=1}^q{u}_r{Y}_{rk}}{\sum _{i=1}^m{v}_i{X}_{ik}}，$$

限定范围为$H_k\in [0,1]$. 现要求在所有可比决策单元的效率均不超过$1$的条件下，调节权重$[v_1,v_2,\dots ,v_m]$和$[u_1,u_2,\dots ,u_q]$，使得被评价的单元$DM{U}_k$效率值最大化。即如下规划模型： $$\begin{array}{rl}max& z=\frac{\sum _{r=1}^q{u}_r{Y}_{rk}}{\sum _{i=1}^m{v}_i{X}_{ik}}\\ s.t.& \frac{\sum _{r=1}^q{u}_r{Y}_{rj}}{\sum _{i=1}^m{v}_i{X}_{ij}}\le 1(v_i\ge 0,u_r\ge 0)\\ & i=1,2,\dots ,m;r=1,2,\dots ,q;j=1,2,\dots ,n\end{array}$$

该模型所确定的权重系数组合是对于$DM{U}_k$最有利的。由于$\sum _{i=1}^m{v}_i{X}_{ik}>0$，模型的规划条件即为 $$s.t.\sum _{r=1}^q{u}_r{Y}_{rj}-\sum _{i=1}^m{v}_i{X}_{ij}\le 0;$$

在上述线性规划模型中，所要求的是投入和产出项目的指标权重，在原始意义上就是$\overrightarrow{u}$和$\overrightarrow{v}$。模型指定了一个可比决策单元$DM{U}_k$，在规划过程中力求这个单元最终的效率评价达到最高。实际应用中，如何选择这个特殊的$DMU$必然会成为问题的关键。或许需要根据决策者的感性判断（或其他手段辅助），选定一个最优的决策单元，或者人为设置一个对照的“理想型”，将其视为效率最高的情形。

投入导向CCR

令$t={\left(\sum _{i=1}^m{v}_i{X}_{ik}\right)}^{-1}$，则模型的目标函数为

$$maxz=t\sum _{r=1}^q{u}_r{Y}_{rk}=\sum _{r=1}^{q}t{u}_r{Y}_{rk},$$

令${\mu }_r=t{u}_r,{\nu }_i=t{v}_i$，根据$t$的定义可知$\sum _{i=1}^m{\nu }_i{X}_{ik}=1$；将规划条件两边同乘$t$，得到

$$s.t.\sum _{r=1}^q{\mu }_r{Y}_{rj}-\sum _{i=1}^m{\nu }_i{X}_{ij}\le 0.$$

从而将非线性规划模型转换为线性规划模型：

$$\begin{array}{rl}max& z=\sum _{r=1}^q{\mu }_r{Y}_{rk}\\ s.t.& \sum _{r=1}^q{\mu }_r{Y}_{rj}-\sum _{i=1}^m{\nu }_i{X}_{ij}\le 0({\mu }_r\ge 0,{\nu }_i\ge 0)\\ & \sum _{i=1}^m{\nu }_i{X}_{ik}=1\\ & i=1,2,\dots ,m;r=1,2,\dots ,q;j=1,2,\dots ,n\end{array}$$

其对偶模型为（对偶模型求解过程见文末）

$$\begin{array}{rl}min& \theta =\sum _{j=1}^n\frac{X_{mj}}{X_{mk}}{\lambda }_j+{\lambda }_{n+1}\\ s.t.& \sum _{j=1}^n{\lambda }_j{X}_{ij}\le (\theta +{\lambda }_{n+1})X_{ik}({\lambda }_j\ge 0)\\ & \sum _{j=1}^n{\lambda }_j{Y}_{rj}\ge Y_{rk}\\ & i=1,2,\dots ,m;r=1,2,\dots ,q;j=1,2,\dots ,n\end{array}$$

以上处理后得到的是投入导向的DEA模型。

这里与常见的处理方式有点差异，因为在转化为标准形式的时候规划模型的自由度减1。

产出导向CCR

同理，若在降低自变量自由度时采用$t={\left(\sum _{r=1}^q{u}_r{Y}_{rk}\right)}^{-1}$，则可演算得到产出导向的DEA模型：

$$\begin{array}{rl}min& z=\sum _{i=1}^m{\nu }_r{X}_{ik}\\ s.t.& \sum _{r=1}^q{\mu }_r{Y}_{rj}-\sum _{i=1}^m{\nu }_i{X}_{ij}\le 0({\mu }_r\ge 0,{\nu }_i\ge 0)\\ & \sum _{r=1}^q{\mu }_r{Y}_{rk}=1\\ & i=1,2,\dots ,m;r=1,2,\dots ,q;j=1,2,\dots ,n\end{array}$$

其对偶模型为

$$\begin{array}{rl}min& \varphi =\sum _{j=1}^n\frac{Y_{qj}}{Y_{qk}}{\lambda }_j+{\lambda }_{n+1}\\ s.t.& \sum _{j=1}^n{\lambda }_j{Y}_{rj}\le (\theta +{\lambda }_{n+1})Y_{rk}({\lambda }_j\ge 0)\\ & \sum _{j=1}^n{\lambda }_j{X}_{ij}\ge X_{rk}\\ & i=1,2,\dots ,m;r=1,2,\dots ,q;j=1,2,\dots ,n\end{array}$$

SBM模型

不带非期望产出的情形

在DEA模型中，线性规划的最终目的是使得$H_k$尽可能大，即$DM{U}_k$的产出效率尽可能高。虽然考虑到了多种产出的情形，但是每一种产出的权重都是正的($u_r\ge 0$) ，所以依然是要尽可能使得每一种产出都越高越好。当产出中存在非期望内容时，我们希望这些非期望的产出尽可能低，同时期望产出尽可能高。Kaoru Tone提出的SBM模型(slacks-based model)²可以处理这类非期望产出的问题。

我们仍然以$X=(x_{ij})\in R^{m\times n}$表示实际投入，$Y=(y_{rj})\in R^{q\times n}$表示实际产出。在不带非期望产出的情况下，有可能存在的生产情形必然包括投入更多、产出更低的情况（即效率更低的情况总是能够实现的）：

$$P=\{(\overrightarrow{x},\overrightarrow{y})|\overrightarrow{x}\ge X\overrightarrow{\lambda },\overrightarrow{y}\le Y\overrightarrow{\lambda },\overrightarrow{\lambda }\ge \overrightarrow{0}\}$$

其中$\overrightarrow{\lambda }$是$R^n$上的一个非负向量。

$X$矩阵的每一列代表了一个可比决策单元，当$X$右乘一个$n$维向量后得到的是一个$m$维向量，$X\overrightarrow{\lambda }={\left[\sum _{j=1}^n{x}_{ij}{\lambda }_j\right]}_{m\times 1}$，可以视为由所有$DMU$按照$\overrightarrow{\lambda }$所规定的配比混合而成的新的生产投入组合。同理，$Y\overrightarrow{\lambda }$代表了所有$DMU$按照$\overrightarrow{\lambda }$所规定的配比混合而成的新的产出组合。$(X\overrightarrow{\lambda },Y\overrightarrow{\lambda })$就是所有实际生产情形线性混合而成的“虚拟生产情形”（基于规模报酬不变）。在这个“虚拟生产情形”的基础之上，更加浪费的生产情形必然是可行的，故采用上述形式定义有可能的生产情形。

在这种定义下，默认了效率最高的生产情形暗含在实际生产情形通过线性组合所形成的参数空间范围之中。即我们不考虑生产效率远远超出现有技术水平的情形。

具体地，用$({\overrightarrow{x}}_0,{\overrightarrow{y}}_0)$来描述某一个$DM{U}_0$（从实际生产情形$X,Y$中而来）： $$\begin{array}{c}{\overrightarrow{x}}_0=X\overrightarrow{\lambda }+s{s}^{-}{\overrightarrow{y}}_0=Y\overrightarrow{\lambda }-s{s}^{+}\end{array}$$ 其中$s{s}^{-}\ge 0,s{s}^{+}\ge 0$分别称为投入过量(input excess)和产出不足(output shortfall)，统一称为和投入和产出的松弛向量(slacks)。据此定义指标：

$$\rho =\frac{1-\frac{1}{m}\sum _{i=1}^m\left(s_i^{-}/x_{i0}\right)}{1+\frac{1}{q}\sum _{r=1}^q\left(s_r^{+}/y_{r0}\right)}$$

该指标满足以下特性：

1. 每个$DMU$都有一个不同的$\rho$
2. 随着投入和产出松弛向量$s{s}^{-},s{s}^{+}$而单调变化
3. 取值范围为$\rho \in (0,1]$

为了表达$({\overrightarrow{x}}_0,{\overrightarrow{y}}_0)$的效率，以$\rho$指标构建SBM模型：

$$\begin{array}{rl}min& \rho =\frac{1-\frac{1}{m}\sum _{i=1}^m\left(s_i^{-}/x_{i0}\right)}{1+\frac{1}{q}\sum _{r=1}^q\left(s_r^{+}/y_{r0}\right)}\\ s.t.& {\overrightarrow{x}}_0=X\overrightarrow{\lambda }+s{s}^{-}\\ & {\overrightarrow{y}}_0=Y\overrightarrow{\lambda }-s{s}^{+}\\ & \overrightarrow{\lambda }\ge 0,s{s}^{-}\ge 0,s{s}^{+}\ge 0\end{array}$$

当$DM{U}_0$达到效率最优(SBM-efficient)时，满足 $${\rho }^{\ast }=1,s{s}^{-}=0,s{s}^{+}=0$$

如何理解$DMU$效率最优等同于松弛变量均为零？在作者的理解中，认为SBM-efficient意味着该决策单元的任何优化都不会引起投入过量和产出不足。

若有唯一的$\overrightarrow{\lambda }$，则确定唯一的$(s{s}^{-},s{s}^{+})\sim ({\overrightarrow{x}}_0,{\overrightarrow{y}}_0)$关系。规划自变量为$(s{s}^{-},s{s}^{+})$，目的是使得$\rho$尽可能小，即$s{s}^{-}$和$s{s}^{+}$都要尽可能大。这可以看成是一种对当下生产效率不信任的态度。

当处于效率最优状态时，生产投入的要素配比恰到好处，此时进行任何调节都无法使得效率进一步提升。

再将$\overrightarrow{\lambda }$也视为自变量，即尝试不同组合情况下的$DM{U}_0$，通过调整$\overrightarrow{\lambda }$可以进一步优化目标函数。

**Charnes-Cooper transformation: **将目标函数的分子和分母都乘$t$，使得分母变为$1$，即

$$t={[1+\frac{1}{q}\sum _{r=1}^q\left(\frac{s_r^{+}}{y_{r0}}\right)]}^{-1},$$

从而目标函数变为

$$min\rho =t-\frac{1}{m}\sum _{i=1}^m\frac{t{s}_i^{-}}{x_{i0}}.$$

定义$S{S}^{-}=ts{s}^{-},S{S}^{+}=ts{s}^{+},\mathrm{\Lambda }\mathrm{\Lambda }=t\overrightarrow{\lambda }$，则原非线性规划中的约束条件转换为 $$\begin{array}{rl}s.t.& 1=t+\frac{1}{q}\sum _{r=1}^q\frac{S_r^{+}}{y_{r0}}\\ & t{\overrightarrow{x}}_0=X\mathrm{\Lambda }\mathrm{\Lambda }+S{S}^{-}\\ & t{\overrightarrow{y}}_0=Y\mathrm{\Lambda }\mathrm{\Lambda }-S{S}^{+}\\ & \mathrm{\Lambda }\ge 0,S{S}^{-}\ge 0,S{S}^{+}\ge 0,t\ge 0\end{array}$$ 最后完整地写一下转换后的线性规划（多项式形式）：

$$\begin{array}{rl}min& \rho =t-\frac{1}{m}\sum _{i=1}^m\frac{t{s}_i^{-}}{x_{i0}}\\ s.t.& 1=t+\frac{1}{q}\sum _{r=1}^q\frac{S_r^{+}}{y_{r0}}\\ & t{x}_{i0}=\sum _{j=1}^n{x}_{ij}{\mathrm{\Lambda }}_j+S_i^{-}\\ & t{y}_{r0}=\sum _{j=1}^n{y}_{rj}{\mathrm{\Lambda }}_j-S_r^{+}\\ & {\mathrm{\Lambda }}_j\ge 0,S_i^{-}\ge 0,S_r^{+}\ge 0,t\ge 0\\ & i=1,2,\dots ,m;r=1,2,\dots ,q;j=1,2,\dots ,n\end{array}$$

通过线性规划的可行解$({\rho }^{\ast },{\overrightarrow{\mathrm{\Lambda }}}^{\ast },t^{\ast },S{S}^{-\ast },S{S}^{+\ast })$，可以解出原型的可行解$({\rho }^{\ast },{\overrightarrow{\lambda }}^{\ast },s{s}^{-\ast },s{s}^{+\ast })$。

若${\rho }^{\ast }=1$，则称该$DMU$为SBM-efficient，此时$s{s}^{-\ast }=0,s{s}^{+\ast }=0$； 若${\rho }^{\ast }<1$，则称该$DMU$为SBM-inefficient，对于一个SBM-inefficient的决策单元（假设为$DM{U}_k$），其投入和产出满足：

$$\begin{array}{c}{\overrightarrow{x}}_k=X{\overrightarrow{\lambda }}^{\ast }+s{s}^{-\ast }\\ {\overrightarrow{y}}_k=Y{\overrightarrow{\lambda }}^{\ast }-s{s}^{+\ast }\end{array}$$

由于得到了投入过量和产出不足的具体数值，因此我们可以对该可比决策单元进行优化，称为SBM投射(SBM-projection)：

$$\begin{array}{c}{\overrightarrow{x}}_k\leftarrow {\overrightarrow{x}}_k-s{s}^{-\ast }\\ {\overrightarrow{y}}_k\leftarrow {\overrightarrow{y}}_k+s{s}^{+\ast }\end{array}$$

带有非期望产出的情形

如果我们考虑到非期望产出的情形³，这些非期望产出在低效情况下的情形并非产出不足，而是“产出过多”，因此需要引入新的松弛向量$s{s}^{{+}^{\prime }}$以规范非期望产出。记$Y$为实际期望产出，$Z$为实际非期望产出，分别有$q_1,q_2$个。生产可能性中也应区分期望产出和非期望产出，即$({\overrightarrow{x}}_0,{\overrightarrow{y}}_0,{\overrightarrow{z}}_0)$，其中 $${\overrightarrow{z}}_0=Z\overrightarrow{\lambda }+s{s}^{{+}^{\prime }}.$$ 同理以$\rho$指标构建SBM模型：

$$\begin{array}{rl}min& \rho =\frac{1-\frac{1}{m}\sum _{i=1}^m\left(s_i^{-}/x_{i0}\right)}{1+\frac{1}{q_1+q_2}[\sum _{r_1=1}^{q_1}\left(s_r^{+}/y_{r0}\right)+\sum _{r=1}^{q_2}\left(s_r^{{+}^{\prime }}/z_{r0}\right)]}\\ s.t.& {\overrightarrow{x}}_0=X\overrightarrow{\lambda }+s{s}^{-}\\ & {\overrightarrow{y}}_0=Y\overrightarrow{\lambda }-s{s}^{+}\\ & {\overrightarrow{z}}_0=Z\overrightarrow{\lambda }+s{s}^{{+}^{\prime }}\\ & \overrightarrow{\lambda }\ge 0,s{s}^{-}\ge 0,s{s}^{+}\ge 0\end{array}$$

当某一个$DMU$达到效率最优时，

$${\rho }^{\ast }=1,s{s}^{-}=0,s{s}^{+}=0,s{s}^{{+}^{\prime }}=0.$$

通过Charnes-Cooper转换，得到线性规划模型： $$\begin{array}{rl}min& \rho =t-\frac{1}{m}\sum _{i=1}^m\frac{S_i^{-}}{x_{i0}}\\ s.t.& 1=t+\frac{1}{q_1+q_2}(\sum _{r=1}^{q_1}\frac{S_r^{+}}{y_{r0}}+\sum _{r=1}^{q_2}\frac{S_r^{{+}^{\prime }}}{y_{r0}})\\ & t{x}_{i0}=\sum _{j=1}^n{x}_{ij}{\mathrm{\Lambda }}_j+S_i^{-}(i=1,2,\dots ,m)\\ & t{y}_{r0}=\sum _{j=1}^n{y}_{rj}{\mathrm{\Lambda }}_j-S_r^{+}(r=1,2,\dots ,q_1)\\ & t{z}_{r0}=\sum _{j=1}^n{z}_{rj}{\mathrm{\Lambda }}_j-S_r^{+}(r=1,2,\dots ,q_2)\\ & {\mathrm{\Lambda }}_j\ge 0,S_i^{-}\ge 0,S_r^{+}\ge 0,t\ge 0\end{array}$$ 由线性规划的可行解$({\rho }^{\ast },{\overrightarrow{\mathrm{\Lambda }}}^{\ast },t^{\ast },S{S}^{-\ast },S{S}^{+\ast },S{S}^{{+}^{\prime }\ast })$，可以解出原始模型的可行解$({\rho }^{\ast },{\overrightarrow{\lambda }}^{\ast },s{s}^{-\ast },s{s}^{+\ast },s{s}^{{+}^{\prime }\ast })$，

$$\lambda {\lambda }^{\ast }=\frac{\Lambda {\Lambda }^{\ast }}{t^{\ast }},s{s}^{-\ast }=\frac{S{S}^{-\ast }}{t^{\ast }},s{s}^{+\ast }=\frac{S{S}^{+\ast }}{t^{\ast }},s{s}^{{+}^{\prime }\ast }=\frac{S{S}^{{+}^{\prime }\ast }}{t^{\ast }}$$

若${\rho }^{\ast }=1$，则称该$DMU$为SBM-efficient，此时$s{s}^{-\ast }=0,s{s}^{+\ast }=0,s{s}^{{+}^{\prime }\ast }=0$；

若${\rho }^{\ast }<1$，则称该$DMU$为SBM-inefficient，对于一个SBM-inefficient的决策单元（假设为$DM{U}_k$），其投入和产出满足：

$$\begin{array}{c}{\overrightarrow{x}}_k=X{\overrightarrow{\lambda }}^{\ast }+s{s}^{-\ast }\\ {\overrightarrow{y}}_k=Y{\overrightarrow{\lambda }}^{\ast }-s{s}^{+\ast }\\ {\overrightarrow{z}}_k=Z{\overrightarrow{\lambda }}^{\ast }+s{s}^{{+}^{\prime }\ast }\end{array}$$

由于得到了投入过量和产出不足的具体数值，因此我们可以对该可比决策单元进行优化，称为SBM投射(SBM-projection)： $$\begin{array}{c}{\overrightarrow{x}}_k\leftarrow {\overrightarrow{x}}_k-s{s}^{-\ast }\\ {\overrightarrow{y}}_k\leftarrow {\overrightarrow{y}}_k+s{s}^{+\ast }\\ {\overrightarrow{z}}_k\leftarrow {\overrightarrow{z}}_k-s{s}^{{+}^{\prime }\ast }\end{array}$$

Python代码实现

为了便于编程，我们将以上规划问题写成矩阵形式， $$\begin{array}{rl}min& \rho =C\cdot XX\\ s.t.& A\cdot XX=b\\ & XX\ge 0\end{array}$$ 其中自变量为 $$XX={[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n,t,S_1^{-},S_2^{-},\dots ,S_m^{-},S_1^{+},S_2^{+},\dots ,S_{q_1}^{+},S_1^{{+}^{\prime }},S_2^{{+}^{\prime }},\dots ,S_{q_2}^{{+}^{\prime }}]}^{\mathrm{T}},$$ 目标函数系数矩阵为 $$C=[0,0,\dots ,0,1,-\frac{1}{m{x}_{10}},-\frac{1}{m{x}_{20}},\dots ,-\frac{1}{m{x}_{m0}},0,\dots ,0],$$ 约束条件方程系数为 $$\begin{array}{c}A=\left[\begin{array}{ccccccccccccc}0& \cdots & 0& 1& 0& \cdots & 0& {\alpha }_1& \cdots & {\alpha }_{q_1}& {\beta }_1& \dots & {\beta }_{q_2}\\ \\ \\ x_{11}& \cdots & x_{1n}& -x_{10}& 1& \cdots & 0& 0& \cdots & 0& 0& \cdots & 0\\ & ⋮& & & & \ddots & & & ⋮& & & ⋮\\ x_{m1}& \cdots & x_{mn}& -x_{m0}& 0& \cdots & 1& 0& \cdots & 0& 0& \cdots & 0\\ \\ \\ y_{11}& \cdots & y_{1n}& -y_{10}& 0& \cdots & 0& -1& \cdots & 0& 0& \cdots & 0\\ & ⋮& & & & ⋮& & & \ddots & & & ⋮\\ y_{q_{1}1}& \cdots & y_{q_{1}n}& -y_{q_{1}0}& 0& \cdots & 0& 0& \cdots & -1& 0& \cdots & 0\\ \\ \\ z_{11}& \cdots & z_{1n}& -z_{10}& 0& \cdots & 0& 0& \cdots & 0& 1& \cdots & 0\\ & ⋮& & & & ⋮& & & ⋮& & & \ddots \\ z_{q_{2}1}& \cdots & z_{q_{2}n}& -z_{q_{2}0}& 0& \cdots & 0& 0& \cdots & 0& 0& \cdots & 1\end{array}\right],\\ {\alpha }_r=\frac{1}{(q_1+q_2)y_{r0}}(r=1,2,\dots ,q_1),\\ {\beta }_r=\frac{1}{(q_1+q_2)z_{r0}}(r=1,2,\dots ,q_2),\end{array}$$

约束条件常数项为 $$b={[1,0,0,\dots ,0]}^{\mathrm{T}}$$

代码结构来自wonder1322的博客

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

import pandas as pd
import scipy.optimize as op
import numpy as np

def sbmeff2(input_variable, desirable_output, undesirable_output, dmu, data, method = 'revised simplex'):
    """用于求解sbm模型
    
    Parameters:
    -----------
    input_variable:
        投入[v1,v2,v3,...]
    desirable_output:
        期望产出[v1,v2,v3,...]
    undesirable_output:
        非期望产出[v1,v2,v3,...] 
    dmu:
        决策单元
    data:
        主数据
    method:
        求解方法.默认'revised simplex',可选'interior-point'

	Return:
	------
	res : DataFrame
		结果数据框[dmu	TE	slack...]
    """
 
    res = pd.DataFrame(columns = ['dmu','TE'], index = data.index)
    res['dmu'] = data[dmu]
    ## lambda有dmu个数个，S有变量个数个
    dmu_counts = data.shape[0]
    ## 投入个数
    m = len(input_variable)
    ## 期望产出个数
    s1 = len(desirable_output)
    ## 非期望产出个数
    s2 = len(undesirable_output)
    ## x[:dmu_counts] 为lambda
    ## x[dmu_counts : dmu_counts+1] 为 t
    ## x[dmu_counts+1 : dmu_counts + m + 1] 为投入slack
    ## x[dmu_counts+ 1 + m : dmu_counts + 1 + m + s1] 为期望产出slack
    ## x[dmu_counts + 1 + m + s1 :] 为非期望产出slack
    total = dmu_counts + m + s1 + s2 + 1
    cols = input_variable + desirable_output + undesirable_output
    newcols = []
    for j in cols:
        newcols.append(j+'_slack')
        res[j+'_slack'] = np.nan
        
    for i in range(dmu_counts):
        ## 优化目标：目标函数的系数矩阵
        c = [0]*dmu_counts + [1] + list(-1 / (m*data.loc[i, input_variable])) + [0] * (s1+s2)
        
        ## 约束条件：约束方程的系数矩阵
        A_eq = [[0]*dmu_counts + [1] + [0]*m +
                list(1/((s1+s2)*data.loc[i, desirable_output])) +
                list(1/((s1+s2)*data.loc[i, undesirable_output]))]
        
        ## 约束条件（1）：投入松弛变量为正
        for j1 in range(m):
            list1 = [0] * m
            list1[j1] = 1
            eq1 = list(data[input_variable[j1]]) + [-data.loc[i ,input_variable[j1]]] + list1 + [0]*(s1+s2)
            A_eq.append(eq1)
        ## 约束条件（2）：期望产出松弛变量为正
        for j2 in range(s1):
            list2 = [0] * s1
            list2[j2] = -1
            eq2 = list(data[desirable_output[j2]]) + [-data.loc[i, desirable_output[j2]]] + [0]*m + list2 + [0]*s2
            A_eq.append(eq2)
        ## 约束条件（3）：非期望产出松弛变量为正
        for j3 in range(s2):
            list3 = [0] * s2
            list3[j3] = 1
            eq3 = list(data[undesirable_output[j3]]) + [-data.loc[i, undesirable_output[j3]]] + [0]*(m+s1) + list3
            A_eq.append(eq3)         
        ## 约束条件：常数向量
        b_eq = [1] + [0]*(m+s1+s2)               
        bounds = [(0, None)]*total ## 约束边界为零
        ## 求解
        op1 = op.linprog(c = c,A_eq=A_eq,b_eq=b_eq,bounds=bounds,method = method)
        res.loc[i, 'TE'] = op1.fun
        res.loc[i, newcols] = op1.x[dmu_counts+1 :]
    return res

补充：求解对偶模型的过程

线性规划的矩阵形式

学习SBM模型的时候遇到了对偶规划这一步骤，之前没有接触过，在B站上找了一个教程对照着推导了一遍。

原始的线性规划模型，自变量是$\overrightarrow{\mu }=[{\mu }_1,{\mu }_2,\dots ,{\mu }_q],\overrightarrow{\nu }=[{\nu }_1,{\nu }_2,\dots ,{\nu }_m]$。

注意：虽然规划模型中有$X$、$Y$，但是这是代表可比决策单元的投入量和产出量的，是已知的数据，所以并不是变量。 $$\begin{array}{rl}max& z=\sum _{r=1}^q{\mu }_r{Y}_{rk}\\ s.t.& \sum _{r=1}^q{\mu }_r{Y}_{rj}-\sum _{i=1}^m{\nu }_i{X}_{ij}\le 0({\mu }_r\ge 0,{\nu }_i\ge 0)\\ & \sum _{i=1}^m{\nu }_i{X}_{ik}=1\\ & i=1,2,\dots ,m;r=1,2,\dots ,q;j=1,2,\dots ,n\end{array}$$ 先转换为矩阵形式，自变量向量记为$XX$， $$\begin{array}{rl}max& Z=C\cdot XX\\ s.t.& A\cdot XX\le b\\ & XX\ge 0\end{array}$$ 其中由于等式$\sum _{i=1}^m{\nu }_i{X}_{ik}=1$，导致自变量的自由度减1，这里剔除${\nu }_m$，即 $$XX={[{\mu }_1,{\mu }_2,\dots ,{\mu }_q,{\nu }_1,{\nu }_2,\dots ,{\nu }_{m-1}]}^{\mathrm{T}}.$$ 目标函数的系数向量为 $$C=[Y_{1k},Y_{2k},\dots ,Y_{qk},0,0,\dots ,0],$$ 其中有$m-1$个$0$。

除了要满足$XX\ge 0$以外，还应保证${\nu }_m\ge 0$，即 $${\nu }_m=\frac{1-\sum _{i=1}^{m-1}{\nu }_i{X}_{ik}}{X_{mk}}\ge 0\iff \sum _{i=1}^{m-1}{\nu }_i{X}_{ik}\le 1.$$ 对于模型中的主要约束条件，消去${\nu }_m$后得到 $$\sum _{r=1}^q{\mu }_r{Y}_{rj}+\sum _{i=1}^{m-1}{\nu }_i(\frac{X_{ik}{X}_{mj}}{X_{mk}}-X_{ij})\le \frac{X_{mj}}{X_{mk}},j=1,2,\dots ,n.$$ 以上两式整合得到$A\cdot XX\le b$，则

$$\begin{array}{c}{A}_{(n+1)\times (q+m-1)}=\left[\begin{array}{cccccccc}{Y}_{11}& Y_{21}& \cdots & Y_{q1}& {\alpha }_{11}& {\alpha }_{21}& \cdots & {\alpha }_{m-1,1}\\ Y_{12}& Y_{22}& \cdots & Y_{q2}& {\alpha }_{12}& {\alpha }_{22}& \cdots & {\alpha }_{m-1,2}\\ & & ⋮& & & & ⋮\\ Y_{1n}& Y_{2n}& \cdots & Y_{qn}& {\alpha }_{1n}& {\alpha }_{2n}& \cdots & {\alpha }_{m-1,n}\\ 0& 0& \cdots & 0& X_{1k}& X_{2k}& \cdots & X_{m-1,k}\end{array}\right],\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\alpha }_{ij}=(\frac{X_{ik}{X}_{mj}}{X_{mk}}-X_{ij}),b={[\frac{X_{m1}}{X_{mk}},\frac{X_{m2}}{X_{mk}},\dots ,\frac{X_{mn}}{X_{mk}},1]}^{\mathrm{T}}\end{array}$$

至此得到矩阵形式的全部系数。

线性规划的对偶形式

这里补充一下关于线性规划问题转换为其对偶问题的数学推导⁴。已知有如下线性规划 $$\begin{array}{rl}PP:max& Z=C\cdot XX\\ s.t.& A\cdot XX\le b\\ & XX\ge 0\end{array}$$ 的对偶形式为 $$\begin{array}{rl}DD:min& W={YY}^{\mathrm{T}}\cdot b\\ s.t.& A^{\mathrm{T}}\cdot YY\ge C^{\mathrm{T}}\\ & YY\ge 0\end{array}$$ $DD$问题的自变量为$YY=[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{n+1}{]}^{\mathrm{T}}$，目标函数展开后为 $$W={YY}^{\mathrm{T}}\cdot b=\sum _{j=1}^n\frac{X_{mj}}{X_{mk}}{\lambda }_j+{\lambda }_{n+1}.$$ 主要约束条件展开后为 $$\begin{array}{rl}& \sum _{j=1}^n{\lambda }_j{Y}_{rj}\ge Y_{rk}\\ & \sum _{j=1}^n{\lambda }_j(\frac{X_{ik}{X}_{mj}}{X_{mk}}-X_{ij})\ge 0\\ & r=1,2,\dots ,q;i=1,2,\dots ,m\end{array}$$ 将目标函数代入第二个不等式，得到 $$\sum _{j=1}^n{\lambda }_j{X}_{ij}\le X_{ik}({\lambda }_{n+1}+W).$$ 以上对偶规划模型的多项式形式如下所示： $$\begin{array}{rl}min& W=\sum _{j=1}^n\frac{X_{mj}}{X_{mk}}{\lambda }_j+{\lambda }_{n+1}\\ s.t.& \sum _{j=1}^n{\lambda }_j{Y}_{rj}\ge Y_{rk}({\lambda }_j\ge 0)\\ & \sum _{j=1}^n{\lambda }_j{X}_{ij}\le X_{ik}({\lambda }_{n+1}+W)\\ & i=1,2,\dots ,m;r=1,2,\dots ,q;j=1,2,\dots ,n\end{array}$$