纳什均衡求解器

1 实验目的

• 掌握求解纳什均衡的相关算法
• 锻炼数学基础能力及编程求解问题的能力

2 实验内容

本次实验要求使用Python语言在给定代码框架下编程求解纳什均衡 (Nash Equilibrium, NE), 包括纯策略 NE 与混合策略 NE, 并提交相应源码、输出文件以及实验报告。

3 实现过程

本次大作业要求我们用Python语言写一个纳什均衡求解器，除 Python 标准库外, 仅允许额外引入Numpy和Scipy两个包。Numpy是一个用于科学计算的基础包，使用它可以提高我们代码运行的效率；Scipy是一个用于执行数学、科学和工程计算的python开源代码，我们会用到它的 optimize模块求解线性规划问题，用到它的 spatial模块求解空间几何问题。下面将简要介绍纳什均衡求解器的实现过程，并主要聚焦于文件读取、纯策略均衡求解和混合策略均衡求解。

3.1 文件读取

每个博弈数据存储在以.nfg为后缀的文件中，我们需要读取这个文件并解析它。经观察发现，文件的最后一行总是payoff列表，倒数第三行总是包含了每个玩家可选的策略数目，根据这两行内容，我们就可以获取关于这个博弈的所有信息。

根据Gambit format的规定，玩家可选的策略数目总是包含在一对花括号里，并用空格分隔，比如{ 3 2 5 }表示玩家1有3种策略，玩家2有2种策略，玩家3有5种策略。知道了它的表示形式后，我们可以使用正则表达式把每个玩家可选的策略数目提取出来，同时我们也相当于知道了玩家的数量。

根据Gambit format的规定，如果总共有N N N个玩家，那么payoff列表中的前N N N个数表示这N N N个玩家都选第一种策略时，每个玩家的收益分别是多少。第一个数对应于玩家1的收益，第二个数对应于玩家2的收益，以此类推。第N + 1 N+1 N +1个数到第2 N 2N 2 N个数表示玩家1选第二种策略，其余玩家选第一种策略时，每个玩家的收益。根据约定，首先是第一个玩家的策略递增，当它递增完自己所有可选择的策略之后，它的策略会”回滚”到1，然后第二个玩家的策略依次递增。在我的实现中，我选择了将每个玩家的收益从总的payoff列表中分离出来。分离的方法也很简单，对于payoff列表中第i i i个数（从0开始编号），它对应于玩家i % N i \% N i %N（玩家从0开始编号）的收益。现在分离出来的每个玩家的收益仍然是列表的形式，我们需要将其表示为矩阵的形式。

对于有N N N个玩家的博弈，每个玩家的收益都应该是一个N N N维的矩阵。按照前面说过的收益列表与策略的对应关系，每个玩家收益列表中的第一个数，对应于所有玩家都选第一种策略的情况，在收益矩阵中的坐标应该是(1, 1, ⋯ \cdots ⋯, 1)。第二个数对应于玩家1选第二种策略，其余玩家选第一种策略的情况，在收益矩阵中的坐标应该是(2, 1, ⋯ \cdots ⋯, 1)。如果玩家1有s 1 s_1 s 1 ​种策略，那么收益列表中的第s 1 + 1 s_{1}+1 s 1 ​+1个数对应于玩家2选策略2，其余玩家选策略1的情况，在收益矩阵中的坐标应该是(1, 2, ⋯ \cdots ⋯, 1)。收益列表中第s 1 + 2 s_{1}+2 s 1 ​+2个数对应于玩家1选策略2，玩家2选策略2，其余玩家选策略1的情况，在收益矩阵中的坐标应该是(2, 2, ⋯ \cdots ⋯, 1)，第s 1 + 3 s_1 + 3 s 1 ​+3个数对应于收益矩阵中的坐标应该是(3, 2, ⋯ \cdots ⋯, 1)。可以看出，如果对收益矩阵按照高维优先访问，那么就可以得到收益列表。同理，从收益列表转换到收益矩阵也需要按照高维优先的原则。在Numpy的reshape函数中，提供了一个名为order的参数，当指定其为C时，它会按照C语言的方式，即低维优先进行转换，当指定其为F时，它会按照Fortran的方式，即高维优先进行转换。因此，我们只需要指定order为F就可以很方便的将收益列表转换为收益矩阵。请注意，任意一个收益矩阵中的第k k k维对应于玩家k k k的策略选择。

至此，我们已经从.nfg文件中读取并解析出了每个玩家可选择的策略数，以及每个玩家的收益矩阵。

