2020年第九届数学建模国际赛小美赛
D题 石头剪刀游戏与合作
原题再现：
小时候你可能至少玩过几次石头剪刀游戏。在这个游戏中，你几乎有三个选择，每一个都有一个项目要打败，一个项目输给。石头打败剪刀，剪刀剪纸和布覆盖岩石。你所要做的就是选择，看看结果如何。一切似乎都很简单。
从宏观的角度看人类，有许多例子表明人们合作形成社会、国家、宗教和其他团体。
然而，在基本的两人水平上，人们往往会背叛对方，就像在囚徒困境这样的社会困境博弈中发现的那样，即使人们之间合作会得到更好的回报。
合作的主题以及人们如何以及何时开始相互信任已经被许多研究者研究过，他们已经用数字的方法解决了这个问题。在《混沌》杂志的一篇论文中，研究人员分析了合作的驱动因素。
根据经典生态学，当两个物种竞争同一资源时，最终成功的物种会胜出，而另一个物种则会灭绝。但这一规则不能解释像亚马逊这样的生态系统，那里有成千上万的树种占据着相似的生态位。芝加哥大学和加州大学圣巴巴拉分校的研究人员在《美国国家科学院院刊》上发表报告称，儿童时期玩石头剪刀的游戏可以解决这个难题。围绕游戏动态设计的数学模型产生了无限生物多样性的潜力，并提出了一些令人惊讶的新生态规则。
石头剪刀规则是一个“不及物性”竞争的例子，参与者不能简单地从最好的到最坏的排序。当两人一组时，胜利者和失败者就会出现：石头打剪刀，布打石头，剪刀打布。但当这三种策略相互竞争时，就陷入了僵局，没有一个要素是无可争议的赢家。
在自然界中，细菌和蜥蜴的三个物种群都观察到了这种石-布-剪关系。但是，科学家们还没有研究与三个以上的参与者建立更复杂的不及物性关系——想想石头、布、剪刀、炸药等等——如何模拟更复杂的社会困境游戏。
1975年，R.M.May和W.J.Leonard首次使用石头布剪刀游戏来模拟三个物种周期性地相互支配的生态场景：一个物种支配第二个物种，第二个物种支配第三个物种，第三个物种支配第一个物种。例如，这个游戏可以很好地模拟具有周期优势的大肠杆菌的不同菌株。
您的团队已被国际合作组织管理层（ICM）聘为顾问。您应完成以下任务：
•建立一个数学模型来模拟与三个以上参与者的关系，当作为May Leonard模型的扩展案例实施时，该模型的工作原理与原始的石头剪纸模型有所不同。
•通过对模型参数的敏感性分析，设计一个自适应生物系统，其中一个物种可以改变相互作用概率以保证其生存，并探索生物相互作用如何平衡物种之间的不平衡关系以及疾病和其他捕食者的影响。
•为许多其他社会困境游戏开发一些分析解决方案，如石头剪刀、性别之战、猎鹿或探索新冠病毒感染动态。该仿真模型可以帮助理解在许多其他社会困境中合作行为的真正驱动因素。
•向ICM写一份两页的备忘录，描述一个现实的全球目标——可实现的最低防疫合作水平、目前必须解决的一系列问题以及可能加速或阻碍实现目标的任何情况。
提交内容应包括：
•一页汇总表
•目录
•两页备忘录
您的解决方案不超过30页，最多34页，包含您的摘要、目录和两页备忘录。
整体求解过程概述(摘要)
一些科学家使用石头剪刀游戏来模拟三个物种周期性地相互支配的生态场景：一个物种支配第二个物种，第二个物种支配第三个物种，第三个物种支配第一个物种。但在自然界中应该有很多物种（超过三个），科学家们还没有研究过如何更复杂的不及物关系与三个以上的参与者。我们试图建立一个模型来模拟三个以上参与者的关系，并利用该模型设计一个自适应的生物系统，然后分析多个角色的两难博弈，提出可用于抗议合作的解决方案。
对于问题一：在模拟三个或多个参与者之间的关系时，我们发现它不仅与参与者的数量有很大的关系，而且数量是偶数还是奇数也会有很大的影响，所以我们将分别进行研究，主要以N=4和N=5为例进行详细分析，我们扩展了Lotka-Volterra模型，确定空间均匀初始条件（ode模型）、扩散驱动的不同参与者的斑块传播、相互作用和循环动力学（pde模型）的动力学。为了更好地适应背景，我们将玩家视为一个物种，因此我们的模型可以很好地适应任何自适应生物系统。
对于问题二：循环竞争博弈模型，特别是石头剪刀模型，在探索空间生态系统中多物种共存问题中起着重要作用。我们提出了一个扩展的石头剪刀博弈模型来模拟五个物种之间的循环相互作用，发现当生物多样性消失时，这五个物种中的两个可以共存，这与石头剪刀博弈不同。我们从两个方面研究了“手指”模型：直接模拟REs、偏微分方程和非线性偏微分方程。一个重要的发现是，循环竞争博弈中的物种数量对生物多样性的出现有影响。
对于问题三：在RPS对策中考虑了传染病传播的平均场动力学，然后将RPS对策模型与SIS传染病模型相结合，得到了平衡状态下三斑块系统的数值解。其次，我们还采用了种群模型。所有个体通过随机行走在斑块之间自由移动，以模拟流行病在未隔离状态下的传播。因此，我们在未感染者、感染者、恢复者之间找到了一个数值平衡，并得到了这个困境的最终解决方案。
问题四：在建立和分析前三个问题模型，特别是第三个问题的基础上，对COVID-19进行了较为详细的分析和解决，结合各国从疫情爆发到现在的现状，认真研究了各国对疫情的态度，最后结合我院向ICM写了一份关于各国防疫现状和全球防疫趋势的备忘录，主要从三个方面进行了详细的分析和介绍：全球防疫可实现的最低目标、目前必须解决的一些问题以及任何可能加速或阻碍实现目标的因素。
模型假设：
影响安全检查的因素很多。我们的模型是在一个相对理想的环境中建立和实现的。因此，我们做出如下假设：
（1） 假设生物环境相对封闭。
（2） 即使经过很长一段时间，玩家之间的互动也不会改变，也就是游戏规则不会改变。
（3） 在讨论玩家关系时，我们只考虑一次失败，假设没有其他情况，例如双向失败。
（4） 对于该流行病，假定人们可以恢复并获得终身免疫。
（5） 假定在疫情期间，无论自然出生和死亡人数如何，总人口保持不变。
（6） 我们假定所参考的数据是真实可靠的。
问题分析：
对于问题1
为了模拟玩家之间的关系，我们开发、分析和模拟了一个May-Leonard模型，以供参考。但是我们的目标是试图在多个玩家（多于三个玩家）之间找到一个平衡点，所以我们扩展了经过深入研究的三人石头剪刀（RPS）系统。然后我们发现，玩家的平衡不仅与玩家的数量有关，还与是否是偶数或奇数有关，这也有着重要的影响。所以我们把游戏者的数量分成偶数和奇数进行分析，用N=4和N=5表示。我们详细分析了4个参与者和5个参与者之间的关系，确定了空间均匀初始条件（即ode模型）的动力学，以及扩散和循环动力学（pde模型）驱动的不同参与者斑块的传播和相互作用。为了更好地适应背景，我们将玩家视为一个物种。
对于问题二
循环竞争博弈模型，特别是石头剪子模型，在探索空间生态系统中多物种共存问题中发挥着重要作用。我们提出了一个扩展的石头剪刀博弈模型来模拟五个物种之间的循环相互作用，发现当生物多样性消失时，这五个物种中的两个可以共存，这与石头剪刀博弈不同。当手指个数为5时，我们将新模型命名为五种群博弈，其中拇指、食指、中指、无名指和小指周期性地支配其后续种群，且由其先前种群支配。我们从两个方面研究了“手指”模型：直接模拟和非线性偏微分方程。一个重要的发现是，循环竞争博弈中的物种数量对生物多样性的出现有影响。在种群规模、繁殖率、选择率、迁移率等变量相同的情况下，岩-布-剪模型有利于维持生物多样性。研究还表明，流动性和繁殖率可以促进或危害生物多样性。
对于问题3
本文详细分析了N=4和N=5种群循环生态系统的确定性模型。我们发现：对于N=4的条件下共存（生物多样性）总是不稳定的。对于N=2n（N≠2）与N=4相同。在N=5的条件下，所有的生物选择态都是不稳定的。当N=2n+1（N≠3）与N=5时，生物多样性是最稳定的。下面我们还以COVID-19感染动力学为例，探讨基于这一结论的社会困境解决方案。
对于问题4
对于ICM的备忘录，我们主要描述了三个层次。一个切实可行的全球目标——可达到的最低防疫合作水平；当前全球防疫工作必须解决的问题；任何可能加速或阻碍实现这一目标的情况下，我们分析了各国的疫情现状，并分析了各自对防疫工作的态度进行了总结，最后回到全球防疫斗争本身得出结论。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
import numpy as npimport matplotlib.pyplot as pltpredlabel = 'Predator Count (Thousands)'preylabel = 'Prey Count (Thousands)'timelabel = 'Time'predcolor = 'red'preycolor = 'blue'class Lotka_Volterra(object): '''Sets up a simple Lotka_Volterra system'''def __init__(self, pdgrow,pddie,pygrow,pydie,tmax,timestep,prey_capacity=100.0,predator_capacity = 100.0):self.n = int(tmax/timestep) self.dt = timestep self.time = np.zeros(self.n) self.prey = np.zeros(self.n) self.predator = np.zeros(self.n) self.preygrow = pygrow self.preydie = pydie self.predgrow = pdgrow self.preddie = pddie self.prey_capacity = prey_capacity self.predator_capacity = predator_capacitydef set_initial_conditions(self,pdzero,pyzero, tzero=0.0): '''set initial conditions''' self.prey[0] = pyzero self.predator[0] = pdzero self.time[0] = tzero def plot_vs_time(self, filename='populations_vs_time.png', plot_capacity=False):fig1 = plt.figure() ax1 = fig1.add_subplot(111) ax2 = ax1.twinx() ax1.set_xlabel('Time', fontsize=22) ax1.set_ylabel(predlabel,fontsize=22, color=predcolor) ax1.tick_params('y', colors=predcolor) ax2.set_ylabel(preylabel,fontsize=22, color=preycolor) ax2.tick_params('y', colors='blue', color=preycolor) ax1.plot(self.time, self.predator, label='Predator', color=predcolor, linestyle='dashed')ax2.plot(self.time, self.prey, label = 'Prey', color = preycolor) if(plot_capacity): ax2.axhline(self.prey_capacity, label= 'Prey carrying capacity', color=preycolor, linestyle='dotted') plt.show() fig1.savefig(filename, dpi=300)def plot_predator_vs_prey(self, filename = 'predator_vs_prey.png'):fig1 = plt.figure() ax1 = fig1.add_subplot(111)ax1.set_xlabel(predlabel,fontsize=22) ax1.set_ylabel(preylabel,fontsize=22) ax1.plot(self.predator,self.prey, color='black') plt.show() fig1.savefig(filename,dpi=300)def plot_both_figures(self):fig1 = plt.figure()ax1 = fig1.add_subplot(211) ax2 = ax1.twinx() ax1.set_xlabel(timelabel) ax1.set_ylabel(predlabel, color=predcolor) ax2.set_ylabel(preylabel, color =preycolor) ax1.plot(self.time, self.predator, label='Predator', color=predcolor) ax2.plot(self.time, self.prey, label = 'Prey', color = preycolor) ax1.legend()ax3 = fig1.add_subplot(212)ax3.set_xlabel(predlabel) ax3.set_ylabel(preylabel) ax3.plot(self.predator,self.prey, color = 'black') plt.show()import numpy as npimport osfrom scipy import specialimport itertoolsimport timeimport warningswarnings.filterwarnings('error')class actor_critic: def __init__(self, dim_theta=3, d=47): self.dim_theta = dim_thetaself.theta = np.array([[1],[-1],[3]]) # good initialization # initialize weight vector (column) for value function approximation self.w = self.init_w(d) self.d = d self.tensor_phi = np.zeros([self.d, self.d, self.dim_theta])self.mat_alpha = np.zeros([self.d, self.d])def train_log(self, vector, filename): f = open(filename, 'a') vector.tofile(f, sep=',', format="%f") f.write("\n") f.close()gradient = self.calc_gradient_vectorized(P, pi) self.theta = self.theta - lr_actor * delta * gradient discount = discount * gamma pi = pi_next total_cost += cost list_cost.append(total_cost) if (episode % consecutive == 0): print("Theta\n", self.theta) print("pi\n", pi) cost_avg = sum(list_cost)/consecutive print("Average cost during previous %d episodes: " % consecutive, str(cost_avg)) list_cost = [] if write_file: self.train_log(self.theta, file_theta) self.train_log(pi, file_pi) self.train_log(np.array([cost_avg]), file_cost)if __name__ == "__main__": ac = actor_critic() t_start = time.time() ac.train(num_episodes=1, gamma=1, lr_critic=0.1, lr_actor=0.1, consecutive=1, write_file=0) t_end = time.time() print("Time elapsed", t_end - t_start)