用shapley法解决3个村庄合作修建电视接收塔的成本分摊博弈的matlab实现

linprog函数

首先，linprog函数主要用来求线型规划中的最小值问题（最大值的镜像问题，求最大值只需要加个“-”）。
算法结构及使用方法：针对约束条件为Ax=b或Ax≤b的问题。

x=linprog(f,A,b)
x=linprog(f,A,b,Aeq,beq)
x=linprog(f,A,b,Aeq,beq,lb,ub)
x=linprog(f,A,b,Aeq,beq,lb,ub,x0)
%linprog函数
%参数简介
%f：目标函数
%A：不等式约束条件矩阵
%b：对应不等式右侧的矩阵
%Aeq：等式约束条件矩阵
%beq：对应等式右侧的矩阵
%lb：x的下界
%ub：x的上界
%x0：设置初始点x0，这个选择项只是对medium-scale算法有效。默认的large-scale算法和简单的算法忽略任何初始点。

**

例1.求解线性规划

输入代码：

f=[-5 -4 -6];
    A=[1 -1 -1 ;3 2 4 ; 3 2 0;];
    b=[20;42;30];
    Aeq=[]; beq=[];
    vlb=[0;0;0;]; vub=[];
[x,fval]=linprog(f,A,b,Aeq,beq,vlb,vub);
disp(x)
disp(-fval)

输出结果：

Optimal solution found.

         0
   15.0000
    3.0000

78

例2.求解线性规划

输入代码：

f=[6 4];
    A=[2 3 ;4 2 ;];
    b=[100;120];
    Aeq=[]; beq=[];
    vlb=[0;0;]; vub=[];
[x,fval]=linprog(-f,A,b,Aeq,beq,vlb,vub);
disp(x)
disp(-fval)

输出结果：

Optimal solution found.

   20.0000
   20.0000

   200

合作博弈模型

此处的合作博弈模型采用shapley向量法对问题进行求解。

例3.shapley法

有三个村庄合作在一个附近山峰上修建了一座电视接收塔，并用通信电缆将信号传送到三个村庄。设电视接收塔所在地为0，三个村庄的所在地分别为1,2,3.任意两地之间架设通信电缆的成本费用标注如下图（单位：万元）。

若暂不考虑电视接收塔和信号中转站的成本费用。只考虑通信电缆的架设成本，我们可以清楚地看到，只须假设（0,1）之间、（1,2）之间和（1,3）之间架电缆，就可以满足传输信号的要求，总成本是25万元。现在问题是，这25万元的总成本应该在三个村庄之间如何分摊？

首先，我们需要构造一个shapley函数方便调用计算：

%% shapley函数
%% shapley函数
function shapley = shapley(s,v)
%% 计算参与者数量
n = log2((length(v)+1));
 
%% 计算w(|s|)
s=s';
s_1=sum(s);
s_2=sum(s~=0, 2);
s_all_1 = unique(s_1);
s_all_2 = unique(s_2);
w=zeros(max(s_all_1),max(s_all_2));
v_1=zeros(max(s_all_1),max(s_all_2));
v_2=zeros(max(s_all_1),max(s_all_2));
[s_row, ~] = size(s);
for i=1:1:s_row                                         %得到|s|矩阵
    A=s;
    A(:, ~any(A(i,:),1)) = [];                          %删除每行的含0列
    A_1=sum(A);                                         %计数分类
    [unique_values, ~, unique_counts] = unique(A_1);    %[计数值，计数值重复行]
    counts = histcounts(unique_counts, 1:length(unique_values)+1);
    A_1_repeated_counts = repelem(counts, counts);
    [~, A_1_repeated_counts_col] = size(A_1_repeated_counts);
    [~, counts_col] = size(counts);
    w(i,:)=(ones(1,A_1_repeated_counts_col)./A_1_repeated_counts)./counts_col;   %得到w矩阵
end
 