3.2 纯策略求解

当有三个及以上的玩家时，我们仅需要求解它所有的纯策略纳什均衡点。在上一小节中我们已经介绍过，每个玩家都有自己的收益矩阵，并且收益矩阵中的第i i i维对应于第i i i个玩家的策略。假如对于玩家i i i，我们想要计算在其余玩家选定策略后，它选择哪个策略可以达到最大收益，只需计算玩家i i i的收益矩阵中第i i i维的最大值即可。Numpy自带了求最大值位置的函数 argmax，要求在第i i i维上的最大值，只需要指定参数axis=i i i即可，但是此函数却并不适用于本次任务。我们的目标是求所有最大值的位置，当有多个最大值时，要全部计算出来，而 argmax只会返回第一个最大值的位置，这显然不符合我们的要求。经过查阅资料，我最终采用了将 max函数和 where函数结合起来的方案。首先用 max函数求第i i i维上的最大值，然后再用 where函数找出所有与最大值相同的元素的位置。

当对所有的玩家i ∈ [ N ] i \in [N]i ∈[N ]，其中N N N表示玩家数量，都找出了它在第i i i维上的最大值的坐标集合后，将所有玩家的坐标集合求交集即为纯策略纳什均衡点。

3.3 混合策略求解

当只有两个人进行博弈时，我们需要求解它所有的纯策略和混合策略纳什均衡点（无穷多个时求极点）。在我的求解中，参考了Algorithmic Game Theory中3.3节的Labeled polytopes算法。下面将先简要介绍此算法，然后强调几个实现过程中的重点。

3.3.1 Labeled polytopes算法理论

在求解混合策略纳什均衡时，我们需要知道均衡策略中的哪些策略的概率不为0，而最优反应多面体（best response polyhedron）就是一种非常直接的求解方法。记x \boldsymbol{x}x和y \boldsymbol{y}y分别表示玩家1和玩家2的混合策略，A \boldsymbol{A}A和B \boldsymbol{B}B分别表示玩家1和玩家2的收益矩阵，u u u和v v v分别表示玩家1和玩家2期望收益的上包络（upper envelope），M M M和N N N分别表示玩家1和玩家2可选择的策略数量。玩家1和玩家2的最优反应多面体P ˉ \bar{P}P ˉ和Q ˉ \bar{Q}Q ˉ​分别被定义为
P ˉ = { ( x , v ) ∈ R M × R ∣ x ≥ 0 , 1 ⊤ x = 1 , B ⊤ x ≤ 1 v } Q ˉ = { ( y , u ) ∈ R N × R ∣ A y ≤ 1 u , y ≥ 0 , 1 ⊤ y = 1 } \begin{aligned} & \bar{P}=\left{(\boldsymbol{x}, v) \in \mathbb{R}^{M}\times \mathbb{R} \mid \boldsymbol{x} \geq \mathbf{0},\ \mathbf{1}^{\top} \boldsymbol{x}=1,\ \boldsymbol{B}^{\top} \boldsymbol{x} \leq \mathbf{1} v\right} \ & \bar{Q}=\left{(\boldsymbol{y}, u) \in \mathbb{R}^{N} \times \mathbb{R} \mid \boldsymbol{A} \boldsymbol{y} \leq \mathbf{1} u,\ \boldsymbol{y} \geq \mathbf{0},\ \mathbf{1}^{\top} \boldsymbol{y}=1\right} \end{aligned}​P ˉ={(x ,v )∈R M ×R ∣x ≥0 ,1 ⊤x =1 ,B ⊤x ≤1 v }Q ˉ​={(y ,u )∈R N ×R ∣A y ≤1 u ,y ≥0 ,1 ⊤y =1 }​
可以看出，P ˉ \bar{P}P ˉ中是一些形如( x 1 , x 2 , ⋯ , x M , v ) (x_1, x_2, \cdots, x_M, v)(x 1 ​,x 2 ​,⋯,x M ​,v )的元组，Q ˉ \bar{Q}Q ˉ​中是一些形如( y M + 1 , y M + 2 , ⋯ , y M + N , u ) (y_{M+1}, y_{M+2}, \cdots, y_{M+N}, u)(y M +1 ​,y M +2 ​,⋯,y M +N ​,u )的元组，其中x i x_i x i ​和y j y_j y j ​分别表示玩家1选择策略i i i和玩家2选择策略j j j的概率。下面根据作业指导手册上的例子，对P ˉ \bar{P}P ˉ和Q ˉ \bar{Q}Q ˉ​中的不等式约束做进一步的说明。玩家1和玩家2的收益矩阵分别如表1和表2所示。

