数据包络分析--CCR模型

支持博主有三宝，点赞、留言、关注公众号

个人公众号：数模派

写这篇文章的目的是分享一下我在学习CCR模型时遇到的一些疑惑以及查阅资料后得到的相应解答。

数据包络分析可以用来解决多指标问题。CCR模型是最早被提出来的数据包络分析方法。该模型的本质是线性规划。模型名称中的三个字母分别代表三个作者的名字的首字母。

1. 多指标问题

什么是多指标问题？《数学建模算法与应用》中给出了这样一个例子：某市教委需要对6所中学进行评价，其相应的指标如表格所示。表1中的生均投入和非低收入家庭百分比是输入指标，生均写作得分和生均科技得分是输出指标。请根据这些指标来评价学校。

该问题的实质是，根据输入、输出指标，给出一个最终的综合评价指标来给学校排序。看到这个例子，小编联想到了2020年华数杯C题(脱贫帮扶绩效评价)。这个题是可以用CCR模型来做的，但是小编当时不知道这个模型。正是因为这个原因，小编才决定深入学习CCR模型。

2. CCR模型

现有nnn个决策单元DMU，每个DMU都有mmm种输入和rrr种输出。Xj=(xj1,…,xji,…,xjm)TX_{j}=\left(x_{j 1}, \ldots, x_{j i}, \ldots, x_{j m}\right)^{T}Xj=(xj1,…,xji,…,xjm)T为第jjj个DMU的输入向量，其中，xjix_{ji}xji表示第jjj个DMU的第iii种输入。Yj=(yj1,…,yjs,…,yjr)TY_{j}=\left(y_{j 1}, \ldots, y_{js}, \ldots, y_{j r}\right)^{T}Yj=(yj1,…,yjs,…,yjr)T为第j个DMU的输出向量，其中，xjsx_{js}xjs表示第jjj个DMU的第sss种输出。为输入向量的权重，为输出向量的权重。CCR模型的数学模型可以表示为：
max⁡uTYj0vTXj0s.t. {uTYjvTXj≤1,j∈[1,n]u>0,v>0(1)\begin{array}{l} \max \quad \frac{u^{T} Y_{j_{0}}}{v^{T} X_{j_{0}}} \\ \text { s.t. }\left\{\begin{array}{c} \frac{u^{T} Y_{j}}{v^{T} X_{j}} \leq 1, j \in[1, n] \\ u>0, v>0 \end{array}\right. \end{array}\tag{1}maxvTXj0uTYj0 s.t. {vTXjuTYj≤1,j∈[1,n]u>0,v>0(1)

CCR模型是对每一个DMU进行评价，而每一个DMU的评分就是目标函数的值。这里的j0j_0j0表示nnn个决策单元DMU中的任意一个。以比率式作为评价指标，更符合实际意义。我们常说“投入和产出不成比例”就是这个意思。由于现实中的任何一项技术都不能使得输入全部转化为输出，所以uTYjvTXj≤1\frac{u^{T} Y_{j}}{v^{T} X_{j}} \leq 1vTXjuTYj≤1。

从第一个约束条件可以看出，在计算任意一个DMU的评分时，每一个DMU的评分都要求小于等于1。换句话说，第一个约束条件的个数实际上是nnn个。分式规划(分式作为目标函数)的缺点是它的解释不唯一。这是因为如果vvv、uuu是解，那么kvkvkv、kukuku也是解。这给求解增加了难度。为了求解方便，通常令目标函数的分母vTXj0=1v^TX_{j0}=1vTXj0=1或者通过Charnes-Cooper变化，将分式规划变为线性规划。我们采用后者。

令ω=tv,μ=tu,t=1vTXj0\omega=t v, \quad \mu=t u, \quad t=\frac{1}{v^{T} X_{j_{0}}}ω=tv,μ=tu,t=vTXj01，公式(1)可以变为：
max⁡μTYj0s.t.{ωTXj−μTYj≥0,j∈[1,n]ωTXj0=1ω>0,μ>0(2)\begin{aligned} &\max \quad \mu^{T} Y_{j_0}\\ &\text {s.t.}\left\{\begin{array}{c} \omega^{T} X_{j}-\mu^{T} Y_{j} \geq 0, j \in[1, n] \\ \omega^{T} X_{j_{0}}=1 \\ \omega>0, \mu>0 \end{array}\right. \end{aligned}\tag{2}maxμTYj0s.t.⎩⎨⎧ωTXj−μTYj≥0,j∈[1,n]ωTXj0=1ω>0,μ>0(2)

对于公式(2)，很多资料只给出了这个结果，本文将详细地说明其转化过程。首先是目标函数的转化过程。针对μ=tu\mu=tuμ=tu，我们先将等式两边同时做转置运算，再将等式两边同时除以uTu^TuT，得到t=μTuTt=\frac{\mu^T}{u^T}t=uTμT。将t=μTuTt=\frac{\mu^T}{u^T}t=uTμT、t=1vTXj0t=\frac{1}{v^{T} X_{j_{0}}}t=vTXj01代入目标函数，化解后就可以把分数规划转为线性规划。具体过程如下：
μ=tu⇒μT=tuT⇒t=μTuTmax⁡uTYj0vTXj0⇒max⁡tuTYj0⇒max⁡uTμTuTYj0⇒max⁡μTYj0\begin{array}{l} \mu=t u \Rightarrow \mu^{T}=t u^{T} \Rightarrow t=\frac{\mu^{T}}{u^{T}} \\ \max \quad \frac{u^{T} Y_{j_{0}}}{v^{T} X_{j_{0}}} \Rightarrow \max \quad t u^{T} Y_{j_{0}} \Rightarrow \max \quad u^{T} \frac{\mu^{T}}{u^{T}} Y_{j_{0}} \Rightarrow \max \quad \mu^{T} Y_{j_{0}} \end{array}μ=tu⇒μT=tuT⇒t=uTμTmaxvTXj0uTYj0⇒maxtuTYj0⇒maxuTuTμTYj0⇒maxμTYj0
接着是约束条件的转化过程。对于公式(2)中的第一个约束条件ωTXj0−μTYj0>0\omega^TX_{j_0}-\mu^TY_{j_0}>0ωTXj0−μTYj0>0的转化过程如下：
ω=tv,μ=tu⇒ωT=tvT,μT=tuTuTYjvTXj≤1⇒vTXj−uTYj≥0⇒ωTtXj−μTtYj≥0⇒ωTXj−μTYj≥0\begin{array}{l} \omega=t v, \mu=t u \Rightarrow \omega^{T}=t v^{T}, \mu^{T}=t u^{T} \\ \frac{u^{T} Y_{j}}{v^{T} X_{j}} \leq 1 \\ \Rightarrow v^{T} X_{j}-u^{T} Y_{j} \geq 0 \\ \Rightarrow \frac{\omega^{T}}{t} X_{j}-\frac{\mu^{T}}{t} Y_{j} \geq 0 \\ \Rightarrow \omega^{T} X_{j}-\mu^{T} Y_{j} \geq 0 \end{array}ω=tv,μ=tu⇒ωT=tvT,μT=tuTvTXjuTYj≤1⇒vTXj−uTYj≥0⇒tωTXj−tμTYj≥0⇒ωTXj−μTYj≥0

对于公式(2)中的第二个约束条件ωTXj0=1\omega^TX_{j_0}=1ωTXj0=1的转化过程如下：
t=1vTXj0⇒tvTXj0=1⇒ωTXj0=1\begin{array}{l} t=\frac{1}{v^{T} X_{j_{0}}} \\ \Rightarrow t v^{T} X_{j_{0}}=1 \\ \Rightarrow \omega^{T} X_{j_{0}}=1 \end{array}t=vTXj01⇒tvTXj0=1⇒ωTXj0=1
对于任何一个线性规划模型，它都存在一个等价的对偶模型。因此，我们可以把公式(2)写成它的对偶形式，如公式(3)所示。
max⁡θs.t. {∑j=1nλjXj≤θXj0∑j=1nλjYj≤Yj0λj≥0,j∈[1,n](3)\begin{aligned} &\max \quad \theta\\ &\text { s.t. }\left\{\begin{array}{l} \sum_{j=1}^{n} \lambda_{j} X_{j} \leq \theta X_{j_{0}} \\ \sum_{j=1}^{n} \lambda_{j} Y_{j} \leq Y_{j_{0}} \\ \lambda_{j} \geq 0, j \in[1, n] \end{array}\right. \end{aligned}\tag{3}maxθ s.t. ⎩⎨⎧∑j=1nλjXj≤θXj0∑j=1nλjYj≤Yj0λj≥0,j∈[1,n](3)
我们对公式(3)做一下说明。首先θ\thetaθ是一个数值，从别人的代码中我了解到，它就是我们对第个DMU的评分。如何从公式中得出这一结论，目前还没有想到。λj\lambda_jλj是未知的权重，它需要通过线性规划来求得。换句话说，我们需要找到一组λj\lambda_jλj，使得θ\thetaθ最小。由于许多资料只给出了公式(3)，没有给出其推导过程，使得大家很困惑，所以本文在这里给出推导过程，希望可以解除大家的困惑。关于线性规划的原始模型和对偶模型的一般式，如公式(4)所示，左侧为原始模型，右侧为对偶模型。
max⁡cTx′min⁡bTy′s.t.{Ax′≥bx′≥0s.t. {ATy′≤cy′≥0(4)\begin{array}{ll} \max\quad {c^{T}} x^{\prime} & \min\quad {b^{T}} y^{\prime} \\ \text {s.t.}\left\{\begin{array}{c} A x^{\prime} \geq b \\ x^{\prime} \geq 0 \end{array}\right. & \text { s.t. }\left\{\begin{array}{c} A^{T} y^{\prime} \leq c \\ y^{\prime} \geq 0 \end{array}\right. \end{array}\tag{4}maxcTx′s.t.{Ax′≥bx′≥0minbTy′ s.t. {ATy′≤cy′≥0(4)

对应公式(3)的目标函数μTYj0\mu^T Y_{j_0}μTYj0，未知量是μT\mu^TμT，而Yj0Y_{j_0}Yj0是已知量。我们对它做如下变化：
max⁡μTYj0⇒min⁡−μTYj0(5)\max \quad \mu^{T} Y_{j_{0}} \Rightarrow \min \quad-\mu^{T} Y_{j_{0}}\tag{5}maxμTYj0⇒min−μTYj0(5)
为了把公式(3)写成公式(4)的形式，我们认为：
c=[0,…,0⏟m,−yj01,…,−yj0]Tx′=[ω1,…,ωm,μ1,…μr]TA=[x11⋯x1m−y11⋯−y1r⋮⋱⋮⋮⋱⋮xn1⋯xnm−yn1⋯−ymr−xj01⋯−xj0m0⋯0]b=[0,…,0⏟n,−1]T(6)\begin{array}{l} c=[\underbrace{0, \ldots, 0}_{m},-y_{j_{0} 1}, \ldots,-y_{j_{0}}]^{T} \\ x^{\prime}=\left[\omega_{1}, \ldots, \omega_{m}, \mu_{1}, \ldots \mu_{r}\right]^{T} \\ A=\left[\begin{array}{cccccc} x_{11} & \cdots & x_{1 m} & -y_{11} & \cdots & -y_{1 r} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ x_{n 1} & \cdots & x_{n m} & -y_{n 1} & \cdots & -y_{m r} \\ -x_{j_{0} 1} & \cdots & -x_{j_{0} m} & 0 & \cdots & 0 \end{array}\right] \\ b =[\underbrace{0, \ldots, 0}_{n},-1]^{T} \end{array}\tag{6}c=[m0,…,0,−yj01,…,−yj0]Tx′=[ω1,…,ωm,μ1,…μr]TA=⎣⎢⎢⎢⎡x11⋮xn1−xj01⋯⋱⋯⋯x1m⋮xnm−xj0m−y11⋮−yn10⋯⋱⋯⋯−y1r⋮−ymr0⎦⎥⎥⎥⎤b=[n0,…,0,−1]T(6)
这里需要注意的是，公式(3)中的第一个约束条件实际上表示nnn个约束条件，把它写成矩阵的形式，就可以得到AAA、bbb；第二个约束条件需要在等式两边同时乘以-1。现在我们就可以得到公式(3)的对偶模型了。
目标函数推导：
max⁡bTy′⇒max⁡[0,…,0⏟n,−1][λ1,…λn,θ]T⇒max⁡−θ⇒min⁡θ\max \quad b^{T} y^{\prime} \Rightarrow \max \quad[\underbrace{0, \ldots, 0}_{n},-1]\left[\lambda_{1}, \ldots \lambda_{n}, \theta\right]^{T} \Rightarrow \max \quad-\theta \Rightarrow \min \quad \thetamaxbTy′⇒max[n0,…,0,−1][λ1,…λn,θ]T⇒max−θ⇒minθ
约束条件推导：
ATy′≤c⇒[x11⋯xn1−xj01⋮⋱⋮⋮x1m⋯xnm−xj0m−y11⋯−yn10⋮⋱⋮⋮−y1r⋯−ynr0][λ1,…λn,θ]T≤[0,…,0⏟m,−yj01,…,−yj0]T⇒{λ1x11+λ2x21+…+λnxn1−θxj01≤0⋮λ1x1m+λ2x2m+…+λnxnm−θxj0m≤0λ1y11+λ2y21+…+λnyn1≥yj01λ1y1r+λ2y2r+…+λnym≥yj0\begin{aligned} &A^{T} y^{\prime} \leq c\\ &\Rightarrow\left[\begin{array}{cccc} x_{11} & \cdots & x_{n 1} & -x_{j_{0} 1} \\ \vdots & \ddots & \vdots & \vdots \\ x_{1 m} & \cdots & x_{n m} & -x_{j_{0} m} \\ -y_{11} & \cdots & -y_{n 1} & 0 \\ \vdots & \ddots & \vdots & \vdots \\ -y_{1 r} & \cdots & -y_{n r} & 0 \end{array}\right]\left[\lambda_{1}, \ldots \lambda_{n}, \theta\right]^{T} \leq\left[\begin{array}{c} \left.\underbrace{0, \ldots, 0}_{m},-y_{j_{0} 1}, \ldots,-y_{j_{0}}\right]^{T} \\ \end{array}\right.\\ &\Rightarrow\left\{\begin{array}{c} \lambda_{1} x_{11}+\lambda_{2} x_{21}+\ldots+\lambda_{n} x_{n 1}-\theta x_{j_{0} 1} \leq 0 \\ \vdots \\ \lambda_{1} x_{1 m}+\lambda_{2} x_{2 m}+\ldots+\lambda_{n} x_{n m}-\theta x_{j_{0} m} \leq 0 \\ \lambda_{1} y_{11}+\lambda_{2} y_{21}+\ldots+\lambda_{n} y_{n 1} \geq y_{j_{0} 1} \\ \lambda_{1} y_{1 r}+\lambda_{2} y_{2 r}+\ldots+\lambda_{n} y_{m} \geq y_{j_{0}} \end{array}\right. \end{aligned}ATy′≤c⇒⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡x11⋮x1m−y11⋮−y1r⋯⋱⋯⋯⋱⋯xn1⋮xnm−yn1⋮−ynr−xj01⋮−xj0m0⋮0⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤[λ1,…λn,θ]T≤[m0,…,0,−yj01,…,−yj0]T⇒⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧λ1x11+λ2x21+…+λnxn1−θxj01≤0⋮λ1x1m+λ2x2m+…+λnxnm−θxj0m≤0λ1y11+λ2y21+…+λnyn1≥yj01λ1y1r+λ2y2r+…+λnym≥yj0