%% 运算得到v_1矩阵 v_2矩阵
for i=1:1:s_row                                         %遍历s矩阵生成v_1                                        
    B=s; 
    B(n+1,:)=v;
    B(:, ~any(B(i,:),1)) = [];
    v_1(i,:)=B(n+1,:);                                  %v_1矩阵
    B_1=B;
    [~, B_col] = size(B);
    B_1(i,:)=zeros(1,B_col);
    for j=1:1:s_row                                     %遍历生成C矩阵做比较
        C=s;
        C(n+1,:)=v;
        C(:, ~any(C(j,:),1)) = [];
        [~, C_col] = size(C);
        for k=1:1:C_col                                 %遍历B矩阵和C矩阵比较
            for l=1:1:C_col
                if isequal(B_1(1:s_row, k), C(1:s_row, l))
                    v_2(i,k)=C(n+1,l);                  %遍历得到v_2
                end
            end
        end
    end
end
%% 计算sum_w(|s|)*v
v_3=v_1-v_2;
x=sum(w.*v_3,2);
shapley=x

此处的s,v分别合作关系与特征函数，函数调用过程如下：

s=[1 0 0
    0 1 0
    0 0 1
    1 1 0
    1 0 1
    0 1 1
    1 1 1];
v=[0 0 0 9 11 10 20];
shapley=shapley(s,v);
chengben=[10 15 20];
result=chengben-shapley'

于是得到最终输出结果：

shapley =

    6.6667
    6.1667
    7.1667


result =

3.3333    8.8333   12.8333

所以1，2，3村分别摊费用为：3.3333万元、8.8333万元、12.8333万元

此外，提供一个多主体的示例

 s=[1 0 0 0 0 0
    0 1 0 0 0 0
    0 0 1 0 0 0
    0 0 0 1 0 0
    0 0 0 0 1 0
    0 0 0 0 0 1
    1 1 0 0 0 0 
    1 0 1 0 0 0
    1 0 0 1 0 0
    1 0 0 0 1 0
    1 0 0 0 0 1
    0 1 1 0 0 0
    0 1 0 1 0 0
    0 1 0 0 1 0
    0 1 0 0 0 1
    0 0 1 1 0 0
    0 0 1 0 1 0
    0 0 1 0 0 1
    0 0 0 1 1 0
    0 0 0 1 0 1
    0 0 0 0 1 1
    1 1 1 0 0 0
    1 1 0 1 0 0
    1 1 0 0 1 0
    1 1 0 0 0 1
    1 0 1 1 0 0
    1 0 1 0 1 0
    1 0 1 0 0 1
    1 0 0 1 1 0
    1 0 0 1 0 1
    1 0 0 0 1 1
    0 1 1 1 0 0
    0 1 1 0 1 0
    0 1 1 0 0 1
    0 1 0 1 1 0
    0 1 0 1 0 1
    0 1 0 0 1 1
    0 0 1 1 1 0
    0 0 1 1 0 1
    0 0 1 0 1 1
    0 0 0 1 1 1
    1 1 1 1 0 0
    1 1 1 0 1 0
    1 1 1 0 0 1
    1 1 0 1 1 0
    1 1 0 1 0 1
    1 1 0 0 1 1
    1 0 1 1 1 0
    1 0 1 1 0 1
    1 0 1 0 1 1
    1 0 0 1 1 1
    0 1 1 1 1 0
    0 1 1 1 0 1
    0 1 1 0 1 1
    0 1 0 1 1 1
    0 0 1 1 1 1
    0 1 1 1 1 1
    1 0 1 1 1 1
    1 1 0 1 1 1
    1 1 1 0 1 1
    1 1 1 1 0 1
    1 1 1 1 1 0
    1 1 1 1 1 1];
v=[0 0 0 0 0 0 0.2 0.2 0.2 0.2 0.2 0.4 0.4 0.4 0.4 0.6 0.6 0.6 0.8 0.8 1 0.6 0.6 0.6 0.6 0.8 0.8 0.8 1 1 1.2 1 1 1 1.2 1.2 1.4 1.4 1.4 1.6 1.8 1.2 1.2 1.2 1.4 1.4 1.6 1.6 1.6 1.8 2 1.8 1.8 2 2.2 2.4 2.8 2.6 2.4 2.2 2 2 3];

shapley =

    0.1667
    0.3267
    0.4767
    0.6100
    0.7100
    0.7100


sum(shapley)

ans =

     3

欢迎留言讨论，交流学习。