玩家1有3种策略，即x = ( x 1 , x 2 , x 3 ) \boldsymbol{x} = (x_1, x_2, x_3)x =(x 1 ​,x 2 ​,x 3 ​)，玩家2有两种策略，即y = ( y 4 , y 5 ) \boldsymbol{y} = (y_4, y_5)y =(y 4 ​,y 5 ​)。对于P ˉ \bar{P}P ˉ来说，x ≥ 0 \boldsymbol{x} \ge \boldsymbol{0}x ≥0可以产生三个不等式，B ⊤ x ≤ 1 v \boldsymbol{B}^{\top} \boldsymbol{x} \leq \mathbf{1} v B ⊤x ≤1 v可以产生两个不等式，因此P ˉ \bar{P}P ˉ的所有不等式约束为
x 1 ≥ 0 x 2 ≥ 0 x 3 ≥ 0 1 ⋅ x 1 + 2 ⋅ x 2 + 2 ⋅ x 3 ≤ v 1 ⋅ x 1 + 3 ⋅ x 2 + 0 ⋅ x 3 ≤ v \begin{aligned} \nonumber x_1 &\ge 0\ x_2 &\ge 0\ x_3 &\ge 0\ 1\cdot x_1 + 2 \cdot x_2 + 2 \cdot x_3 &\le v\ 1\cdot x_1 + 3\cdot x_2 + 0\cdot x_3 &\le v \end{aligned}x 1 ​x 2 ​x 3 ​1 ⋅x 1 ​+2 ⋅x 2 ​+2 ⋅x 3 ​1 ⋅x 1 ​+3 ⋅x 2 ​+0 ⋅x 3 ​​≥0 ≥0 ≥0 ≤v ≤v ​

对于Q ˉ \bar{Q}Q ˉ​来说，A y ≤ 1 u \boldsymbol{A} \boldsymbol{y} \leq \mathbf{1} u A y ≤1 u可以产生三个不等式，y ≥ 0 \boldsymbol{y} \ge \boldsymbol{0}y ≥0可以产生两个不等式，因此Q ˉ \bar{Q}Q ˉ​的所有不等式约束为
1 ⋅ y 4 + 1 ⋅ y 5 ≤ u 0 ⋅ y 4 + 0 ⋅ y 5 ≤ u 0 ⋅ y 4 + 2 ⋅ y 5 ≤ u y 4 ≥ 0 y 5 ≥ 0 \begin{aligned} 1\cdot y_4 + 1\cdot y_5 &\le u\ 0\cdot y_4 + 0\cdot y_5 &\le u\ 0\cdot y_4 + 2 \cdot y_5 &\le u\ y_4 &\ge 0\ y_5 &\ge 0 \end{aligned}1 ⋅y 4 ​+1 ⋅y 5 ​0 ⋅y 4 ​+0 ⋅y 5 ​0 ⋅y 4 ​+2 ⋅y 5 ​y 4 ​y 5 ​​≤u ≤u ≤u ≥0 ≥0 ​

下面引入标记的概念。总的来说，如果某个点被标记为l l l，那么这个点能使第l l l个不等式中的等号成立。因为P ˉ \bar{P}P ˉ和Q ˉ \bar{Q}Q ˉ​中不等式的顺序不同，具体来说，在P ˉ \bar{P}P ˉ中，是先x ≥ 0 \boldsymbol{x} \ge \boldsymbol{0}x ≥0，后B ⊤ x ≤ 1 v \boldsymbol{B}^{\top} \boldsymbol{x} \leq \mathbf{1} v B ⊤x ≤1 v，而在Q ˉ \bar{Q}Q ˉ​中，是先A y ≤ 1 u \boldsymbol{A} \boldsymbol{y} \leq \mathbf{1} u A y ≤1 u，后y ≥ 0 \boldsymbol{y} \ge \boldsymbol{0}y ≥0，所以标记的含义也有所区别。令M = { 1 , 2 , 3 } M = {1, 2, 3}M ={1 ,2 ,3 }表示玩家1的标记集合，N = { 4 , 5 } N = {4, 5}N ={4 ,5 }表示玩家2的标记集合。对于某个点( y , u ) ∈ Q ˉ (\boldsymbol{y}, u) \in \bar{Q}(y ,u )∈Q ˉ​，它的标记k ∈ M ∪ N k \in M \cup N k ∈M ∪N，如果k = i ∈ M k=i \in M k =i ∈M，那么说明第i i i个不等式等号成立，即∑ j ∈ N a i j y j = u \sum_{j \in N} a_{i j} y_{j}=u ∑j ∈N ​a ij ​y j ​=u；如果k = j ∈ N k=j \in N k =j ∈N，那么说明第j j j个不等式等号成立，即y j = 0 y_j = 0 y j ​=0。对于某个点( x , v ) ∈ P ˉ (\boldsymbol{x}, v) \in \bar{P}(x ,v )∈P ˉ，它的标记k ∈ M ∪ N k \in M \cup N k ∈M ∪N，如果k = i ∈ M k = i \in M k =i ∈M，那么x i = 0 x_i = 0 x i ​=0，如果k = j ∈ N k = j \in N k =j ∈N，那么∑ i ∈ M b i j x i = v \sum_{i \in M} b_{i j} x_{i}=v ∑i ∈M ​b ij ​x i ​=v。

