shapley计算示例

目的

我们设置1个模拟电量租户根据提供的电量，应该获取收益的计算。

1、电价和电量是满足一定的关系的，主要是市场提供的电量增多，电价下降。具体查看公式calculate_y。

2、租户（电厂）提供电量给电网，电网根据电量还有电价，按照shapley的计算方法，计算出每个租户应得到的利益。

3、函数calculate_shapley输入X是租户供给的电量，输出是每个租户应得利益。

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
86
87
88
89
90
91
92
93
94
import itertools
import math
import bisect
import sys

#参考： https://github.com/susobhang70/shapley_value/blob/master/shapley.py

# 电价和电量的关系的斜率值
XIE_LV = -2
# 初始单价
INIT_PRICE = 500

# 定义函数
def calculate_y(x):
    return x * (XIE_LV*x + INIT_PRICE)

def calcuate_total_value(X=[1, 2, 3, 4, 5, 6, 7, 8]):
    # 定义x的 possible values
    # 找出所有可能的组合
    combinations = []
    for i in range(1, len(X) + 1):
        # if i != 8: continue
        combinations.extend(itertools.combinations(X, i))

    total = []
    # 根据每个组合计算y的值
    for combination in combinations:
        sum_x = sum(combination)
        y_values = calculate_y(sum_x)
        print(f"For combination {combination}, y values are {y_values}")
        total.append(y_values)
    return total

def power_set(List):
    PS = [list(j) for i in range(len(List)) for j in itertools.combinations(List, i+1)]
    return PS

def calculate_shapley(X=[1, 2, 3, 4, 5, 6, 7, 8]):
    """

    Args:
        X (): 代表每个用户提供的电量

    Returns:

    """
    # 用户数量
    X_length = len(X)
    # 计算所有可能的组合, zuhe: 1个组合 [[0], [1], [2], [3], [0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3], [0, 1, 2], [0, 1, 3], [0, 2, 3], [1, 2, 3], [0, 1, 2, 3]]
    tempList = list([i for i in range(X_length)])
    zuhe = power_set(tempList)
    # 计算每个组合的收益,[-228, -375, -627, -3000, -1323, -1767, -5148, -2142, -5775, -6667, -4002, -8643, -9727, -10582, -14362]
    total_prices = calcuate_total_value(X)

    assert len(zuhe) == len(total_prices), "组合的数量和组合的收益的数量一定相等"
    shapley_values = []
    for i in range(X_length):  #遍历每个用户的索引
        shapley = 0
        for j in zuhe: #遍历每个组合
            if i not in j:
                cmod = len(j)  #计算组合里面的用户数量
                Cui = j[:]   #Cui是j的复制 [1, 3], 目前i:0
                bisect.insort_left(Cui,i)  #把i 插入到cui左边,[0, 1, 3]
                l = zuhe.index(j)  #找出目前[1,3]时的索引, eg:8
                k = zuhe.index(Cui)   #找出插入之后[0, 1, 3]的索引， eg:10
                # temp = float(float(total_prices[k]) - float(total_prices[l])) *\
                #            float(math.factorial(cmod) * math.factorial(X_length - cmod - 1)) / float(math.factorial(X_length))
                # 表示增加1个用户后的收益差值, 具体查看shaply计算公式
                chazhi = float(float(total_prices[k]) - float(total_prices[l]))
                # 阶乘的计算,分子
                # jiecheng_jisuan = float(math.factorial(cmod) * math.factorial(X_length - cmod - 1))
                # 修订版本,
                jiecheng_jisuan = float(math.factorial(cmod-1) * math.factorial(X_length - cmod))
                # 除数阶乘，分母
                chushu_jiecheng = float(math.factorial(X_length))
                # 临时权重
                temp = chazhi * jiecheng_jisuan / chushu_jiecheng
                # 总得权重，对应求和公式
                shapley += temp
                # if i is 0:
                #     print j, Cui, cmod, n-cmod-1, characteristic_function[k], characteristic_function[l], math.factorial(cmod), math.factorial(n - cmod - 1), math.factorial(n)
        cmod = 0
        Cui = [i]
        k = zuhe.index(Cui)
        temp = float(total_prices[k]) * float(math.factorial(cmod) * math.factorial(X_length - cmod - 1)) / float(math.factorial(X_length))
        shapley += temp
        shapley_values.append(shapley)
    print(shapley_values)  # 每个人应该得到的利益


if __name__ == '__main__':
    # print(calcuate_total_value(X=[12,15,19,40]))
    calculate_shapley(X=[12,15,19,40,90])
    # calculate_shapley(X=[2,4])

