学習目標
CCRモデルにおける効率性 \theta の定義を理解する．
CCRモデルを分数計画問題として定式化できる．
分数計画問題を線形計画問題に変換する理由と方法を理解する．
最適値 \theta^* と最適な重みから，DMUが効率的か非効率的かを判定できる．
参照集合の意味を理解し，求めることができる．
データと重み
CCRモデルは，Charnes, Cooper, and Rhodesが1978年に提案した，DEAの基本となるモデルである(Charnes ほか 1978年 ) ．複数の入力と出力をそれぞれ重み付き和にまとめ，その比によってDMUの相対的な効率性を評価する．
それぞれ m 個の入力と s 個の出力を持つ n 個のDMUを考える．\text{DMU}_j （j=1,2,\ldots,n ）の第 i 入力を x_{ij} ，第 r 出力を y_{rj} と表す．入力データと出力データを，それぞれ次の行列にまとめる．
\mathbf{X} =
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1n} \\
x_{21} & x_{22} & \cdots & x_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
x_{m1} & x_{m2} & \cdots & x_{mn}
\end{pmatrix}
\mathbf{Y} =
\begin{pmatrix}
y_{11} & y_{12} & \cdots & y_{1n} \\
y_{21} & y_{22} & \cdots & y_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
y_{s1} & y_{s2} & \cdots & y_{sn}
\end{pmatrix}
ここで，各行は入力または出力の項目に対応し，各列は1つのDMUに対応する．したがって，\mathbf{X} は m \times n 行列，\mathbf{Y} は s \times n 行列である．
入力と出力に対する非負の重みを，それぞれ次のベクトルで表す．
\mathbf{v} =
\begin{pmatrix}
v_1 \\ v_2 \\ \vdots \\ v_m
\end{pmatrix},\quad
\mathbf{u} =
\begin{pmatrix}
u_1 \\ u_2 \\ \vdots \\ u_s
\end{pmatrix}
評価対象のDMUを \text{DMU}_o とする．与えられた重みのもとで，その効率性 \theta を，入力の重み付き和に対する出力の重み付き和の比として定義する．すなわち，
\theta = \frac{\sum_{r=1}^{s} u_r y_{ro}}{\sum_{i=1}^{m} v_i x_{io}} = \frac{u_1 y_{1o} + u_2 y_{2o} + \cdots + u_s y_{so}}{v_1 x_{1o} + v_2 x_{2o} + \cdots + v_m x_{mo}}
である．ただし，分母である入力の重み付き和は正とする．このように重みを用いることで，複数の入力と出力から1つの効率性の指標を得ることができる．
例 5.1 (与えられた重みによる効率性の計算) 表 5.1 に示す，2入力2出力の3つのDMUを考える．
各DMUのデータを列ごとに並べると，入力行列 \mathbf{X} と出力行列 \mathbf{Y} は次のようになる．
\mathbf{X} =
\begin{pmatrix}
20 & 19 & 25 \\
151 & 131 & 160
\end{pmatrix}, \quad
\mathbf{Y} =
\begin{pmatrix}
100 & 150 & 160 \\
90 & 50 & 55
\end{pmatrix}
入力と出力に対する重みベクトルは，それぞれ
\mathbf{v} =
\begin{pmatrix}
v_1 \\ v_2
\end{pmatrix}, \quad
\mathbf{u} =
\begin{pmatrix}
u_1 \\ u_2
\end{pmatrix}
である．この重みを用いたときの各DMUの効率性を \theta_1 ，\theta_2 ，\theta_3 とすると，\text{DMU}_1 については，
\begin{align*}
\theta_1 &= \frac{\sum_{r=1}^{s} u_r y_{r1}}{\sum_{i=1}^{m} v_i x_{i1}} \\
&= \frac{u_1 y_{11} + u_2 y_{21}}{v_1 x_{11} + v_2 x_{21}} \\
&= \frac{100 u_1 + 90 u_2}{20 v_1 + 151 v_2}
\end{align*}
となる．同じ重みを用いると，\text{DMU}_2 と \text{DMU}_3 の効率性は，それぞれ次のように表される．
\begin{align*}
\text{DMU}_2: \quad & \theta_2 = \frac{\sum_{r=1}^{s} u_r y_{r2}}{\sum_{i=1}^{m} v_i x_{i2}} = \frac{150 u_1 + 50 u_2}{19 v_1 + 131 v_2} \\
\text{DMU}_3: \quad & \theta_3 = \frac{\sum_{r=1}^{s} u_r y_{r3}}{\sum_{i=1}^{m} v_i x_{i3}} = \frac{160 u_1 + 55 u_2}{25 v_1 + 160 v_2}
\end{align*}
例えば，重みを
\mathbf{u} =
\begin{pmatrix}
0.003 \\ 0.007
\end{pmatrix}, \quad
\mathbf{v} =
\begin{pmatrix}
0.005 \\ 0.006
\end{pmatrix}
とすると，\text{DMU}_1 の効率性は次のように計算できる．
\begin{align*}
\theta_1
&= \frac{100 \times 0.003 + 90 \times 0.007}
{20 \times 0.005 + 151 \times 0.006} \\
&= \frac{0.930}{1.006} \\
&\approx 0.9245
\end{align*}
この値は，ここで与えた重みに対する効率性であり，CCRモデルによって重みを最適化した結果ではない．
効率性の値は，重みの選び方によって変わる．CCRモデルでは，すべてのDMUに同じ重みを適用したとき，どのDMUの効率性も1を超えないという条件のもとで，評価対象の \text{DMU}_o の効率性を最大化する重み \mathbf{u} と \mathbf{v} を求める．つまり，他のDMUとの比較による制約を満たす範囲で，評価対象に最も有利な重みを選ぶのである．評価対象を別のDMUに変えるときは，そのDMUについて改めて重みを最適化する．
CCRモデル
評価対象の \text{DMU}_o の効率性を最大化する目的関数は，次のように表される．
\text{maximize} \quad \theta = \frac{\sum_{r=1}^{s} u_r y_{ro}}{\sum_{i=1}^{m} v_i x_{io}}
\tag{5.1}
ここでは，次の仮定を置く．
すべての入力データと出力データが正である．
入力と出力の重みはいずれも非負とし，入力の重みについては，少なくとも1つが正であるとする．
この仮定により，すべてのDMUについて，分母となる入力の重み付き和は正となる．
一方，同じ重みを用いたときの各DMUの効率性が1以下となるように，次の制約を課す．
\frac{\sum_{r=1}^{s} u_r y_{rj}}
{\sum_{i=1}^{m} v_i x_{ij}} \leq 1,
\quad j=1,\ldots,n
\tag{5.2}
以上をまとめると，上記の仮定のもとで，CCRモデルは次の分数計画問題として定式化できる．
\begin{aligned}
\text{maximize} \quad
& \theta = \frac{\sum_{r=1}^{s} u_r y_{ro}}
{\sum_{i=1}^{m} v_i x_{io}} \\
\text{subject to} \quad
& \frac{\sum_{r=1}^{s} u_r y_{rj}}
{\sum_{i=1}^{m} v_i x_{ij}} \leq 1,
&& j=1,\ldots,n, \\
& u_r\geq 0,
&& r=1,\ldots,s, \\
& v_i\geq 0,
&& i=1,\ldots,m.
\end{aligned}
この問題の最適値を \theta^* ，最適な重みを \mathbf{u}^* ，\mathbf{v}^* と表す．評価対象の \text{DMU}_o 自身も制約に含まれるため，\theta^*\leq 1 である．
なお，入力と出力の重みを同じ正の定数倍しても，各DMUの効率性は変わらない．実際，任意の \alpha>0 に対して，
\frac{\sum_{r=1}^{s}(\alpha u_r)y_{rj}}
{\sum_{i=1}^{m}(\alpha v_i)x_{ij}}
=
\frac{\sum_{r=1}^{s}u_r y_{rj}}
{\sum_{i=1}^{m}v_i x_{ij}}
が成り立つ．したがって，(\mathbf{u}^*,\mathbf{v}^*) が最適解であれば，(\alpha\mathbf{u}^*,\alpha\mathbf{v}^*) も同じ最適値を与える最適解となる．このように，重みの尺度には任意性があるため，この分数計画問題には無限に多くの最適解が存在する．
線形計画問題への変換
重みの尺度の任意性を利用して，評価対象の入力の重み付き和が1となるように重みを正規化すると，分数計画問題は線形計画問題に変換できる．
定理 5.1
\begin{aligned}
\text{maximize} \quad
& \theta=\sum_{r=1}^{s}\mu_r y_{ro} \\
\text{subject to} \quad
& \sum_{i=1}^{m}\nu_i x_{io}=1, \\
& \sum_{r=1}^{s}\mu_r y_{rj}
-\sum_{i=1}^{m}\nu_i x_{ij}\leq 0,
&& j=1,\ldots,n, \\
& \mu_r\geq 0,
&& r=1,\ldots,s, \\
& \nu_i\geq 0,
&& i=1,\ldots,m.
\end{aligned}
また，この線形計画問題の最適解を (\boldsymbol{\mu}^*,\boldsymbol{\nu}^*) とすると，\mathbf{u}=\boldsymbol{\mu}^* ，\mathbf{v}=\boldsymbol{\nu}^* は分数計画問題の最適解である．
証明 . (\mathbf{u},\mathbf{v}) に対して，
t=\frac{1}{\sum_{i=1}^{m}v_i x_{io}}>0,\qquad
\mu_r=t u_r,\qquad
\nu_i=t v_i
とおく．このとき，
\sum_{i=1}^{m}\nu_i x_{io}
=t\sum_{i=1}^{m}v_i x_{io}
=1
となり，目的関数は
\theta
=\frac{\sum_{r=1}^{s}u_r y_{ro}}{\sum_{i=1}^{m}v_i x_{io}} = t \sum_{r=1}^{s}u_r y_{ro} = \sum_{r=1}^{s}\mu_r y_{ro}
と表される．また，各DMUの入力の重み付き和は正であるため，式 5.2 は
\sum_{r=1}^{s}u_r y_{rj}
-\sum_{i=1}^{m}v_i x_{ij}\leq 0,
\quad j=1,\ldots,n
と書き換えられる．さらに，両辺に t>0 を掛けると，
\sum_{r=1}^{s}\mu_r y_{rj}
-\sum_{i=1}^{m}\nu_i x_{ij}\leq 0,
\quad j=1,\ldots,n
となる．したがって，(\boldsymbol{\mu},\boldsymbol{\nu}) は線形計画問題の実行可能解であり，その目的関数値は元の分数計画問題の目的関数値に等しい．
逆に，線形計画問題の任意の実行可能解 (\boldsymbol{\mu},\boldsymbol{\nu}) に対して，\mathbf{u}=\boldsymbol{\mu} ，\mathbf{v}=\boldsymbol{\nu} とおく．\sum_{i=1}^{m}\nu_i x_{io}=1 より入力の重みの少なくとも1つは正であるから，各DMUの入力の重み付き和は正となる．したがって，線形計画問題の制約の両辺をこの和で割ると 式 5.2 が得られ，(\mathbf{u},\mathbf{v}) は元の分数計画問題でも実行可能となる．また，目的関数の分母は1であるから，目的関数値も一致する．
以上より，両問題の最適値は一致する．さらに，線形計画問題の最適解は，分数計画問題の実行可能解として同じ目的関数値，すなわち最適値を与えるから，分数計画問題の最適解である．
以下では，定理 5.1 の線形計画問題の最適解，すなわち正規化された最適な重みを \mathbf{u}^* ，\mathbf{v}^* と表す．ただし，正規化によって最適な重みが一意に定まるとは限らない．
効率性の判定
最適値 \theta^* と最適な重みに基づいて，CCR効率性を次のように定義する．
定義 5.1
\theta^*=1 であり，かつ \mathbf{u}^*>\mathbf{0} ，\mathbf{v}^*>\mathbf{0} を満たす最適解が少なくとも1つ存在するとき，\text{DMU}_o はCCR効率的 （CCR-efficient）である．それ以外の場合，\text{DMU}_o はCCR非効率的 （CCR-inefficient）である．
ここで，\mathbf{u}^*>\mathbf{0} ，\mathbf{v}^*>\mathbf{0} は，入力と出力のすべての重みが正であることを意味する．したがって，\theta^*=1 であっても，すべての最適解において少なくとも1つの重みが0となる場合には，CCR非効率的と判定される．ただし，得られた1つの最適解に0の重みが含まれているだけでは，CCR非効率的とは断定できない．すべての重みが正となる別の最適解が存在する可能性があるためである．
参照集合
\theta^*<1 の場合を考える．このとき，少なくとも1つの j\neq o について，最適解における制約が等号で成り立つ．すなわち，
\sum_{r=1}^{s}u_r^* y_{rj}
-\sum_{i=1}^{m}v_i^* x_{ij}=0
となるDMUが存在する．
もしすべてのDMUについて制約が厳密な不等号で成り立つならば，入力の重みを固定したまま出力の重みを一斉にわずかに増やすことで，すべての制約を満たしながら目的関数値を増加させることができる．これは最適性に反する．また，\theta^*<1 であるため，評価対象の \text{DMU}_o 自身の制約は等号では成り立たない．したがって，等号で成り立つ制約に対応するDMUは，評価対象以外のDMUである．
定義 5.2 \theta^*<1 である \text{DMU}_o について，1つの最適解 (\mathbf{u}^*,\mathbf{v}^*) を固定し，制約が等号で成り立つDMUの添字の集合を
\mathcal{R}'_o
=
\left\{
j\in\{1,\ldots,n\}
\;\middle|\;
\sum_{r=1}^{s}u_r^*y_{rj}
-\sum_{i=1}^{m}v_i^*x_{ij}=0
\right\}
とする．さらに，\mathcal{R}'_o のうち，CCR効率的なDMUに対応する添字からなる部分集合 \mathcal{R}_o を，\text{DMU}_o の参照集合 （reference set）と呼ぶ．
参照集合に含まれるDMUは，評価対象の \text{DMU}_o の非効率性を示す比較基準となる． ただし，制約が等号で成り立つことだけでは，そのDMUがCCR効率的であるとは限らないため， 一般に \mathcal{R}'_o と \mathcal{R}_o は区別する必要がある． 一方，定義で用いた最適解のすべての重みが正であれば，次のように両者は一致する．
PythonによるCCRモデルの実装
以下に，PythonでCCRモデルを解く関数を示す．Gurobiをソルバーとして使用している．Gurobiのインストールと基本的な使い方は チャプター 3
import gurobipy as gpfrom gurobipy import GRBimport numpy as npdef ccr(inputs, outputs, dmu_index):= inputs.shape[1 ]= inputs.shape[0 ]= outputs.shape[0 ]= gp.Model("DEA_CCR" )= 0 # ログを表示しない # 変数の定義 = m.addVars(num_outputs, name= "u" )= m.addVars(num_inputs, name= "v" )# 目的関数の定義 * outputs[r, dmu_index] for r in range (num_outputs)),# 制約条件の定義 * inputs[i, dmu_index] for i in range (num_inputs)) == 1 ,"input_normalization" ,= m.addConstrs(* outputs[r, j] for r in range (num_outputs))- gp.quicksum(v[i] * inputs[i, j] for i in range (num_inputs))<= 0 for j in range (num_dmu)"theta_constraint" ,# 最適化の実行 if m.Status == GRB.OPTIMAL:= m.ObjVal= [u[r].X for r in range (num_outputs)]= [v[i].X for i in range (num_inputs)]# 制約が等号で成り立つDMU（数値誤差を考慮し，スラックが十分小さい制約を等号とみなす） = []for j in range (num_dmu):if abs (theta_constraint[j].Slack) < 1e-6 :return theta, output_weights, input_weights, bindingelse :return None ccr()関数は，最適値 \theta^* ，最適な重み \mathbf{u}^* ，\mathbf{v}^* と，得られた最適解において制約が等号で成り立つDMUの添字のリストを返す．\theta^*<1 のとき，このリストは \mathcal{R}'_o である．最適解が複数存在する場合，このリストは得られた最適解に依存する．
ccr()関数を用いて，すべてのDMUについてCCRモデルを解く関数を以下に示す．
def solve_all_ccr(inputs, outputs):= np.atleast_2d(inputs)= np.atleast_2d(outputs)= inputs.shape[1 ]= [chr (65 + i) for i in range (num_dmu)]= []= []= []= []for i in range (num_dmu):= ccr(inputs, outputs, i)for j in binding])for i in range (num_dmu):= dmu_names[i]= results[i]= u_list[i]= v_list[i]= binding_sets[i]= ", " .join(f" { w:.4f} " for w in u_val)= ", " .join(f" { w:.4f} " for w in v_val)= ", " .join(binding)print (f"DMU { name} " )print (f" theta : { theta:.4f} " )print (f" Output Weights u : [ { formatted_u} ]" )print (f" Input Weights v : [ { formatted_v} ]" )print (f" Binding DMUs : [ { formatted_binding} ]" )print ("-" * 50 )
1入力1出力のCCRモデル
1入力1出力の分数計画問題において，目的関数と制約は
\begin{align*}
\frac{u y_o}{v x_o} &= \frac{u}{v} \frac{y_o}{x_o}, \\
\frac{u y_j}{v x_j} &= \frac{u}{v} \frac{y_j}{x_j} \leq 1, \quad j=1,\ldots,n
\end{align*}
となる．制約より u/v \leq 1/\max_j (y_j/x_j) であり，この上限は評価対象によらない．目的関数は u/v が大きいほど大きいから，u/v = 1/\max_j (y_j/x_j) のとき最大となり，最適値は
\theta^* = \frac{y_o/x_o}{\max_j (y_j/x_j)}
となる．したがって，1入力1出力のCCRモデルの最適化問題を解くことは，すべてのDMUについて y_j/x_j を計算し，その最大値で割ることと同じである．
まとめ
すべてのDMUの効率性が1以下であるという制約を課すことで，\theta は [0, 1] に収まる．
分数計画問題の最適解は正の定数倍しても最適解であるため，解が一意に定まらない．分母を1に固定することで，線形計画問題に変換できる．
\theta^* = 1 であり，かつ \mathbf{u}^* > \mathbf{0} ，\mathbf{v}^* > \mathbf{0} を満たす最適解が存在するとき，\text{DMU}_o はCCR効率的である．\theta^* < 1 の場合と，\theta^* = 1 であってもすべての最適解に0の重みが含まれる場合は，CCR非効率的である．\theta^* < 1 のとき，最適解において制約が等号で成り立つDMUのうち，CCR効率的なDMUの集合が参照集合であり，\text{DMU}_o を非効率的にしているDMUを表す．
文献案内
DEAは，Charnes ほか (1978年 ) が提案したCCRモデルに始まる．
本書で扱わないモデルや応用例については，Cooper ほか (2007年 ) が体系的な教科書であり，Cooper ほか (2011年 ) には各分野の専門家による解説が収められている．
補足
定理 5.1 の後で述べたように，正規化しても最適な重みが一意に定まるとは限らない．これは \theta^*<1 の場合も同様であり，非効率的なDMUについても最適解が複数存在することがある．定義 5.2 では1つの最適解を固定して \mathcal{R}'_o を定めるから，最適解が複数存在する場合には，どの最適解を固定するかによって \mathcal{R}'_o が変わることがあり，その結果，参照集合 \mathcal{R}_o も変わることがある．
練習問題
練習 5.1 表 5.5 は，3つのDMUの入力と出力を示している．
入力と出力に対する重みが，それぞれ次のように与えられている．
\mathbf{v} =
\begin{pmatrix}
0.1 \\ 0.2
\end{pmatrix}, \quad
\mathbf{u} =
\begin{pmatrix}
0.05 \\ 0.1
\end{pmatrix}
この重みを用いて，以下の問いに答えよ．
各DMUについて，入力の重み付き和と出力の重み付き和を求めよ．
各DMUの効率性 \theta を求めよ．
与えられた重みのもとで，効率性が最も高いDMUはどれか．
Charnes, A, W W Cooper, と E Rhodes. 1978年. 「Measuring the efficiency of decision making units」 . Eur. J. Oper. Res. 2 (6): 429–44.
Cooper, William W, Lawrence M Seiford, と Kaoru Tone. 2007年. Data envelopment analysis: A comprehensive text with models, applications, references and DEA-Solver software . 2番目 版. Springer.
Cooper, William W, Lawrence M Seiford, と Joe Zhu, 編. 2011年. Handbook on data envelopment analysis . 2番目 版. Springer.