对于( ( x , v ) , ( y , u ) ) ∈ P ˉ × Q ˉ ((\boldsymbol{x}, v), (\boldsymbol{y}, u)) \in \bar{P} \times \bar{Q}((x ,v ),(y ,u ))∈P ˉ×Q ˉ​而言，当它是被完全标记时，( x , y ) (\boldsymbol{x}, \boldsymbol{y})(x ,y )是纳什均衡点。完全标记意味着对于任意一个标记k ∈ M ∪ N k \in M \cup N k ∈M ∪N，k k k要么出现在( x , v ) (\boldsymbol{x}, v)(x ,v )中，要么出现在( y , u ) (\boldsymbol{y} ,u)(y ,u )中（可以同时出现）。

从上面的算法描述中可以看出，计算P ˉ \bar{P}P ˉ和Q ˉ \bar{Q}Q ˉ​需要知道u u u和v v v，而当收益矩阵A \boldsymbol{A}A和B ⊤ \boldsymbol{B}^{\top}B ⊤非负时，我们可以简化掉u u u和v v v。得到的新的最优反应多面体为
P = { x ∈ R M ∣ x ≥ 0 , B ⊤ x ≤ 1 } Q = { y ∈ R N ∣ A y ≤ 1 , y ≥ 0 } \begin{aligned} &P=\left{\boldsymbol{x} \in \mathbb{R}^{M} \mid \boldsymbol{x} \geq \mathbf{0}, \boldsymbol{B}^{\top} \boldsymbol{x} \leq \mathbf{1}\right} \ &Q=\left{\boldsymbol{y} \in \mathbb{R}^{N} \mid \boldsymbol{A} \boldsymbol{y} \leq \mathbf{1}, \quad \boldsymbol{y} \geq \mathbf{0}\right} \end{aligned}​P ={x ∈R M ∣x ≥0 ,B ⊤x ≤1 }Q ={y ∈R N ∣A y ≤1 ,y ≥0 }​

和原先的相比，去掉了等值约束，并将上包络正则化为1，此时集合中的元素为相应玩家的混合策略（未正则化）。在使用P P P和Q Q Q时，我们需要去掉其中的0向量，对于非零的向量x ∈ P \boldsymbol{x} \in P x ∈P或y ∈ Q \boldsymbol{y} \in Q y ∈Q，我们可以通过除以1 ⊤ x \boldsymbol{1}^{\top} \boldsymbol{x}1 ⊤x或1 ⊤ y \boldsymbol{1}^{\top} \boldsymbol{y}1 ⊤y将其正则化为概率向量。

上面的P P P和Q Q Q要求收益矩阵A \boldsymbol{A}A和B ⊤ \boldsymbol{B}^{\top}B ⊤非负，而当不满足这个条件时，可以通过平移收益矩阵，即给矩阵中的每一个元素加上一个大于0的常数，使其满足非负性的条件，并且这样做不会改变纳什均衡的求解结果。因此，对于任意的收益矩阵，都可以用于计算简化后的最优反应多面体。

