2016年第五届数学建模国际赛小美赛
C题 对超级细菌的战争
原题再现：
最近有很多关于我们抗生素耐药性危机的讨论。进化出的能够抵抗抗生素的细菌每年杀死70万人，越来越强大的细菌正在世界各地传播。研究人员担心，我们将进入一个后抗生素时代，在这个时代里，我们被细菌感染，这些细菌可以击败药物提供的每一种药物。下周，联合国将召开一次高级别会议，协调全球打击这些无形敌人的斗争。
巴顿和其他志同道合的科学家们在60多年前就警告说，抗生素危机即将来临，尽管他们今天大多被遗忘了。他们是对的，但他们被忽视了。他们的失败为今天的新十字军提供了一些重要的教训。
然而，就在第二次世界大战结束一年后，青霉素的发现者警告说，它可能会变得毫无用处。1945年，亚历山大弗莱明获得诺贝尔奖，他指出了失败的原因。弗莱明说：“无知的人很容易给自己剂量不足，这是一种危险，通过将他的微生物暴露在非致命数量的药物中，使其产生耐药性。”
这并不是说其他科学家否认了进化论的真实性。这似乎并不重要。当时的科学家发现了更多的抗生素。如果细菌对青霉素产生抗药性，它们总是能够转换成另一种武器。
在一场持续了十多年的斗争中，改革者推动食品和药品管理局对抗生素进行更严格的监管。他们取得了成功，从市场上撤下了大量违法药品，并制定了更严格的规则，允许新药品进入市场。
与此同时，抗生素耐药性的真正威胁越来越明显。在20世纪60年代早期，科学家们发现，一旦一种细菌中进化出一种抗性基因，微生物就可以将其捐赠给其他细菌。微生物可以将这些捐赠的基因装载在一块DNA上，进一步加速耐药性的传播。
1966年，Look杂志上的一篇文章用近乎天启的修辞描述了抗生素耐药性的新观点。“细菌是不是赢得了与人类的战争？”标题写道。文章用一个新的绰号来提及这些细菌：超级细菌。
然而，事实证明，组织一场对抗抗生素耐药性的斗争比对抗无效或危险的药物要困难得多。早在20世纪50年代，世界卫生组织就组织了关于抗生素耐药性的会议，但最终失败了。参加会议的专家们陷入了关于如何衡量抵抗力以及考虑对公共健康的威胁程度的争论之中。
同时也很难弄清楚如何对抗抗生素耐药性。为了使劣药退出市场，改革者只有一个目标：FDA。为了减少良药的使用，改革者必须同时达到许多不同的目标：医生、医院管理人员、病人、政府、制药公司——甚至是那些开始给牲畜喂食抗生素以使牲畜长大的农民。
更糟糕的是，规范抗生素的运动让医生们很生气。他们说，“FDA是谁把这些药拿走的？我已经用了30年了，我的病人似乎越来越好。”
联合国已经要求你的团队，ICM-FDA帮助你更好地了解抗生素耐药性危机的相关因素。
任务1：建立一个模型，为抗生素的使用和各方的敏感性提供利益链。在建模过程中，可能需要考虑影响供需的因素的动态性质。
任务2：在完全竞争市场下预测危机的发展趋势。
任务3：设计一个便携式抗生素使用管理系统是超级细菌的当务之急。ICM-FDA已被要求参加一个政策策略会议，要求您的团队就您的模型编写一份报告并提出一套政策。
任务4：制定奖励政策，鼓励饲养者在不使用抗生素的情况下饲养牛。
整体求解过程概述(摘要)
最近有很多关于我们抗生素耐药性危机的讨论。超级细菌的出现使抗生素耐药危机成为亟待解决的问题。
针对课题一，通过供需动态因素分析，建立抗生素使用利益链，分析参与各方对抗生素使用的敏感性，并考虑抗生素对易感人群的影响。首先，通过供需链引出关于抗生素使用的利益链。对于敏感度分析，采用单因素敏感度分析方法对各相关方的敏感度进行分析，得出患者对抗生素使用的敏感度最高，政府对抗生素使用的敏感度最低的结论。
针对任务二，结合任务中的供需关系中的利益链模型，从供给因素、需求因素、细菌因素三个方面，分析抗生素耐药的相关因素。以两个指标共7项为参考数据，利用BP神经网络建立了抗生素耐药危机发展趋势预测模型，并对图像进行拟合，即对发展趋势进行拟合，从而得出结论：在完全竞争的市场下，抗生素耐药危机将向更为严重的方向发展。
针对任务三，结合任务一的利益链和任务二的危机预测，创建了便携式抗生素使用管理系统。首先确定管理系统的参与者是相关部门，然后创建系统用例模型，最终生成一个综合分类管理系统。最后，根据任务一和任务二，建立管理制度，提出合理的政策，遏制抗生素滥用。
针对任务四，主要根据前三个任务的研究结果，综合考虑供需关系、抗生素耐药危机及相关部门等因素，平衡各方利益，提出建立渔牧兽医服务站的建议，它将为农民饲养的牛提供各种服务，包括对牛的防疫和治疗以及监测抗生素的使用。禁止超过抗生素使用标准的黄牛进入市场，对于不使用抗生素饲养黄牛的农户，将获得市场价格3%的奖励，然后利用奖励矩阵建立政策的可行性分析模型，并通过进行详细分析，验证了政策的可行性。
模型假设：
（1） 对于同一疾病，我们总能找到一种仿制药代替抗生素，但没有明显的抗菌效果。
（2） 只考虑一般疾病，针对特定疾病予以清除。
（3） 牲畜中抗生素残留检测标准的存在。
（4） 排除抗生素生产难度等无关因素。
问题重述：
任务一的陈述与分析
建立一个模型，为抗生素的使用和各方的敏感性提供利益链。在建模过程中，应考虑影响供给和需求的因素的动态性质。对于第一个任务，首先要明确利益链的定义，利益链就是利益链上下各方的利益链。对于第一个问题，首先要明确利益链的定义，利益链是链条上人与下游人之间的利益关系。首先，通过对标题的分析得出一些关键角色：医生、医院管理者、政府、患者、制药公司和养殖户。首先对供应链进行分析，得到抗生素生产和使用过程中的供需关系，得到利益链模型，并对各参与方的利益链模型从供应侧（抗生素生产商）和需求侧（易感人群）两个方面进行敏感性分析。
任务二的陈述与分析
在完全竞争市场下预测危机的发展趋势。对于第二项任务，首先必须明确市场完全竞争的概念和危机，即在政府等外部因素缺位的情况下，市场完全由“看不见的手”控制抗生素耐药性危机。为此，总结了抗生素生产和使用过程的流程图，找出影响抗生素耐药性的六个因素，利用灰色神经网络系统建立危机预测模型，通过数据采集和归一化处理，利用模型预测发展趋势。
任务三的重述与分析
设计一个便携式抗生素使用管理系统是当务之急。我们应该就我们的模式写一份报告，并提出一套政策。对于任务三，抗生素的滥用已经成为一个日益严重的问题，为了避免超级细菌的产生，我们需要对抗生素的使用进行管理和监测。首先明确便携式管理系统的用户是相关部门，将系统的参与者确定为相关部门、抗生素医疗机构和抗生素农业机构，建立系统用例模型，定义每个参与者的功能，最后生成一个层次化的管理系统。针对不同类型的参与者和不同方面的措施，从医疗机构、农户、饲料生产企业、知识普及等方面提出了具体的政策建议。
任务4的重述和分析
我们被要求为农民制定一项奖励政策，鼓励他们饲养没有抗生素的牛。对于第四项任务，我们需要利用前三项任务的成果得出具体的激励政策，从市场供求、抗生素耐药危机的发展趋势、相关部门的具体措施等方面，结合农民群体的特殊性和自身利益，提出一些激励政策，使农民不使用抗生素喂养。然后利用报酬矩阵建立可行性验证模型，并利用模型分析得出激励政策是否可行。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
from scipy import stats
import numpy as np
from numpy.random import rand
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import chisquare
from scipy.optimize import curve_fit
pre_rst = [0.1, 0.12, 0.2, 0.15, 0.15, 0.17, 0.14, 0.14, 0.16,
0.4, 0.37, 0.38, 0.28, 0.28]
pre_use = [153, 136, 125, 117, 116, 111, 117, 121, 120, 113, 111, 105,
104, 102]
def cal_r(t, x0, x1, x2):r = (1/t)*math.log((1/x0-1/x1)/(1/x1-1/x2))return r
def cal_n(r, t, x0, x1):n = (1-math.e**(-r*t))/(1/x1-math.e**(-r*t)/x0)return n
def cal_rst(t, r, n):return n/(1+(n/0.1-1)*math.e**(-r*t))
r=cal_r(2.0, 0.14, 0.28, 0.37)
n=cal_n(r, 2.0, 0.14, 0.37)
pdt_rst = []
years = len(pre_rst)-2
for i in range(1,years,1):pdt_rst.append(cal_rst(i,r,n))
# pre_use = pre_use[:len(pre_use)-3]
pre_use = sorted(pre_use)
pdt_rst = sorted(pre_rst)
# z1 = np.polyfit(pre_use, pdt_rst, 3)#用 3 次多项式拟合
# p1 = np.poly1d(z1)
# 模拟医院，医生，病人的动作
def func(x,a,b):return a*np.log(b*x)
popt, pcov = curve_fit(func, pre_use, pdt_rst)
a=popt[0]#popt 里面是拟合系数，读者可以自己 help 其用法
b=popt[1]
print a,b
print 'func:',func(147,a,b)
print 'optimal:',func(120,a,b)
class Patient:# p_id 为病人的病号，drug 为当前用药,resistant 为抗药性，recover 为恢复状
况，cure_days 为已经治疗的天数def __init__(self, pid, drug, resistant, recover, cure_days):self.pid = pidself.drug = drugself.resistant = resistantself.recover = recoverself.cure_days = cure_days
RESISTANT = 0.5
T1 = 5
T2 = 10
DAY_PATIENTS = 30 # 每天 30 个病人
p_id = 0
sim_cure_days = []
sim_d1_sale = []
pre_use = [153, 136, 125, 117, 116, 111, 117, 121, 120, 113, 111, 105,
104, 102]
mean_use = np.mean(pre_use)
# 描述进入医院的病人
def ent_hpt(patients):# day_rst = stats.poisson.pmf(np.arange(DAY_PATIENTS), RESISTANT) # 每
天病人的抗药性分布day_rst = rand(DAY_PATIENTS)tmp = np.mean(day_rst)for i in range(len(day_rst)):day_rst[i] = day_rst[i]*RESISTANT/0.5day_ptn = map(lambda x: Patient(p_id, 1, x, 0.0, 0.0), day_rst)patients += day_ptnreturn patients
def cure(patients, rec_ptn, drug1_sale, drug2_sale, use, total_cure_days):tmp = []for p in patients:if p.drug == 1 and p.cure_days > 0:# cdp = use*((0.5*(1-p.recover/p.cure_days) +
0.5*p.cure_days/T1)**2)# cdp = (use*(1-p.recover/p.cure_days))**2cdp = 0.3*((1-use)*(1 - p.recover / p.cure_days))if rand() < cdp:p.drug = 2rec_ratio = 0if p.drug == 1:rec_ratio = (float)((1-p.resistant))/T1drug1_sale += 1if p.drug == 2:rec_ratio = 1.0/T2drug2_sale += 1p.recover += rec_ratiop.cure_days += 1if p.recover >= 1:tmp.append(p)rec_ptn.append(p)total_cure_days.append(p.cure_days)for p in tmp:patients.remove(p)return patients, rec_ptn, drug1_sale, drug2_sale, total_cure_days
def sim_hosp(use):all_patients = []rec_patients = []cure_days = []d1_sale = 0d2_sale = 0for i in range(30):
patients = ent_hpt(all_patients)while all_patients:patients, recover_patients, d1s, d2s, cure_days = cure(patients,
rec_patients, d1_sale, d2_sale, use, cure_days)d1_sale = d1sd1_count = 0for p in rec_patients:if p.drug == 1:d1_count += 1sim_cure_days.append(np.mean(cure_days))sim_d1_sale.append(d1s)# print '平均治疗时间: ', np.mean(total_cure_days)# print '抗生素销售: ', d1_sale# print '二线药物销售: ', d2_sale# print '坚持使用抗生素病人数: ', d1_count
# 市场自由发展
next_cure_days = []
next_d1_sale = []
max_use = np.mean(pre_use)
RESISTANT = func(max_use,a,b)
curr_use = max_use
pdt_rst = [0.1792628046926503, 0.2943181818181818, 0.43063362673704436,
0.558124635993011, 0.6535518026442507, 0.713738475784049, 0.7476984370881284,
0.7656659431653637, 0.7748492466183726, 0.7794600206940622,
0.7817543053571258]
for rst in pdt_rst:print rstRESISTANT = rstsim_hosp(curr_use/200)
print sim_cure_days
print sim_d1_sale
sim_uses = np.arange(80,200,1)
min_use = np.min(sim_uses)
max_use = np.max(sim_uses)
for use in sim_uses:RESISTANT = func(use,a,b)sim_hosp(use)print use,'==>',RESISTANT,'==>',sim_cure_days[len(sim_cure_days)-1]
print sim_cure_days
print sim_d1_sale
d1s_min = np.min(sim_d1_sale)
d1s_max = np.max(sim_d1_sale)
sim_d1_sale = map(lambda x: x*1.0/d1s_max,sim_d1_sale)
cd_min = np.min(sim_cure_days)
cd_max = np.max(sim_cure_days)
# sim_cure_days = map(lambda x: x*1.0/cd_max,sim_cure_days)
profits = []
for i in range(len(sim_cure_days)):profits.append(sim_cure_days[i]+10.0/sim_d1_sale[i])
print profits
uses = sim_uses
min_pft = np.min(profits)
max_pft = np.max(profits)
print max_pft,min_pft
plt.figure(figsize=(30,9),dpi=80)
plt.yticks(np.linspace(min_pft-1,max_pft+1,num=(int)(max_pftmin_pft)+1,endpoint=True))
plt.xticks(np.linspace(min_use,max_use,num=(max_usemin_use)/5+1,endpoint=True))
plt.ylim(min_pft-1,max_pft+1)
plt.xlim(min_use-0.2,max_use+0.2)
plt.ylabel('Profit',fontsize=15)
plt.xlabel('Drug use',fontsize=15)
plt.plot(uses,profits,label='Profit',marker='d',linewidth=2.5,markersize=6,
color='fuchsia')
#
plt.plot(X,fit_vals,label='Fit',marker='s',linewidth=2.5,markersize=10,colo
r='cornflowerblue')
plt.legend(loc='upper right')
plt.show()
plt.figure(figsize=(30,9),dpi=80)
plt.yticks(np.linspace(min_pft-1,max_pft+1,num=(int)(max_pftmin_pft)+1,endpoint=True))
plt.xticks(np.linspace(min_use,max_use,num=(max_use-min_use)/5+1,endpoint=True))
plt.ylim(min_pft-1,max_pft+1)
plt.xlim(min_use-0.2,max_use+0.2)
plt.ylabel('Profit',fontsize=15)
plt.xlabel('Drug use',fontsize=15)
plt.plot(uses,profits,label='Profit',marker='d',linewidth=2.5,markersize=6,
color='fuchsia')
#
plt.plot(X,fit_vals,label='Fit',marker='s',linewidth=2.5,markersize=10,colo
r='cornflowerblue')
plt.legend(loc='upper right')
plt.show()
import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import chisquare
from scipy.optimize import curve_fit
pre_rst = [0.1, 0.12, 0.2, 0.15, 0.15, 0.17, 0.14, 0.14, 0.16,
0.4, 0.37, 0.38, 0.28, 0.28]
pre_use = [153, 136, 125, 117, 116, 111, 117, 121, 120, 113, 111,
105, 104, 102]
def cal_r(t, x0, x1, x2):r = (1/t)*math.log((1/x0-1/x1)/(1/x1-1/x2))return r
def cal_n(r, t, x0, x1):n = (1-math.e**(-r*t))/(1/x1-math.e**(-r*t)/x0)return n
def cal_rst(t, r, n):return n/(1+(n/0.1-1)*math.e**(-r*t))
# pre_rs =
[cal_r(1,0.1,0.12,0.2),cal_r(1,0.14,0.16,0.4),cal_r(2,0.14,0.16,0.37)]
# pre_ns =[cal_n(pre_rs[0],1,0.1,0.2),cal_n(pre_rs[1],1,0.14,0.4),cal_n(pre_rs[
2],2,0.14,0.37)]
#
# print pre_ns
r=cal_r(2.0, 0.12, 0.15, 0.17)
n=cal_n(r, 2.0, 0.12, 0.17)
print 'r=',r
print 'n=',n
pdt_rst = []
years = len(pre_rst)-2
for i in range(1,years,1):pdt_rst.append(cal_rst(i,r,n))
next_rsts = []
for i in range(12,22,1):next_rsts.append(cal_rst(i,r,n))
print 'next_years:',next_rsts
data_names = range(2002,2013,1)
X = range(len(data_names))
print 'pdt_rst', pdt_rst
plt.figure(figsize=(30,9),dpi=80)
plt.yticks(np.linspace(0,1,num=11,endpoint=True))
plt.xticks(range(len(data_names)),data_names)
X = range(len(data_names))
plt.ylim(0,1.2)
plt.xlim(-0.2,len(data_names)+0.2)
plt.xlabel('Drug use',fontsize = 15)
plt.ylabel('Resistant',fontsize=15)
plt.plot(X,pdt_rst,label='Predict',marker='d',linewidth=2.5,markersiz
e=6,color='fuchsia')
plt.plot(X,pre_rst[3:],label='Real',marker='s',linewidth=2.5,markersi
ze=10,color='cornflowerblue')
plt.legend(loc='upper right')
plt.show()
# pre_use = pre_use[:len(pre_use)-3]
# pre_use = sorted(pre_use)
# # print pre_use
# pdt_rst = sorted(pre_rst[3:])
# def func(x,a,b):
# return a*np.log(b*x)
# popt, pcov = curve_fit(func, pre_use, pdt_rst)
# a=popt[0]
# b=popt[1]
# print a,b
plt.plot(X,pdt_rst,label='Real',marker='d',linewidth=2.5,markersize=6
,color='fuchsia')
#
plt.plot(X,fit_vals,label='Fit',marker='s',linewidth=2.5,markersize=1
0,color='cornflowerblue')
# plt.legend(loc='upper right')
# plt.show()
pdt_rst = [0.1792628046926503, 0.2943181818181818, 0.43063362673704436,
0.558124635993011, 0.6535518026442507, 0.713738475784049,
0.7476984370881284, 0.7656659431653637, 0.7748492466183726,
0.7794600206940622, 0.7817543053571258]
data_names = range(2013,2024,1)
plt.figure(figsize=(30,9),dpi=80)
plt.yticks(np.linspace(0,1,num=11,endpoint=True))
plt.xticks(range(len(data_names)),data_names)
X = range(len(data_names))
plt.ylim(0,1.2)
plt.xlim(-0.2,len(data_names)+0.2)
plt.xlabel('Year',fontsize = 15)
plt.ylabel('Resistant',fontsize=15)
plt.plot(X,pdt_rst,label='Prediction',marker='d',linewidth=2.5,marker
size=6,color='fuchsia')
#
plt.plot(X,fit_vals,label='Fit',marker='s',linewidth=2.5,markersize=1
0,color='cornflowerblue')
plt.legend(loc='upper right')
plt.show()
