前言

本文源于笔者的《系统工程》课程的小组作业，笔者尝试运用DEMATEL-ISM方法来进行分析，建模求解，但在网络上并没有找到相应的，特别是集合DEMATEL-ISM方法的代码。因此自己码了DEMATEL-ISM模型的Python代码，并作为第一个博客发布~

参考文献中，笔者主要参考了李广利等¹的研究，本文也将依此论文进行方法解读和代码复现。

网上查找资料的过程中，笔者发现了一个MCDA方法的python代码库²，里面有很多多准则决策分析模型的相关代码，其中就有DEMATEL的代码，笔者做了一定的参考。

DEMATEL-ISM分析方法

方法简介

DEMATEL（Decision Making Trial and Evaluation Laboratory），全称为“决策试验和评价实验法”，是一种运用图论与矩阵工具进行系统要素分析的方法，通过分析系统中各要素之间的逻辑关系与直接影响关系，可以判断要素之间关系的有无及其强弱评价。ISM（Interpretative Structural Modelling）法全称为“解释结构模型”，其特点是把复杂的系统分解为若干子系统（要素），通过代数运算将系统构造成一个多级递阶的结构模型。

DEMATEL 模型可利用矩阵运算求出因素间的因果关系和影响强度，通过可视化因素间的因果关
系，得以揭示复杂问题中的关键影响因素及影响程度；但该方法无法有效识别系统中因素的层级结
构。ISM 法则通过分析构成系统的各子系统( 因素或要素) 之间的直接二元相关关系，基于布尔代数运算等，构造多级递阶有向拓扑图，但无法确定要素对系统的影响程度。

将两种方法结合，可以识别系统中关键要素及其影响程度，并构建要素的层级结构。DEMATEL-ISM方法的过程如下：

步骤

明确系统要素

明确分析系统所构成的要素，将构成系统的要素标记为
x
1
x_1
x1​,
x
2
x_2
x2​,
x
3
x_3
x3​, $\dots$,
x
n
x_n
xn​。

确定直接影响矩阵

采用专家打分法，比较
x
i
x_i
xi​对
x
j
x_j
xj​的影响，由于因素与自身比较为没有影响，直接影响矩阵的对角线值为0。通过比较得到直接影响矩阵$A$。

A
=
[
x
12
⋯
x
1
n
x
21
⋯
x
2
n
⋮
⋮
⋱
⋮
x
m
1
x
m
2
…
]
\begin{align} A= \begin{bmatrix} 0&x_{12}&\cdots&x_{1n}\\ x_{21}&0&\cdots&x_{2n}\\ \vdots &\vdots&\ddots &\vdots \\ x_{m1}&x_{m2}&\dots &0 \end{bmatrix} \end{align}
A=​0x21​⋮xm1​​x12​0⋮xm2​​⋯⋯⋱…​x1n​x2n​⋮0​​​​
式中因素
x
i
j
(
i
=
1
,
2
,
…
,
m
;
j
=
1
,
2
,
…
,
n
;
i
≠
j
)
x_{ij}(i=1,2, \ldots, m;j=1,2, \ldots,n;i\neq j)
xij​(i=1,2,…,m;j=1,2,…,n;i=j)表示因素
x
i
x_i
xi​对
x
j
x_j
xj​的直接影响。$i=j$时，
x
i
j
=
x_{ij}=0
xij​=0。

规范影响矩阵

归一化原始关系矩阵得到规范影响矩阵。归一化方法有很多种，这里可以采用行最大值法进行归一化，即将矩阵$A$的每一行求和，在这些值中取最大值，将矩阵$A$中元素除以最大值，得到规范直接影响矩阵$B$。

B
=
x
i
j
m
a
x
(
∑
j
=
1
n
x
i
j
)
\begin{align} B=\frac{x_{ij}}{max(\sum\limits_{j=1}^n x_{ij})} \end{align}
B=max(j=1∑n​xij​)xij​​​​

计算综合影响矩阵

综合系统矩阵体现系统中各个元素间的影响的综合效应，规范影响矩阵不断自乘后矩阵的所有值会趋近于0，即
lim
⁡
k
→
∞
B
k
=
\lim\limits_{k \to \infty} B^k=0
k→∞lim​Bk=0。计算综合影响矩阵时，得到