至此，Labeled polytopes的主要算法就介绍完了。当计算出了P P P和Q Q Q之后，我们可以枚举这两个集合中的元素（x = 0 \boldsymbol{x} = \boldsymbol{0}x =0或y = 0 \boldsymbol{y} = \boldsymbol{0}y =0的情况除外）来找到均衡点。对于x ∈ P − { 0 } \boldsymbol{x} \in P – {\boldsymbol{0}}x ∈P −{0 }和y ∈ Q − { 0 } \boldsymbol{y} \in Q – {\boldsymbol{0}}y ∈Q −{0 }，如果( x , y ) (\boldsymbol{x}, \boldsymbol{y})(x ,y )是被完全标记的，那么( x / ( 1 ⊤ x ) , y / ( 1 ⊤ y ) ) \left(\boldsymbol{x} / (\boldsymbol{1}^{\top} \boldsymbol{x}), \boldsymbol{y} /(\boldsymbol{1}^{\top} \boldsymbol{y}) \right)(x /(1 ⊤x ),y /(1 ⊤y ))即为混合策略纳什均衡。

; 3.3.2 Labeled polytopes算法实现

对于一个二人博弈，我们编写代码的重点是如何求解P P P和Q Q Q，以及这两个集合里每个元素的标记。显然，对于P P P或Q Q Q的求解是一个线性规划问题。为了能方便的使用现有的库函数，我们对P P P和Q Q Q的定义做一个等价转化，即都转换成小于等于的形似：
P = { x ∈ R M ∣ − x ≤ 0 , B ⊤ x ≤ 1 } Q = { y ∈ R N ∣ A y ≤ 1 , − y ≤ 0 } \begin{aligned} &P=\left{\boldsymbol{x} \in \mathbb{R}^{M} \mid -\boldsymbol{x} \le \mathbf{0}, \boldsymbol{B}^{\top} x \leq \mathbf{1}\right} \ &Q=\left{\boldsymbol{y} \in \mathbb{R}^{N} \mid \boldsymbol{A} \boldsymbol{y} \leq \mathbf{1}, \quad -\boldsymbol{y} \le \mathbf{0}\right} \end{aligned}​P ={x ∈R M ∣−x ≤0 ,B ⊤x ≤1 }Q ={y ∈R N ∣A y ≤1 ,−y ≤0 }​

仍然以表1和表2表示的博弈为例，此时P P P中的不等式约束变为
− x 1 ≤ 0 − x 2 ≤ 0 − x 3 ≤ 0 1 ⋅ x 1 + 2 ⋅ x 2 + 2 ⋅ x 3 ≤ 1 1 ⋅ x 1 + 3 ⋅ x 2 + 0 ⋅ x 3 ≤ 1 (1) \begin{aligned} -x_1 &\le 0\ -x_2 &\le 0\ -x_3 &\le 0\ 1\cdot x_1 + 2 \cdot x_2 + 2 \cdot x_3 &\le 1\ 1\cdot x_1 + 3\cdot x_2 + 0\cdot x_3 &\le 1 \end{aligned} \tag{1}−x 1 ​−x 2 ​−x 3 ​1 ⋅x 1 ​+2 ⋅x 2 ​+2 ⋅x 3 ​1 ⋅x 1 ​+3 ⋅x 2 ​+0 ⋅x 3 ​​≤0 ≤0 ≤0 ≤1 ≤1 ​(1 )
Q Q Q中的不等式约束变为
1 ⋅ y 4 + 1 ⋅ y 5 ≤ 1 0 ⋅ y 4 + 0 ⋅ y 5 ≤ 1 0 ⋅ y 4 + 2 ⋅ y 5 ≤ 1 − y 4 ≤ 0 − y 5 ≤ 0 (2) \begin{aligned} 1\cdot y_4 + 1\cdot y_5 &\le 1\ 0\cdot y_4 + 0\cdot y_5 &\le 1\ 0\cdot y_4 + 2 \cdot y_5 &\le 1\ -y_4 &\le 0\ -y_5 &\le 0 \end{aligned} \tag{2}1 ⋅y 4 ​+1 ⋅y 5 ​0 ⋅y 4 ​+0 ⋅y 5 ​0 ⋅y 4 ​+2 ⋅y 5 ​−y 4 ​−y 5 ​​≤1 ≤1 ≤1 ≤0 ≤0 ​(2 )
为了表述方便，下面以P P P为例来介绍求解过程。P P P里面的不等式约束构成了一个空间，现在我们要计算此空间的极点。在我的实现中，极点的计算是由 scipy.spatial模块下的HalfspaceIntersection函数完成的。这个函数有两个必要的参数，第一个是halfspaces，它是一个二维的矩阵，形状为(nineq, ndim+1)，它是M x + b ≤ 0 \boldsymbol{M}\boldsymbol{x}+\boldsymbol{b}\le \boldsymbol{0}M x +b ≤0这种不等式以[ M ; b ] [\boldsymbol{M};\boldsymbol{b}][M ;b ]这种形式表示的系数，其中nineq表示不等式的数目，它等于M \boldsymbol{M}M的第一维的大小，ndim表示x \boldsymbol{x}x的长度，它等于M \boldsymbol{M}M的第二维的大小。比如对于P P P中的不等式约束(1)，它的halfspaces为
[ − 1 0 0 0 0 − 1 0 0 0 0 − 1 0 1 2 2 − 1 1 3 0 − 1 ] \begin{bmatrix} -1 & 0 & 0 & 0\ 0 & -1 & 0 & 0\ 0 & 0 & -1 & 0\ 1 & 2 & 2 & -1\ 1 & 3 & 0 & -1 \end{bmatrix}​−1 0 0 1 1 ​0 −1 0 2 3 ​0 0 −1 2 0 ​0 0 0 −1 −1 ​​