上面代码没有考虑只有自身的情况，修改代码如下

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
86
87
88
89
90
91
92
93
94
95
96
97
98
import itertools
import math
import bisect

#参考： https://github.com/susobhang70/shapley_value/blob/master/shapley.py
#进行了bug修复

# 电价和电量的关系的斜率值
XIE_LV = -2
# 初始单价
INIT_PRICE = 500

# 定义函数
def calculate_y(x):
    return x * (XIE_LV*x + INIT_PRICE)

def calcuate_total_value(X=[1, 2, 3, 4, 5, 6, 7, 8]):
    # 定义x的 possible values
    # 找出所有可能的组合
    combinations = []
    for i in range(1, len(X) + 1):
        # if i != 8: continue
        combinations.extend(itertools.combinations(X, i))

    total = []
    # 根据每个组合计算y的值
    for combination in combinations:
        sum_x = sum(combination)
        y_values = calculate_y(sum_x)
        print(f"For combination {combination}, y values are {y_values}")
        total.append(y_values)
    return total

def power_set(List):
    PS = [list(j) for i in range(len(List)) for j in itertools.combinations(List, i+1)]
    return PS

def calculate_shapley(X=[1, 2, 3, 4, 5, 6, 7, 8]):
    """

    Args:
        X (): 代表每个用户提供的电量

    Returns:

    """
    # 用户数量
    X_length = len(X)
    # 计算所有可能的组合, zuhe: 1个组合 [[0], [1], [2], [3], [0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3], [0, 1, 2], [0, 1, 3], [0, 2, 3], [1, 2, 3], [0, 1, 2, 3]]
    tempList = list([i for i in range(X_length)])
    zuhe = power_set(tempList)
    # 计算每个组合的收益,[-228, -375, -627, -3000, -1323, -1767, -5148, -2142, -5775, -6667, -4002, -8643, -9727, -10582, -14362]
    total_prices = calcuate_total_value(X)

    assert len(zuhe) == len(total_prices), "组合的数量和组合的收益的数量一定相等"
    shapley_values = []
    for i in range(X_length):  #遍历每个用户的索引
        shapley = 0
        for j in zuhe: #遍历每个组合
            if i not in j:
                Cui = j[:]   #Cui是j的复制 [1, 3], 目前i:0
                bisect.insort_left(Cui,i)  #把i 插入到cui左边,[0, 1, 3]
                l = zuhe.index(j)  #找出目前[1,3]时的索引, eg:8
                total_l_price = float(total_prices[l])
            elif i in j and len(j) == 1:
                Cui = j[:]
                total_l_price = 0
            else:
                continue
            k = zuhe.index(Cui)   #找出插入之后[0, 1, 3]的索引， eg:10
            cmod = len(Cui)  #计算组合里面的用户数量
            # 表示增加1个用户后的收益差值, 具体查看shaply计算公式
            chazhi = float(float(total_prices[k]) - total_l_price)
            # 阶乘的计算,分子
            # jiecheng_jisuan = float(math.factorial(cmod) * math.factorial(X_length - cmod - 1))
            # 修订版本,
            jiecheng_jisuan = float(math.factorial(cmod-1) * math.factorial(X_length - cmod))
            # 除数阶乘，分母
            chushu_jiecheng = float(math.factorial(X_length))
            # 临时权重
            temp = chazhi * jiecheng_jisuan / chushu_jiecheng
            # 总得权重，对应求和公式
            shapley += temp
                # if i is 0:
                #     print j, Cui, cmod, n-cmod-1, characteristic_function[k], characteristic_function[l], math.factorial(cmod), math.factorial(n - cmod - 1), math.factorial(n)
        cmod = 0
        Cui = [i]
        k = zuhe.index(Cui)
        temp = float(total_prices[k]) * float(math.factorial(cmod) * math.factorial(X_length - cmod - 1)) / float(math.factorial(X_length))
        shapley += temp
        shapley_values.append(shapley)
    print(shapley_values)  # 每个人应该得到的利益


if __name__ == '__main__':
    # print(calcuate_total_value(X=[12,15,19,40]))
    calculate_shapley(X=[12,15,19,40,90])
    # calculate_shapley(X=[2,4])
shapley计算示例
https://johnson7788.github.io/2024/03/04/shapley%E8%AE%A1%E7%AE%97%E7%A4%BA%E4%BE%8B/
作者
Johnson
发布于
2024年3月4日
许可协议