T
=
(
B
+
B
2
+
⋯
+
B
k
)
=
∑
k
=
1
∞
B
k
=
B
(
I
−
B
)
−
1
\begin{align} T=(B+B^2+\dots +B^k)=\sum\limits_{k=1}^\infty B^k=B(I-B)^{-1} \end{align}
T=(B+B2+⋯+Bk)=k=1∑∞​Bk=B(I−B)−1​​
式中为$I$单位矩阵。

计算各个要素的影响度、被影响度、中心度和原因度

影响度指矩阵$T$中各行值和，表示各行要素对其他所有要素的综合影响值，记作
D
i
D_i
Di​，有：

D
i
=
∑
j
=
1
n
x
i
j
,
(
i
=
1
,
2
,
…
,
n
)
\begin{align} D_i=\sum\limits_{j=1}^n x_{ij},(i=1,2,\dots,n) \end{align}
Di​=j=1∑n​xij​,(i=1,2,…,n)​​
被影响度指矩阵$T$中各列值和，表示各列要素对其他所有要素的综合影响值，记作
C
i
C_i
Ci​，有：

C
i
=
∑
j
=
1
n
x
j
i
,
(
i
=
1
,
2
,
…
,
n
)
\begin{align} C_i=\sum\limits_{j=1}^n x_{ji},(i=1,2,\dots,n) \end{align}
Ci​=j=1∑n​xji​,(i=1,2,…,n)​​
中心度表示因素在评价体系中的位置及其所起作用的大小，某要素的中心度为其影响度于被影响度之和，记作
M
i
M_i
Mi​，有：

M
i
=
D
i
+
C
i
\begin{align} M_i=D_i+C_i \end{align}
Mi​=Di​+Ci​​​
原因度由某要素的影响度和被影响度相减得到，记作，有：

R
i
=
D
i
−
C
i
\begin{align} R_i=D_i-C_i \end{align}
Ri​=Di​−Ci​​​
当原因度大于0时，表示该要素对其他要素的影响程度大，称其为原因要素，当原因度小于0时，该要素为结果要素。

绘制因果图

将中心度
M
i
M_i
Mi​为横坐标，将原因度
R
i
R_i
Ri​作为纵坐标，绘制因果关系图。该图可以直观体现因果关系。

确定整体影响矩阵

整体影响矩阵由综合影响矩阵和单位矩阵相加得到。

确定可达矩阵

确定可达矩阵时需要引入一个阈值$\lambda$以剔除因素之间影响程度较小的关系，从而明确层级结构的划分。引入阈值$\lambda$对整体影响矩阵$E$进行处理，可以得到可达矩阵$F$。则有：