在由收益矩阵生成halfspaces时需要注意，P P P和Q Q Q的生成方式不同。P P P是混合策略（即x \boldsymbol{x}x）小于等于0的约束在前，而Q Q Q是混合策略（即y \boldsymbol{y}y)小于等于0的约束在后。因此，不能使用完全相同的代码来生成P P P和Q Q Q的halfspaces，而是要根据情况来生成。

另一个参数是interior_point，它表示已经确定在halfspaces所定义的区域内部的点，也被称为可行点，可以通过线性规划求得。在HalfspaceIntersection的官方文档中，提供了一种求interior_point的{方法}。考虑具有M x + b ≤ 0 \boldsymbol{M}\boldsymbol{x} + \boldsymbol{b} \le \boldsymbol{0}M x +b ≤0形式的半空间，解决线性规划问题
max ⁡ y s . t . M x + y ∣ ∣ M i ∣ ∣ ≤ − b \max_{y} \quad s.t.\quad \boldsymbol{M}\boldsymbol{x} + y ||\boldsymbol{M}_i|| \le – \boldsymbol{b}y max ​s .t .M x +y ∣∣M i ​∣∣≤−b

其中M i \boldsymbol{M}_i M i ​表示M \boldsymbol{M}M的第i i i行，即每一行的范式。通过这种方法，可以产生一个凸多面体最内部的点x \boldsymbol{x}x。

上面的线性规划求极值问题由scipy.optimize模块的linprog函数完成。这个函数用于在等式和不等式约束下求一个线性函数的最小值。形式化的来说，对于如下形式的问题，
min ⁡ x c T x s . t . A u b x ≤ b u b A e q x = b e q l ≤ x ≤ u \begin{aligned} \min {x}\,& c^{T} x \ s.t.\quad &A{u b} x \leq b_{u b} \ &A_{e q} x=b_{e q} \ &l \leq x \leq u \end{aligned}x min ​s .t .​c T x A u b ​x ≤b u b ​A e q ​x =b e q ​l ≤x ≤u ​

其中x x x是决策变量，c , b u b , b e q , l c, b_{ub}, b_{eq}, l c ,b u b ​,b e q ​,l和u u u都是向量，A u b A_{ub}A u b ​和A e q A_{eq}A e q ​都是矩阵，它可以求c ⊤ x c^{\top} {x}c ⊤x的最小值。如果想要求一个线性函数的最大值，只需要将其加个负号就可以转化成极小值的问题。在HalfspaceIntersection的官方文档中，提供了求解interior_point的代码，只需要稍作修改就可以应用于纳什均衡求解器。

求得了约束空间的极点之后，接下来需要计算每一个非0极点的标记。对于halfspaces为[ M , b ] [\boldsymbol{M}, \boldsymbol{b}][M ,b ]的非0极点z \boldsymbol{z}z，M z − b \boldsymbol{M}\boldsymbol{z} – \boldsymbol{b}M z −b会得到一个向量，这个向量中值为0的元素所在的位置就是z \boldsymbol{z}z的标记。在实际编写代码时发现，由于浮点数计算精度的问题，有时候理论上为0的元素实际上并不为0，而是一个非常小的数，比如e − 17 e^{-17}e −17，因此，在我的代码中使用了np.isclose(M z \boldsymbol{M}\boldsymbol{z}M z, − b – \boldsymbol{b}−b)取代M z \boldsymbol{M}\boldsymbol{z}M z == − b – \boldsymbol{b}−b。

现在我们已经计算出了P P P和Q Q Q中的极点及其标记，接下来就是枚举这两个集合，找到完全标记的组合，此部分实现较为简单，不再详述。

