全国大学生数学建模竞赛2004优秀论文:C、D题().pdf
C 题之一(全国一等奖)酒精在人体内的分布与排除优化模型桂林工学院,袁孟强,王哲,张莉指导教师:数模辅导组摘要:酒精进入机体后,随血液运输到各个器官和组织,不断的被吸收,分布,代谢,最终排除体外。为了研究酒精在体内吸收,分布和排除的动态过程,以及这些过程与人体反应的定量关系,本文建立了一个酒精在人体内的分布与排除优化模型,在药物动力学的一室模型的基础上,进行优化,改进,分别建立了酒精在人体内分布的房室模型 和房室模型 ,以及酒精在人体内的静态排除模型 和动态排除模型 ,导出模型的体液酒精浓度的状态函数,用常数交叉拟合方法,采用 VB编写程序,得到两个重要系数01k 和10k。根据此模型,计算的体液酒精浓度理论值与实验值十分相符,并很好地解释了给出的所有问题,得到一些有价值的结论。关键词:房室模型,排除模型,体液酒精浓度,动态和静态的转换酒精在人体内的分布与排除优化模型一、问题的重述国家质量监督检验检疫局2004 年 5 月 31 日发布了新的 车辆驾驶人员血液、呼气酒精含量阈值与检验国家标准,新标准规定,车辆驾驶人员血液中的酒精含量大于或等于20毫克百毫升,小于 80 毫克百毫升为饮酒驾车(原标准是小于100 毫克百毫升),血液中的酒精含量大于或等于80 毫克百毫升为醉酒驾车(原标准是大于或等于100 毫克百毫升)。大李在中午12 点喝了一瓶啤酒,下午6 点检查时符合新的驾车标准,紧接着他在吃晚饭时又喝了一瓶啤酒,为了保险起见他呆到凌晨2 点才驾车回家,又一次遭遇检查时却被定为饮酒驾车,这让他既懊恼又困惑,为什么喝同样多的酒,两次检查结果会不一样呢?参考下面给出的数据(或自己收集资料)建立饮酒后血液中酒精含量的数学模型,并讨论以下问题:1.对大李碰到的情况做出解释;2.在喝了3 瓶啤酒或者半斤低度白酒后多长时间内驾车就会违反上述标准,在以下情况下回答:1)酒是在很短时间内喝的;2)酒是在较长一段时间(比如2 小时)内喝的。3.怎样估计血液中的酒精含量在什么时间最高。4.根据你的模型论证:如果天天喝酒,是否还能开车?5.根据你做的模型并结合新的国家标准写一篇短文,给想喝一点酒的司机如何驾车提出忠告。参考数据1.人的体液占人的体重的65%至 70%,其中血液只占体重的7%左右;而药物(包括酒精)在血液中的含量与在体液中的含量大体是一样的。2.体重约 70kg 的某人在短时间内喝下2 瓶啤酒后,隔一定时间测量他的血液中酒精含量(毫克百毫升),得到数据如下(表 1):时间(小时)0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5 酒精含量30 68 75 82 82 77 68 68 58 51 50 41 时间(小时)6 7 8 9 10 11 12 13 14 15 16 酒精含量38 35 28 25 18 15 12 10 7 7 4 二、模型假设1、酒精的转移速率,及向体外排除的速率,与该室的血酒浓度成正比。2、酒精的转移速率,及向体外的排除速率,与时间有关,与空间(人体的各个部分)无关。3、中心室与体外有酒精交换,及酒精从体外进入中心室,最后又从中心室排出体外。与转移和排除的数量相比,酒精的吸收可以忽略。三、模型建立与求解房室模型(在短时间内喝下酒精量为0D)在短时间内喝下酒精量为0D,酒精进入胃,人体吸收酒精,然后排除出体外。吸收酒精的过程相当于酒精进入体液(中心室)的过程,全过程可以简化为下图:)(txkf0010)(110txkfout排除体外建模过程:0D短时间内进入胃的酒精;01k为胃室(吸收室)进入中心室的转移速率系数(由人体机能确定的常数);)(0tx是 t 时刻胃室的酒精;其微分方程为:0000100Dxtxktx())(1tx是 t 时刻进入中心室的酒精,其微分方程为:tVctxtftxktx1101101(2)酒精进入中心室的速率为:)(0010txkf(3)将方程(1)的解代入(3)得:tkekDtf010100(4)房室模型(在较长一段时间内喝酒)假设在较长的一段时间内喝下的酒是匀速进入胃室,则简化如下图:inf常数)(txkf0010)(110txkfout胃室)(0tx中心室)(1tx中心室)(1tx排除胃室)(0tx建模过程:inf为酒精进入胃的速率:tDfin0,t为喝酒时间。outf为酒精从中心室排除体内的速率0f 为酒精进入中心室的速率01k 为胃室进入中心室的转移速率(由人体机能确定的常数)10k 为是酒从中心室向外排除的速率系数。)(0tx是 t 时刻胃室的酒精,微分方程为:00)(00010 xftxktxin(5))(0010txkf(6))(1tx是 t 时刻进入中心室的酒精将方程(5)的解代入(6)得:tkinekftx011010(7)tkineftf0110(8)静态排除模型与房室模型 I 配套的静态酒精排除模型)(1tc中心室的血酒浓度;V人体体液量和人体血液量;酒精进入中心室的速率:tkekDtf010100)(1tx中心室的酒精量;微分方程为:tVctxtftxktx1101101(9)10k酒精从中心室向体外排除的速率系数(由人体机能确定的常数)由方程(9)得:Vtftcktc01101(10)对应的通解为:cdteVtfetctktk101001微分方程的解为:cdteVtfetctktk101001cdteVkDetkktk011010010cekkeVkDtkktk01101001100101tktkceekkVkD100101100101令0)0(1c得:10011kkc.tktkeekkVkDtc1001011001011.根据参考数据表 1,已知:短时间内进入胃的酒精0D,人体体液量 V 和一批实验数据(it,)(1itc)(231i)。用交叉常数拟合原理在VB环境中编写程序2,利用该程序算出两个重要系数01k 和10k,给出模型的状态函数.若初始值设为)0(00cc,则011001001kkkDVcctktktkececekkVkDtc10010100011001011动态排除模型 与房室模型配套建立动态酒精排除模型.)(1tc中心室的血酒浓度;V人体体液量和人体血液量;酒精进入中心室的速率:tkineftf0110)(1tx中心室的酒精量;微分方程为:tVctxtftxktx1101101(11)tVctx11(12)10k酒精从中心室向体外排除的速率系数(由人体机能确定的常数).由方程(11)得:Vtftcktc01101(13)对应的通解为:cdteVtfetctktk101001.将0)0(1c及tkineftf0110代入得:tktkinekkkekkkVftc100110011001101011111(240t)(14)四、酒过程的描述1、在短时间内喝下酒精量为0D用静态排除模型 描述:tktkeekkVKDtc1001011001011.2、在较长一段时间内喝酒0D用动态排除模型 描述喝酒的过程,用静态排除模型 描述酒后的过程.即先用状态函数tktkinekkkekkkVftc100110011001101011111描述喝酒的过程,然后用状态函数描述酒后的过程.3、天天喝下酒精量为0D用动态排除模型 描述第一天喝酒(喝一小时)的过程,用静态排除模型 描述酒后 23 小时内的的过程,用动态排除模型 描述第二天喝酒(喝一小时)的过程,用静态排除模型 描述第二天酒后 23 小时内的的过程,再用动态排除模型 描述第三天喝酒(喝一小时)的过程,用静态排除模型 描述第三天酒后 23 小时内的的过程,tktkinekkkekkkVftc100110011001101011111第一天tktkceekkVkD100101100101第二天tktkeekkVKDtc1001011001011tktkinekkkekkkVftc100110011001101011111第二天tktkceekkVkD100101100101第三天,五、参数的选择一瓶啤酒的酒精量4:)24192%8.4(8.0)(6300mgmlD()比重.人体体液:)()(百毫升比重体重4301.1%5.67)(701kgV速率系数5:2.083.11001kk两小时内慢慢喝下两瓶啤酒的输入率:hmgfin241922224192六、拟合效果.中心室中酒精含量(毫克百毫升)浓度状态函数)(1tc拟合效果图效果图显示拟合程度极高,说明参数的选择与客观情况相符合.七、问题分析问题 1 的分析:大李在中午喝了一瓶啤酒下午6 点检查时,利用静态排除模型:.;1110001100101011001tVctxeDtxceekkVkDtctktktk;.6.5483)6(1.81801)6(/20/0234.19)6(101毫克毫克;百毫升(标准);毫克百毫升毫克xxc由浓度状态函数在6t时的浓度:200234.19)6(1c,可知大李此时符合新的驾车标准.紧接着他在吃晚饭时又喝了一瓶啤酒,为了保险起见他呆到凌晨2 点才驾车回家.检测时距晚饭喝的那瓶啤酒已过了八小时,其胃中酒精及体液中的酒精含量分别为:毫克);毫克);百毫升(标准);毫克(33.89888(6.54838x/209031.208101xc由浓度状态函数百毫升(标准)毫克/209031.20)8(1c可知,检查时被定为饮酒驾车.之所以被定为饮酒驾车,关键是此时方程Vtftcktc01101中的初始条件0234.1901c而不是第一次喝酒的001c.问题 2 的分析:在很短时间内喝了3 瓶啤酒,问多长时间内驾车就会违反新的驾车标准.根据静态排除模型:3 瓶啤酒的酒精量:毫克)(725760D;浓度状态函数为tktkeekkVkDtc1001011001011.代入数据得:tteetc83.120.0149076.189.(15)若违反新的驾车标准,则:百毫升毫克201tc,根据人体体液酒精浓度曲线图可知,喝酒后约在t=0.07 小时(42 秒)与 t=11.24 小时之内驾车会违反新标准。将 t=0.07,t=11.24 分别代入方程(15)检验得出:200122.2024.11201490.2007.011cc通过验证,证明观测值基本接近实际值。问题 3 的分析:在较长一段时间(比如 2 小时)内喝 3 瓶啤酒,多长时间内驾车就会违反新标准.假设匀速喝酒,则tDfin0此过程分两阶段:(1)在喝酒过程中,多少时间后驾车会违反新标准。(2)喝酒之后,多长时间内驾车会违反新标准。根据动态模型:tktkinekkkekkkVftc100110011001101011111tktkekkkekkktVD100110010101101001111(16)(1)由人体体液酒精浓度曲线观测出:当 t(20t)大于 0.62 小时(约37.2分 钟)时,体 液 酒 精 浓 度 大 于201tc毫 克,把 已 知 酒 精 数 据)(毫克725760D,2.0,83.11001kk,),(百毫升430V代入方程(15)检验得:20125.2062.01c通过数据验证,证明观测值基本接近实际值。(2)由阶段(1)的已知数据算出,当 t=2 小时(即停止喝酒时)人体酒精为:百毫升毫克7371.10521c此时,根据胃里剩余的酒精方程:tkeDtx0100得出:)(38638)2(0mgx以 2 小时作为零时刻,设状态初始值为0c,则取210cc,此时静态排除模型为:tktkekkkDVcekkVkDtc10010110010101100101121根据此模型,可得曲线:可观测出:当201tc时,t(t2)的值约在 11.67 附近。再将 t=11.67 代入tc1进行检测得:0228.2067.111c(毫克/百毫升)通过数据验证,证明观测值基本接近实际值。结论:若连续 2 小时均匀喝下 3 瓶啤酒,则在 0.62 与 13.67 小时之内驾车会违反新标准。问题 4 的分析:估算出血液中酒精含量何时达到峰值。最高值的估计要分两种情况讨论:一段时间喝完)据动态排除模型(很长喝完的情况)据静态排除模型(一次1、据静态排除模型(一次喝完的情况)酒精进入机体之前,胃里的酒精浓度为0,当喝酒后,胃里的酒精浓度逐渐增大,因为存在着酒精进入中心室的过程,根据人体机理可知,outinff。当胃里的酒精浓度等于中心室的酒精浓度时,outinff,此时血液中酒精含量达到最高值。由静态排除模型得:tktkeekkVkDtc1001011001011若要使tc1的值达到最高,则要求01tc即:01001100101100101tktkekekkkVkDtc代入已知数据得:tteDeDtc2.0083.1019.700366.043063.183.183.1)(算出:t=1.36 小时即:喝酒 1.36小时后,血液中酒精含量最高。(1)根据动态排除模型(酒是在很长一段时间内喝完)喝酒时,由于酒量在不断的增加,根据人体机理可知,胃中的酒精浓度始终是大于中心室的酒精浓度,即outinff,易知血液中的酒精浓度tc1是一个递增的函数。那么在这阶段中,血液的酒精含量的最高值为停止喝酒时血液中的酒精含量。运用动态排除模型:tktkinekkkekkkVftc100110011001101011111可计算出最值。喝完酒后,胃中的酒精浓度仍然大于中心室的酒精浓度,此时,outinff.随着时间的推移,胃中的酒精浓度在不断的减少,而血液的酒精浓度不断增加,这就使得当outinff时,血液的酒精浓度达到最高值(即01tc),若代入具体数据就可算出最大值。问题 5 的分析:根据我们建立的模型论证:如果天天喝酒,是否还能开车?天天喝一瓶酒的数学模型见图一个人天天喝酒,如果每天只喝一瓶啤酒第二天即12 个小时后体内的血酒浓度tc1远远小于 20 毫克/百毫升所以仍能开车;天天喝两瓶酒的数学模型见图:一个人天天喝酒,如果每天只喝两瓶啤酒第二天即12 个小时后体内的血酒浓度tc1远远小于 20 毫克/百毫升所以仍能开车;天天喝三瓶酒的数学模型见图:如果每天喝三瓶啤酒,12 小时后体内的血酒浓度为18.920 毫克/百毫升,比新的国家标准低但很接近,此时仍可以开车却比较危险附表(1)如下:一瓶啤 酒的酒精 含量24192(mg)两瓶啤酒的酒精含量48384(mg)三 瓶啤酒的 酒精含量72576(mg)四 瓶啤酒的 酒精含量96768(mg)五 瓶 啤酒 的 酒精含量120960(mg)12 5.7301 12.6 18.9 25.2 28.604 24 0.5198 1.396 1.5595 2.0793 2.5991 注:据资料表明人体酒精浓度达到200(毫克/百毫升)300(毫克/百毫升)会使人昏迷,一般人喝五瓶啤酒(含酒精量 120960毫克)后,(据静态排出模型)在一小时左右浓度达到213.6729 毫克/百毫升已使人昏迷,因而研究五瓶啤酒以上的酒精量已没有意义。结论:每天喝三瓶以下的啤酒第二天不违反新的国家标准仍可以开车,而喝四瓶到五瓶就不能再开车了,五瓶以上大多数人都会昏迷摄入的酒精量体内的血酒浓度tc1时间(t)问题 6 的分析:通过以上论述,特给爱喝点酒的司机提个醒,希望他们能平安。科学研究表明,人体内每百毫升血液中酒精含量达到20 毫克,会出现头晕脑胀、兴奋健谈;达 60 80 毫克时,感情容易冲动,步态不稳;达120160毫克时,神经进入抑制状态,开始昏睡;达200400 毫克时,意识朦胧,呈木僵状态,如果达 400500 毫克时,就可能导致脑损伤和呼吸麻痹,从而死亡。根据状态函数与静态排除模型的分析方法可得出:喝完一瓶啤酒经 5.74 小时后,血液的酒精浓度仍为20.0383,超过新的驾车标准,为了安全,建议司机饮酒八小时后驾车。.喝完两瓶啤酒经 9.2 小时后,血液的酒精浓度仍为20.063,饮酒的司机在 11小时内不宜驾车。喝完四瓶啤酒经 12.66 小时后,血液的酒精浓度仍为20.086,饮酒的司机在14 个小时内不宜驾车。综上所述,为安全起见,对于爱饮酒的司机,每天的饮酒量不宜超过两瓶啤酒。若饮酒量超过四瓶,第二天司机不宜驾车。八、模型评价1、本模型绘出状态函数后,极好的解决了喝酒过程的全程数学描述,定量的解决了所提出的问题。2、本模型微分方程的待定系数求解方案没有利用现成的方法。而是自行编程解决,给此类模型提供了一个利用的程序。3、本模型通俗易懂,为大学生所接受,但反复运用状态函数,很有特色,尤以初值问题多次循环反复,妙用至极。4、有了本模型,人们对喝酒过程有一个较为精确的定量认识。从事有一定危险工作的人们应该如何正确的喝酒的依据。5、本模型可推广到有毒物质在人体的分布与排除。这对在有毒物场所工作的人起到一定的提醒作用。参考文献:1、ISBN 7-04-004505-2 姜启源,数学模型(第二版),北京:高等教育出版社,1993年 8 月2、ISBN 7-04-011943-9 徐全智 杨晋浩,数学建模,北京:高等教育出版社,2003年七月输入,D,V,M,t,N,1001,KKKKi=1 011000111001()()()IIKtKtiD Kc teeV KKiN i=i+1 1()()()iiic tM tM ti=1 开始iN i=i+1 输出结果10011,()iKKc t101010KKK100K010101KKK010K1010KKK1010010122KKKKKKKK10011001KKKKKK附件常数拟合流程图Y N N N Y Y Y N N 说明:是精度,D 是酒精量,V 是人体体液量,M 是实验的酒精含量数据,t 是实验的时间数据,1001,KKKK是代定常数初值,是步长值。源程序Dim T(23)As Single 实验值,时间(小时)Dim M(23)As Integer 酒精含量Private Sub Form_Load()T(1)=0.25 T(2)=0.5 T(3)=0.75 T(4)=1 T(5)=1.5 T(6)=2 T(7)=2.5 T(8)=3 T(9)=3.5 T(10)=4 T(11)=4.5 T(12)=5 T(13)=6 T(14)=7 T(15)=8 T(16)=9 T(17)=10 T(18)=11 T(19)=12 T(20)=13 T(21)=14 T(22)=15 T(23)=16 M(1)=30 M(2)=68 M(3)=75 M(4)=82 M(5)=82 M(6)=77 M(7)=68 M(8)=68 M(9)=58 M(10)=51 M(11)=50 M(12)=41 M(13)=38 M(14)=35 M(15)=28 M(16)=25 M(17)=18 M(18)=15 M(19)=12 M(20)=10 M(21)=7 M(22)=7 M(23)=4 End Sub Public Function C1(ByVal D0 As Single,ByVal V As Double,ByVal K01 As Single,ByVal K10 As Single,ByVal T As Single)As Single C1=D0*K01*(Exp(-K01*T)-Exp(-K10*T)/V/(K10-K01)End Function Private Sub cmdOK_Click()Const N=23 Const nStep=0.01 如果设得太小,计算时间会很长。当然,精确度会更高。Dim i As Integer Dim MaxK01 As Single Dim MaxK10 As Single Dim initK01 As Single Dim D0 As Single Dim V As Double Dim K01 As Single Dim K10 As Single Dim Ei As Single Dim IsPass As Boolean Dim nTimer As Single txtV.Text=Trim(txtV.Text)If Not IsNumeric(txtV.Text)Then MsgBox 请输入数值。,vbOKOnly+vbInformation,数据拟合 txtV.SetFocus Exit Sub End If If txtV.Text=0 Then MsgBox 请输入大于0 的数值。,vbOKOnly+vbInformation,数据拟合 txtV.SetFocus Exit Sub End If txtD0.Text=Trim(txtD0.Text)If Not IsNumeric(txtD0.Text)Then MsgBox 请输入数值。,vbOKOnly+vbInformation,数据拟合 txtD0.SetFocus Exit Sub End If If txtD0.Text=0 Then MsgBox 请输入大于0 的数值。,vbOKOnly+vbInformation,数据拟合 txtD0.SetFocus Exit Sub End If txtK01.Text=Trim(txtK01.Text)If Not IsNumeric(txtK01.Text)Then MsgBox 请输入数值。,vbOKOnly+vbInformation,数据拟合 txtK01.SetFocus Exit Sub End If If txtK01.Text=0 Then MsgBox 请输入大于0 的数值。,vbOKOnly+vbInformation,数据拟合 txtK01.SetFocus Exit Sub End If txtK10.Text=Trim(txtK10.Text)If Not IsNumeric(txtK10.Text)Then MsgBox 请输入数值。,vbOKOnly+vbInformation,数据拟合 txtK10.SetFocus Exit Sub End If If txtK10.Text=0 Then MsgBox 请输入大于0 的数值。,vbOKOnly+vbInformation,数据拟合 txtK10.SetFocus Exit Sub End If txtE.Text=Trim(txtE.Text)If Not IsNumeric(txtE.Text)Then MsgBox 请输入数值。,vbOKOnly+vbInformation,数据拟合 txtE.SetFocus Exit Sub End If If txtE.Text=0 Then MsgBox 请输入大于0 的数值。,vbOKOnly+vbInformation,数据拟合 txtE.SetFocus Exit Sub End If txtTime.Text=Trim(txtTime.Text)If Not IsNumeric(txtTime.Text)Then MsgBox 请输入数值。,vbOKOnly+vbInformation,数据拟合 txtTime.SetFocus Exit Sub End If If txtTime.Text=0 Then MsgBox 请输入大于0 的数值。,vbOKOnly+vbInformation,数据拟合 txtTime.SetFocus Exit Sub End If On Error GoTo Err1 txtC1.Text=txtC1.Refresh txtK01v.Text=txtK01v.Refresh txtK10v.Text=txtK10v.Refresh D0=txtD0.Text V=txtV.Text MaxK01=txtK01.Text MaxK10=txtK10.Text Ei=txtE.Text initK01=MaxK01 K01=MaxK01 K10=MaxK10 nTimer=Timer Me.MousePointer=11 While Timer-nTimer 0 K01=MaxK01 While K01 0 IsPass=True For i=1 To N If K01=K10 Then IsPass=False Exit For End If If Abs(C1(D0,V,K01,K10,T(i)-M(i)/M(i)=Ei Then IsPass=False Exit For End If Next i If IsPass=True Then GoTo Pass:Else K01=K01-nStep End If Wend K10=K10-nStep Wend MaxK01=MaxK01+1 MaxK10=MaxK10+1 K01=MaxK01 K10=MaxK10 End If While K10 MaxK10-1 K01=MaxK01 While K01 0 IsPass=True For i=1 To N If K01=K10 Then IsPass=False Exit For End If If Abs(C1(D0,V,K01,K10,T(i)-M(i)/M(i)=Ei Then IsPass=False Exit For End If Next i If IsPass Then GoTo Pass:Else K01=K01-nStep End If Wend K10=K10-nStep Wend While K10 0 K01=MaxK01 While K01 MaxK01-1 IsPass=True For i=1 To N If K01=K10 Then IsPass=False Exit For End If If Abs(C1(D0,V,K01,K10,T(i)-M(i)/M(i)=Ei Then IsPass=False Exit For End If Next i If IsPass Then GoTo Pass:Else K01=K01-nStep End If Wend K10=K10-nStep Wend If K10=d-p)cj(i)=4-j;else j=j+1;cj(i)=4-j;d=d-p;end end%a=4 4 3 3;4 3 4 2;3 4 1 2;4 3 3 3;3 4 3 2;3 1 4 3;4 3 2 3;3 4 4 2;3 3 4 3;1 3 4 2;1 2 3 4;4 3 2 4;3 2 1 4;1 3 4 3;4 3 2 3;3 4 3 2;a=a,cj;d(1:7)=0;b=3 4 2 4;4 3 3 2;4 3 3 2;2 2 4 4;2 2 4 4;2 3 3 4;2 3 3 4;b=b,d;%下面的程序求应聘者与部门要求的能力的相对差距矩阵t1(1:7,1:5)=0;for i=1:16 for j=1:7 t1(j,:)=a(i,:)-b(j,:);t(i-1)*7+j,:)=t1(j,:);end end a=t;%tt1(1:112)=0;tt2(1:112)=0;for j=1:112 for i=1:5 if i=0 tt1(j)=tt1(j)+landa2*a(j,i);else tt2(j)=tt2(j)+landa2*abs(a(j,i);end else if a(j,i)=0 tt1(j)=tt1(j)+landa1*a(j,i);else tt2(j)=tt2(j)+landa1*abs(a(j,i);end end end end d(1:16,1:7)=0;for i=1:16 for j=1:7 d(i,j)=tt1(i-1)*7+j);end end d1=d;d(1:16,1:7)=0;for i=1:16 for j=1:7 d(i,j)=tt2(i-1)*7+j);end end d2=d;d=d1-d2%上面的程序求应聘者与部门要求的能力的差距矩阵!下面是 Lingo 程序,求解招聘方案的解。!d2004.lg4 MODEL:SETS:w/1.16/:c;v/1.7/:d;