f
i
j
=
{
1
e
i
j
≥
λ
(
i
,
j
=
1
,
2
,
…
,
n
)
e
i
j
<
λ
(
i
,
j
=
1
,
2
,
…
,
n
)
\begin{equation} f_{ij}= \begin{cases} 1 & & {e_{ij}\geq\lambda(i,j=1,2,\dots,n)}\\ 0 & & {e_{ij}<\lambda(i,j=1,2,\dots,n)} \end{cases} \end{equation}
fij​={10​​eij​≥λ(i,j=1,2,…,n)eij​<λ(i,j=1,2,…,n)​​​

划分层级

由可达矩阵$F$的第$i$行上值为1的列对应的因素求得可达集
R
(
x
1
)
=
{
x
i
∣
F
i
j
=
1
}
R(x_1)=\{x_i|F_{ij} =1\}
R(x1​)={xi​∣Fij​=1}，表示从因素
x
i
x_i
xi​出发可以达到的全部因素的合集，由可达矩阵$F$的第$i$列上值为1的行对应的因素求得先行集
S
(
x
1
)
=
{
x
i
∣
F
j
i
=
1
}
S(x_1)=\{x_i|F_{ji} =1\}
S(x1​)={xi​∣Fji​=1},表示可以达到因素的全部因素的合集。如果
R
(
x
1
)
R(x_1)
R(x1​)和
S
(
x
1
)
S(x_1)
S(x1​)满足
R
(
x
1
)
∩
S
(
x
1
)
=
R
(
x
1
)
R(x_1)\cap S(x_1)=R(x_1)
R(x1​)∩S(x1​)=R(x1​)，则表示中对应的元素均能在中找到前因，将该元素称为高层级的元素。然后从可达矩阵中去除对应的行和列，再从矩阵中抽取最高级的因素，不断重复该过程，直到所有的行和列均被去除。

绘制因素之间的递阶层次结构

根据去除因素的顺序，绘制系统要素间多级递阶有向拓扑图。

实例与代码

DEMATEL-ISM计算

这里笔者直接给出一组数据，10项因素，直接影响矩阵$A$如下。

下面是DEMETAL-ISM计算的Python代码：

#导入所需库
import numpy as np
import pandas as pd
#直接影响矩阵A
A = np.array([[0, 2, 0, 2, 3, 2, 3, 2, 3, 2],
       [0, 0, 0, 1, 3, 2, 0, 1, 1, 2],
       [0, 1, 0, 0, 2, 1, 0, 0, 0, 1],
       [0, 2, 0, 0, 3, 2, 0, 0, 0, 1],
       [1, 3, 0, 2, 0, 1, 1, 0, 1, 1],
       [0, 2, 0, 2, 1, 0, 1, 1, 0, 2],
       [2, 2, 0, 0, 2, 1, 0, 1, 2, 1],
       [0, 2, 0, 0, 1, 1, 1, 0, 1, 2],
       [3, 2, 0, 1, 2, 2, 2, 2, 0, 1],
       [0, 2, 0, 0, 1, 2, 1, 1, 2, 0]])
#行和最大值归一化，得到规范影响矩阵B
row_sum = np.sum(A, axis = 1)
max_sum = np.max(row_sum)
B = A/max_sum
#综合影响矩阵T
T = np.matmul (B, np.linalg.inv(np.identity(A.shape[0]) - B))
#计算影响度D，被影响度C，中心度M，原因度R
D = np.sum(T, axis = 1)
C = np.sum(T, axis = 0)
M = D + C
R = D - C
#判断输出结果要素和原因要素
causal_factors = ''
result_factors = ''
for i in range(len(R)):
    if R[i] > 0:
            causal_factors = causal_factors + 'x' + str(i + 1) + ' '
    elif R[i] < 0:
            result_factors = result_factors + 'x' + str(i + 1) + ' '
print("原因要素:" + causal_factors)
print("结果要素:" + result_factors)
#使用相应的lambda值计算可达矩阵F，并输出节点度排序
E = T + np.identity(A.shape[0])
lbd = 0.20      # 相应的lambda值
F = E
for i in range(A.shape[0]):
    for j in range(A.shape[1]):
        if E[i, j] > lbd:
            F[i,j] = 1
        elif E[i, j] < lbd:
            F[i, j] = 0
node_degree = np.sum(F, axis = 0) + np.sum(F, axis = 1)
print(sorted(node_degree, reverse=True))

上面代码计算出的一些矩阵和值比较重要，如综合影响矩阵$T$，被影响度$C$、影响度$D$、原因度$M$、中心度$R$等，数据在下面的画图也有用的，笔者在此没有输出，在做研究的过程中可以用pandas库将相应的值输出到excel表中。
$\lambda$的取值也很关键，具体取何值下文有说明，这里取0.20，得到最后的可达矩阵如下图所示。

原因-结果图

接下来运用上文计算的结果，绘制原因-结果图，这张图可以看到各个因子之间互相影响的大小，以及对研究内容影响的大小，如下图所示。

第1类为强原因因子集（第Ⅰ区），这类因子对研究对象的形成具有非常显著的影响，且对其他结果型因子有较大的影响；第2类为弱原因因子集（第Ⅱ区），这类因子对研究对象的形成也具有重要影响，对其他结果型因子也有一定的影响；第3类为弱结果因子集（第Ⅲ区），这类因子是其他原因型因子综合作用的结果，对研究对象的形成具有一定的影响；第4类为强结果因子集（第Ⅳ区），这类因子也是其他原因型因子综合作用的结果，但是对研究对象的形成具有非常重要的影响。DEMATEL计算所得因子中心度越大，其影响程度也就越大，因此，需要重点关注第Ⅰ区的强原因因子集及第Ⅳ区的强结果因子集。

下面是原因-结果图的代码，为了图形美观，笔者代码水平也很有限，所以做了大量的调整，大家可以看着自己来。

import matplotlib as mpl
import matplotlib.pyplot as plt
#设置图形格式
config = {
            "font.family": 'serif',
            "font.size": 14,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
            'axes.unicode_minus': False
          }
mpl.rcParams.update(config)
#x轴为上文计算出的中心度
x = [2.9456079466124003,
 3.243636237536544,
 0.5565493379272154,
 1.9429001271455832,
 3.205862930468748,
 2.627384244237508,
 2.4431692291267364,
 1.891879956077375,
 3.017054245988732,
 2.609503379763571]
y = [1.5128414543332274,
 -0.9939990519336557,
 0.5565493379272154,
 -0.1661886145842798,
 -0.8173828786366384,
 -0.6568575842419181,
 0.3753242725311403,
 -0.015046149922503416,
 0.6625117023194962,
 -0.4577524877920831]
#因子名
factors_name = [r'$x_1$', r'$x_2$', r'$x_3$', r'$x_4$', r'$x_5$',
                r'$x_6$', r'$x_7$', r'$x_8$', r'$x_9$', r'$x_{10}$']
#画散点图，并增加相应名称、线段和调整大小和位置
plt.scatter(x, y, s=3, c='k')
plt.xlabel('中心度')
plt.ylabel('原因度')
for i in range(len(x)):
    if i == 0:
        plt.text(x[i]+0.025, y[i]-0.085, factors_name[i], fontsize=17)
    elif i == 3:
        plt.text(x[i]+0.025, y[i]-0.1, factors_name[i], fontsize=17)
    elif i == 6:
        plt.text(x[i]-0.18, y[i]-0.1, factors_name[i], fontsize=17)
    elif i == 7:
        plt.text(x[i]-0.15, y[i]-0.12, factors_name[i], fontsize=17)
    else:
        plt.text(x[i]+0.025, y[i]+0.025, factors_name[i], fontsize=17)
plt.vlines(sum(x)/len(x), -1.55, 1.55, colors='k', linestyles='dashed')
plt.hlines(0, sum(x)/len(x)-2, sum(x)/len(x)+2, colors='k', linestyles='dashed')
plt.xlim(sum(x)/len(x)-2, sum(x)/len(x)+2)
plt.ylim(-1.55, 1.55)
plt.text(sum(x)/len(x)+2-0.2, 1.4-0.1, 'Ⅰ')
plt.text(sum(x)/len(x)-2+0.1, 1.4-0.1, 'Ⅱ')
plt.text(sum(x)/len(x)-2+0.1, -1.55+0.1, 'Ⅲ')
plt.text(sum(x)/len(x)+2-0.2, -1.55+0.1, 'Ⅳ')
plt.show()

节点度图

可达矩阵中每个因素所在行与所在列的和称为该因素的节点度，将各阈值对应的节点度由大到小排列，可得到不同阈值下疫情下大学生焦虑情绪形成因子的节点度衰减散点图。
为了得到可达矩阵，需要引入阈值剔除影响程度较小的关系，$\lambda$通常通过经验取值比较。这里分别取0.18、0.20、0.24和0.27

相应代码如下：

import matplotlib as mpl
import matplotlib.pyplot as plt
#设置图片格式
config = {
            "font.family": 'serif',
            "font.size": 14,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
            'axes.unicode_minus': False
         }
mpl.rcParams.update(config)
#x轴为递增的整数，y轴分别为不同lambda值时候的节点度，节点度在上文代码中可以算出
x = list(range(1, 11))
y_1 = [13.0, 12.0, 12.0, 11.0, 10.0, 10.0, 9.0, 7.0, 6.0, 2.0]   # lbd = 0.15
y_2 = [12.0, 11.0, 11.0, 8.0, 7.0, 7.0, 6.0, 5.0, 5.0, 2.0]      # lbd = 0.18
y_3 = [11.0, 8.0, 8.0, 7.0, 5.0, 4.0, 4.0, 4.0, 3.0, 2.0]        # lbd = 0.20
y_4 = [8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 3.0, 2.0, 2.0, 2.0]         # lbd = 0.24
#画图
plt.plot(x, y_1, marker='o', linestyle='-', label='λ=0.15', linewidth=1, color='black')
plt.plot(x, y_2, marker='x', linestyle='-', label='λ=0.18', linewidth=1, color='black')
plt.plot(x, y_3, marker='*', linestyle='-', label='λ=0.20', linewidth=1, color='black')
plt.plot(x, y_4, marker='D', linestyle='-', label='λ=0.24', linewidth=1, color='black')
plt.legend()
plt.xlabel("序号")
plt.ylabel("节点度")
plt.show()

ISM图

通过最后的步骤可以做出系统要素间多级递阶有向拓扑图，这里就不展示了，根据步骤即刻得到之间的关系并作图。

后记

DEMATEL-ISM模型，更多只是矩阵间的相互计算，想要得出结果还是比较容易的。
它可以用来分析因子间的相互关系，做出阶梯式的关系模型图，但是也有主观性强等的缺点。
因为笔者水平有限，所以分析和代码可能不尽人意，多多包含~