4 实验结果

以examples文件夹下的.nfg文件为样例输入，.ne文件为样例输出，对代码进行了测试。测试结果显示，程序的输出结果和样例输出完全一致，这说明了代码的正确性。

为了测试代码的运行效率，以input文件夹下的.nfg文件为输入，重复运行程序10次，取运行时间的平均值来作为度量。对于32个纳什均衡求解问题，最终测得的平均运行时间为25.17秒。考虑到写文件操作可能会比较耗时，无法真正的展现求解纳什均衡的运行时间，因此上述结果不包含写文件的耗时。

5 源代码

import re
import numpy as np
from scipy.optimize import linprog
from scipy.spatial import HalfspaceIntersection
import os

def write_rst(equilibria, out_path):
"""
    将求得的纳什均衡点写入文件中
    :param equilibria: 所有纳什均衡点
    :param out_path: 要写入的文件路径
    :return:
"""
    with open(out_path, 'wt') as f:
        for eq in equilibria:
            write_content = ''
            for idx, value in enumerate(eq):
                if type(value) is int:
                    value = str(abs(value))
                else:
                    value = str(abs(value))

                if write_content == '':
                    write_content = value
                else:
                    write_content += ',' + value
            f.write(write_content + '\n')

def read_payoff(in_path):
"""
    读取.nfg文件，并将列表形式的收益解析成矩阵形式
    :param in_path: .nfg文件的路径
    :return: 每个玩家矩阵形式的收益，每个玩家可选择的策略数量
"""
    with open(in_path) as fin:
        all_lines = fin.readlines()
        payoff = all_lines[-1].strip().split()
        player_info = all_lines[-3]
    res = re.search('{ (\d+ )*}', player_info)
    strategy_space = [int(num) for num in res.group()[1: -1].split()]
    player_num = len(strategy_space)
    all_payoff = []
    for player in range(player_num):
        payoff_matrix = [int(pay) for idx, pay in enumerate(payoff) if idx % player_num == player]
        payoff_matrix = np.array(payoff_matrix).reshape(strategy_space, order='F')
        all_payoff.append(payoff_matrix)
    return all_payoff, strategy_space

def create_half_spaces(payoff_matrix, is_p=False):
"""
    根据收益矩阵创建half_spaces
    :param payoff_matrix: 收益矩阵，正是根据它来创建half_spaces
    :param is_p: 根据参考书上的算法，polyhedron P和Q中half_spaces的顺序不一样，此变量用于确定是哪一种情况
    :return: 创建好的half_spaces，加入payoff_matrix的维度维d1 * d2，那么half_spaces的维度就是(d1+d2) * (d2 + 1)
"""
    d1, d2 = payoff_matrix.shape
    if is_p:
        b = np.append(np.zeros(d2), -np.ones(d1))
        payoff_matrix = np.append(-np.eye(d2), payoff_matrix, axis=0)
    else:
        b = np.append(-np.ones(d1), np.zeros(d2))
        payoff_matrix = np.append(payoff_matrix, -np.eye(d2), axis=0)

    half_spaces = np.column_stack((payoff_matrix, b.transpose()))
    return half_spaces

def find_feasible_point(half_spaces):
"""
    求线性规划空间内部的点，参考了HalfspaceIntersection官方文档，
    地址：https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.HalfspaceIntersection.html?highlight=halfspaceintersection#scipy.spatial.HalfspaceIntersection
    :param half_spaces:
    :return: half_spaces所表示的线性空间的内部的点
"""
    norm_vector = np.reshape(
        np.linalg.norm(half_spaces[:, :-1], axis=1), (half_spaces.shape[0], 1)
    )
    c = np.zeros((half_spaces.shape[1],))
    c[-1] = -1
    A_ub = np.hstack((half_spaces[:, :-1], norm_vector))
    b_ub = -half_spaces[:, -1:]
    res = linprog(c, A_ub=A_ub, b_ub=b_ub)
    return res.x[:-1]

def compute_vertices_and_labels(half_spaces):
"""
    根据half_spaces来计算vertices，需要注意结果中不能含有0向量
    :param half_spaces:
    :return:
"""
    feasible_point = find_feasible_point(half_spaces)
    hs = HalfspaceIntersection(half_spaces, feasible_point)
    hs.close()
    rst = []
    for v in hs.intersections:
        if not np.all(np.isclose(v, 0)):

            b = half_spaces[:, -1]
            A = half_spaces[:, :-1]
            labeled_v = np.where(np.isclose(np.dot(A, v), -b))[0]
            rst.append((v, labeled_v))
    return rst

def mixed_nash_equilibrium(all_payoff):
"""
    计算两个玩家的混合策略纳什均衡（含纯策略）
    :param all_payoff: 两个玩家的收益矩阵
    :return: 求得的所有纯策略纳什均衡点
"""
    payoff_mat1, payoff_mat2 = all_payoff

    if np.min(payoff_mat1) < 0:
        payoff_mat1 = payoff_mat1 + abs(np.min(payoff_mat1))
    if np.min(payoff_mat2) < 0:
        payoff_mat2 = payoff_mat2 + abs(np.min(payoff_mat2))

    p1_strategy_num, p2_strategy_num = payoff_mat1.shape
    all_strategy_num = p1_strategy_num + p2_strategy_num

    p1_half_spaces = create_half_spaces(payoff_mat2.transpose(), is_p=True)
    p2_half_spaces = create_half_spaces(payoff_mat1, is_p=False)
    rst = []
    p_vertices = compute_vertices_and_labels(p1_half_spaces)
    q_vertices = compute_vertices_and_labels(p2_half_spaces)
    for vec_x, x_l in p_vertices:
        all_labels = np.zeros(all_strategy_num, dtype=int)
        all_labels[x_l] = 1
        for vec_y, y_l in q_vertices:
            all_labels[y_l] += 1
            if np.all(all_labels):
                rst.append(np.append(vec_x / sum(vec_x), vec_y / sum(vec_y)))
            all_labels[list(y_l)] -= 1
    return rst

def compare_eq(computed_NE, expected_NEs):
    '''
    用于测试求解结果是否正确
    :param computed_NE: 通过程序计算出来的均衡点
    :param expected_NEs: 正确的纳什均衡点
    :return: 如果计算出来的NE是正确的，则返回True，否则返回False
    '''
    for one_NE in expected_NEs:
        diff = np.abs(computed_NE - one_NE)
        if np.sum(diff)  1e-4:
            return True
    return False

def pure_nash_equilibrium(all_payoff):
"""
    计算三个及以上的玩家的混合策略纳什均衡
    :param all_payoff: 是一个list，第i个元素是第i个玩家的收益矩阵
    :return: 求得的所有混合策略纳什均衡点
"""
    player_num = len(all_payoff)
    strategy_space = all_payoff[0].shape
    max_payoff_coordinate = set()
    for player, payoff_matrix in enumerate(all_payoff):

        temp_max_payoff_coordinate = set()
        max_pos = np.where(payoff_matrix == np.max(payoff_matrix, axis=player, keepdims=True))
        max_pos = list(max_pos)

        for i in range(len(max_pos[0])):
            temp_list = []
            for j in range(player_num):
                temp_list.append(max_pos[j][i])
            temp_max_payoff_coordinate.add(tuple(temp_list))
        if player == 0:
            max_payoff_coordinate = temp_max_payoff_coordinate
        else:
            max_payoff_coordinate = max_payoff_coordinate.intersection(temp_max_payoff_coordinate)

    equilibria = []
    for equilibrium_coord in max_payoff_coordinate:
        equilibrium = [0] * (sum(strategy_space))
        pre_len = 0
        for pos, strategy_num in zip(equilibrium_coord, strategy_space):
            equilibrium[pre_len + pos] = 1
            pre_len += strategy_num
        equilibria.append(equilibrium)
    return equilibria

def nash(in_path, out_path):
"""
    求解纳什均衡，对于两个玩家的博弈，需要求解所有的纯策略均衡点和混合策略均衡点，对于三个及以上玩家的博弈，只需求解纯策略均衡点
    :param in_path: .nfg文件的路径，此文件包含了与博弈有关的所有信息
    :param out_path: 均衡点输出文件的路径
    :return:
"""
    all_payoff, strategy_space = read_payoff(in_path)
    if len(strategy_space) == 2:

        equilibria = mixed_nash_equilibrium(all_payoff)
    else:

        equilibria = pure_nash_equilibrium(all_payoff)

    write_rst(equilibria, out_path)

def main():
    for f in os.listdir('input'):
        if f.endswith('.nfg') and not f.startswith('._'):
            nash('input/' + f, 'output/' + f.replace('nfg', 'ne'))

if __name__ == '__main__':
    main()

实验中用到的测试文件：Github

Original: https://blog.csdn.net/weixin_43074474/article/details/126170200
Author: 王森ouc
Title: 纳什均衡求解器