1、大学数学实验大学数学实验Mathematical Experiments 实验实验14 14 数学建模与数学实验数学建模与数学实验(综合综合/补充补充)案例：隧道阀门够结实吗？案例：隧道阀门够结实吗？阀门是否足够结实？阀门是否足够结实？阀阀门门爆破点爆破点湖湖水水山体山体hL问题初步分析问题初步分析爆破后隧道内湖水和空气的相互作用过程爆破后隧道内湖水和空气的相互作用过程0 x阀门L湖水加速进入隧道湖水加速进入隧道,然后减速直到停止然后减速直到停止可大致认为隧道内湖水和空气有清晰的接触界面可大致认为隧道内湖水和空气有清晰的接触界面隧道内的空气完全处于被不断压缩的过程中隧道内的空气完全处于被不断压
2、缩的过程中隧道内空气对隧道内湖水的阻力越来越大隧道内空气对隧道内湖水的阻力越来越大隧道内的湖水与隧道内的空气紊乱地混合隧道内的湖水与隧道内的空气紊乱地混合隧道内湖水和空气不再有一个清晰的接触界面隧道内湖水和空气不再有一个清晰的接触界面直到最后隧道内的空气以气泡形式从爆破口排出直到最后隧道内的空气以气泡形式从爆破口排出受力分析受力分析0F3xF1F1F2阀门L受力分析受力分析)()(000LAPVPVPAxLPtttxLLPPt00F3xF1F1F2阀门L受力分析受力分析0F3xF1F1F2阀门LgvRlhf242长长l、速度、速度v的流体两端压强之差（的流体两端压强之差（水头损失水头损失,mm
3、水柱水柱）gxDxgvRlgtP224)(22R为水力半径为水力半径(对圆管满流，对圆管满流，R=D/4)22322)()(xxDAxDxAtPAtF隧道内湖水的基本运动方程隧道内湖水的基本运动方程 热力学第一定律（能量转化和守恒定律）：热力学第一定律（能量转化和守恒定律）：当运动形式发生变化时，能量也从一种形式转化为当运动形式发生变化时，能量也从一种形式转化为另一种形式，从一个系统传递给另一个系统；在转另一种形式，从一个系统传递给另一个系统；在转化和传递中总能量始终不变。化和传递中总能量始终不变。假设隧道内湖水的内能没有变化，也没有热交换发生假设隧道内湖水的内能没有变化，也没有热交换发生 d
4、tdWdtdE隧道内湖水的质量为隧道内湖水的质量为)(tAxm222121xAxmvExxxxAdtdE221隧道内湖水的基本运动方程隧道内湖水的基本运动方程 dtdWdtdExxxxAdtdE221200312)(xxDAxLLAPAPghxFdtdxdtdWii2200212xxxDxLLPPghxx),(xxG0)0()0(),(xxxxGxxODE初值问题初值问题模型求解模型求解 2200212xxxDxLLPPghxx),(xxG0)0()0(),(xxxxGxx参数：重力加速度参数：重力加速度g=9.807m/s2，水的密度，水的密度=1kg/m3，P0=1atm(1大气压大气压)
5、=101325Pa，空气比热比，空气比热比=1.4;据水力学知识，沿程损失系数通常不超过据水力学知识，沿程损失系数通常不超过0.1,取取0.1;取取D=2m,h=100m，L=1000m0)0()0(2212122122100221xxxDxxxLLPPghxxx模型求解模型求解 0)0()0(),(xxxxGxxMatlab程序程序0)0()0(2212122122100221xxxDxxxLLPPghxxx显示显示“Warning:Divide by zero.”0)0()0(2)sgn(2121222122100221xxxDxxxxLLPPghxxx232sgn)(xxDAxtF仍然显
6、示仍然显示“Warning:Divide by zero.”模型求解模型求解0)0()0(),(xxxxGxx模型无解？模型无解？220021sgn2xxxxDxLLPPgh需要进行修改！需要进行修改！保持初始条件，修改基本方程保持初始条件，修改基本方程 保持基本方程，修改初始条件保持基本方程，修改初始条件也许存在也许存在x(t)0 满足方程？满足方程？考虑到考虑到x(t)及其导数在及其导数在t=0时刻的连续性，时刻的连续性，模型无解模型无解!修改基本方程修改基本方程 想象有一些爆破点附近的湖水被掀开，远离了隧想象有一些爆破点附近的湖水被掀开，远离了隧道的爆破点。记隧道外的额外长度为道的爆破点
7、。记隧道外的额外长度为d 0dx阀门阀门L0)0()0(),()(xxxxGxdxxdxxxAdtdE)(212)(dxAm222121xAxmvEtd=Dd=0.1Dd=0.01DxPxPxP001.000001.000001.00001.000032.87491.047935.48301.051935.86191.05252.000060.45281.091262.48791.094562.78431.09503.000083.48451.129885.16971.132785.41901.1332107.0000 833.412712.2942833.396912.2926833.396
8、512.2926108.0000 833.639512.3177833.610912.3147833.608412.3145109.0000 833.667112.3206833.625812.3163833.621312.3158110.0000833.496212.3029833.442212.2973833.435712.2966111.0000833.130412.2651833.064212.2583833.055712.2574112.0000832.577612.2085832.499912.2005832.489612.1995计算结果计算结果修改初始条件修改初始条件 可以想象
9、有一些爆破点附近的湖水已经进入了隧道可以想象有一些爆破点附近的湖水已经进入了隧道 记隧道内的额外长度为记隧道内的额外长度为d，x(0)=d 水如果从深度为水如果从深度为h的容器中自然流出，速度是的容器中自然流出，速度是 ghxdxxxGxx2)0(,)0(),(0Ldx阀门gh2计算结果计算结果td=Dd=0.1Dd=0.01DxPxPxP02.00001.00280.20001.00030.20001.00031.000037.37751.054836.05841.052843.27401.06392.000064.01731.097062.94211.095384.15111.13103.
10、000086.49211.135085.55521.1334123.09971.2019106.0000833.011112.2529832.988712.2506107.0000833.411612.2941833.398012.2927108.0000833.613512.3150833.608712.3145109.0000833.616212.3153833.620312.3157110.0000833.419812.2950833.432812.2963111.0000833.024812.2543833.046512.2565参数参数分析分析0.0550029.0000 435.9
11、43717.75710.05100079.0000 848.358014.23390.052000217.00001667.282612.31790.1050039.0000424.145014.01460.101000109.0000833.625712.31630.102000303.00001650.388511.4927LtmaxxmaxPmax0.0550038.0000 389.25848.25130.051000103.0000 751.17347.01040.052000282.00001469.99756.41880.1050051.0000375.57397.00940.1
12、01000141.0000735.00126.41890.102000392.00001452.16046.12810.0550021.0000 466.113443.30310.05100057.0000 913.533830.79040.052000158.00001801.394025.36600.1050029.0000456.526330.55210.10100079.0000900.698425.36650.102000221.00001785.849622.8263h=50h=100h=200Pmax减少h减少L增加 增加D可类似分析可类似分析yLghLghyghdtdsdsyd
13、ghdtxdx2222.2222无量纲化无量纲化Pmax减少h减少L增加 增加ghLtLxcc2/tLghttsLxxxycc2,sghLtLyx2,yghdsdyghdtdxx22.0)0()0(sgn11112220yyyyyyDLyghPyLdy0)0()0(),()(xxxxGxdx这说明影响这说明影响 y 的参数是的参数是ghPDL/,/0无量纲化无量纲化1111000yghPxLLghPghPPtPt与与 y 的关系的关系)/,/(0DLghP在坐标系在坐标系 O下，下，下，可画出下，可画出Pmax的等值线的等值线(单位：大气压单位：大气压)DL/ghP/0120140200Pma
14、x(atm)进一步研究进一步研究n推广推广 隧道不是水平的，怎么办？隧道不是水平的，怎么办？隧道不是直的，而是一条任意（下降的）曲线，隧道不是直的，而是一条任意（下降的）曲线，怎么办？怎么办？n 理论上分析基本模型与修正模型的性质理论上分析基本模型与修正模型的性质(解的解的存在性与唯一性，近似处理的误差，存在性与唯一性，近似处理的误差，etc.)2000年全国大学生数学建模竞赛年全国大学生数学建模竞赛B题题 (CUMCM-2000B)CUMCM-2000B)钢管订购和运输钢管订购和运输CUMCM-2000B CUMCM-2000B 钢管订购和运输钢管订购和运输由钢管厂订购钢管，经铁由钢管厂订购
15、钢管，经铁路、公路运输，铺设一条路、公路运输，铺设一条钢管管道钢管管道1521AAAA1325801010312012427010881070627030202030450104301750606194205201680480300220210420500600306195202720690520170690462160320160110290115011001200A2A3A4A5A6A7A8A9A10A11A12A13A14A15S1S2S3S4S5S6S7管道铁路公路S1S7 钢管厂火车站450 里程（km)（沿管道建有公路）钢 厂 i1234567产 量 上 限 is800800100
16、02000200020003000销 价ip（万 元）160155155160155150160钢厂的产量和销价（1单位钢管=1km管道钢管）钢厂产量的下限：500单位钢管里程(km)300301350351400401450451500运价(万元)2023262932里程(km)5016006017007018008019009011000运价(万元)37445055601单位钢管的铁路运价1000km以上每增加1至100km运价增加5万元1单位钢管的公路运价：0.1万元/km（不足整公里部分按整公里计）601=300+30144 20+23?（1）制定钢管的订购和运输计划，使总费用最小）制
17、定钢管的订购和运输计划，使总费用最小.（2）分析对购运计划和总费用影响：哪个钢厂钢管销价的）分析对购运计划和总费用影响：哪个钢厂钢管销价的变化影响最大；哪个钢厂钢管产量上限的变化影响最大？变化影响最大；哪个钢厂钢管产量上限的变化影响最大？A1325801010312012427010881070627030202030450104301750606194205201680480300220210420500600306195202720690520170690462160320160110290115011001200A2A3A4A5A6A7A8A9A10A11A12A13A14A15S1S2
18、S3S4S5S6S7A16130A17A18A19A20A21190260100（3）讨论管道为树形图的情形）讨论管道为树形图的情形问题问题1的基本模型和解法的基本模型和解法总费用最小的优化问题总费用：订购，运输（由各厂Si经铁路、公路至各点Aj,i=1,7;j=1,15），铺设管道Aj Aj+1（j=1,14)由Si至Aj的最小购运费用路线及最小费用cij 由Si至Aj的最优运量xij由Aj向Aj Aj-1段铺设的长度yj及向Aj Aj+1段铺设的长度zj最优购运计划最优购运计划约束约束条件条件钢厂产量约束：上限和下限（如果生产的话）运量约束：xij对i求和等于zj 加yj；zj与 yj+1
19、之和等于Aj Aj+1段的长度ljyj zjAj基本模型基本模型由Aj向Aj Aj-1段铺设的运量为 1+yj=yj(yj+1)/2由Aj向Aj Aj+1段铺设的运量为 1+zj=zj(zj+1)/2)6(0,0)5(15,2,1,7,2,10,0,0)4(14,2,1)3(15,2,1)2(7,2,1,5000.)1()1()1(21.0min15117115171151151zyjiyzxjlyzjyzxisxtsyyzzxcjjijjjjjjiijijijijjjjjjijij二次规划求解步骤求解步骤1）求由Si至Aj的最小购运费用路线及最小费用cij 难点：公路运费是里程的线性函数，而
20、铁路运费是里程的分段阶跃函数，故总运费不具可加性。因而计算最短路常用的Dijkstra算法、Floyd算法失效。7010881070627030202030300220210420500170690462160320160110290A10A11A12A13A14A15S4S5S6S7需要对铁路网和公路网进行预处理，才能使用常用需要对铁路网和公路网进行预处理，才能使用常用算法，得到最小购运费用路线。算法，得到最小购运费用路线。-至少求至少求3次最短路次最短路如S7至A10的最小费用路线先铁路1130km，再公路70km,运费为77（万元）先公路（经A15）40km,再铁路1100km，再公路7
21、0km,运费为76（万元）的的处处理理约约束束条条件件)7,2,1(,5000)2151 isxijij问问题题求求解解。，分分解解为为上上述述形形式式的的子子的的那那些些求求解解，再再对对解解中中满满足足先先松松弛弛为为ixisxbjijijij5000)7,1(0)151151 个个子子问问题题共共和和分分解解为为71511512)7,1(5000)isxxaijijjij实际上只有S4和S7需要分解成子问题求解每个子问题是标准的二次规划，决策变量为xij,yj,zj,不超过135个。0,0,0.)(05.0min1511711512151271151zyzyxlyzzyxsxtsyyzz
22、xcjjijjjjjijijjiijjjjjjijijijfi表示钢厂表示钢厂i是否使用；是否使用；xij是从钢厂是从钢厂i运到节点运到节点j的钢管量的钢管量yj是从节点是从节点j向左铺设的钢管量；向左铺设的钢管量；zj是向右铺设的钢管量是向右铺设的钢管量 c)比较好的方法：引入比较好的方法：引入0-10-1变量变量.7,.,1,1,0,0.14,.,1.15,.,1,.7,.,1,500.)1()1(21.0)(151171151151,ifzyjbzyjzyxifSxftszzyyxcpMinijjjjjiijiijijijjjjjjiijiji LINDO/LINGO得到的结果比得到的结
23、果比matlab得到的好得到的好cumcm2000b.lg4yj zjj问题问题1的其它模型和解法的其它模型和解法1）运输问题的0-1规划模型将全长5171km的管道按公里分段，共5171个需求点，钢厂为7个供应点，构成如下的运输问题5171,1,7,1,1,05171,1,17,1,500,0.min71517117151711 jixjxisxtsxcijiijijijijijijcij为从供应点i到需求点j的最小购运费xij=1表示从点i到点j购运1单位钢管求解时要针对规模问题寻求改进算法2）最小费用网络流模型）最小费用网络流模型SourceS1S2S7A1A2A15P11P1l1P21
24、Sink(si,pi)(+,cij)(1,1),(1,li)(1,0)SourceS1S2S7A1A2A15P1P2Sink(si,pi)(+,cij)(li,f(f+1)/2)(li,0)线性费用网络线性费用网络(只有产量上限只有产量上限)非线性费用网络非线性费用网络(只有产量上限只有产量上限)边的标记（流量上限，单位费用）用标准算法(如最小费用路算法)求解无单位费用概念(f(f+1)/2),需修改最小费用路算法2）最小费用网络流模型）最小费用网络流模型产量有下限ri时的修正SourceSiSi(si-ri,pi)(ri,0)(+,0)得到的结果应加上 iiipr 才是最小费用的的处处理理同
25、同前前约约束束条条件件)7,2,1(,5000151 isxijij注：该模型获当年的惟一最高奖（网易杯）注：该模型获当年的惟一最高奖（网易杯）S1S2S3S6S5S1S2S2S3S3S5S5S63)最小面积模型A1A2A3A4A5A6A7A8A9A10A11A12A13A14A15cx作图：Si到管道x单位钢管的最小购运费用c由各条Si首尾相连(横坐标)组成的一条折线对应一个购运方案，折线下面的面积对应方案的费用在产量约束下找面积最小的折线问题问题2:分析对购运计划和总费用影响分析对购运计划和总费用影响(哪个钢厂销价哪个钢厂销价变化影响最大；哪个钢厂产量上限变化影响最大变化影响最大；哪个钢厂
26、产量上限变化影响最大)规划问题的灵敏度分析问题问题3：管道为树形图：管道为树形图701088107062300220210170690462160320160A10A11A12S4S5S6130A17A18A19A20190260100,0,500,0.)(05.0min71)(2112211)(71211 jkijjkkjjkijkEjkijjiijjkjkjEjkijijijyxbyyyxsxtsyyxc(jk)是连接Aj,Ak的边，E是树形图的边集，ljk是(jk)的长度，yjk是由Aj沿(jk)铺设的钢管数量论文中发现的主要问题论文中发现的主要问题1）针对题目给的数据用凑的方法算出结果
27、，没有解决这类问题的一般模型2）局部最优，如将管道分为左右两段，分别寻求方案；如将问题分为购运和铺设两部分，分别寻优（会导致每段管道都从两端铺到中点）的的处处理理不不正正确确约约束束条条件件,5000)3151ijijsx 4）由Si至Aj的最小购运费用路线及最小费用cij 不对5）数字结果相差较大(如最小费用应127.5至128.2亿元)1995年全国大学生数学建模竞赛年全国大学生数学建模竞赛A题题(CUMCM-1995A)一个飞行管理问题一个飞行管理问题 在约在约10000m高空的某边长高空的某边长160km的正方形区域内的正方形区域内,经常有若干架飞机作水平飞行经常有若干架飞机作水平飞行
28、,区域内每架飞机区域内每架飞机的位置和速度向量均由计算机记录其数据的位置和速度向量均由计算机记录其数据,以便以便进行飞行管理进行飞行管理.当一架欲进入该区域的飞机到达当一架欲进入该区域的飞机到达边界区域边缘时边界区域边缘时,记录其数据后记录其数据后,要立即要立即 计算并判计算并判断是否会与其区域内的飞机发生碰撞断是否会与其区域内的飞机发生碰撞.如果会碰如果会碰撞撞,则应计算如何调整各架则应计算如何调整各架(包括新进入的包括新进入的)飞机飞飞机飞行的方向角行的方向角,以避免碰撞以避免碰撞.现假设条件如下现假设条件如下:一个飞行管理问题一个飞行管理问题 1)不碰撞的标准为任意两架飞机的距离大于不碰
29、撞的标准为任意两架飞机的距离大于8km;2)飞机飞行方向角调整的幅度不应超过飞机飞行方向角调整的幅度不应超过30度度;3)所有飞机飞行速度均为每小时为所有飞机飞行速度均为每小时为800km;4)进入该区域的飞机在到达区域边缘时进入该区域的飞机在到达区域边缘时,与区域内与区域内飞机的距离应在飞机的距离应在 60km以上以上;5)最多考虑最多考虑6架飞机架飞机;6)不必考虑飞机离开此区域后的状况不必考虑飞机离开此区域后的状况;请你对这个避免碰撞的飞行管理问题建立请你对这个避免碰撞的飞行管理问题建立数学模型数学模型.列出计算步骤列出计算步骤,对以下数据进行计算对以下数据进行计算(方方向角误差不超过向
30、角误差不超过0.01度度),要求飞机飞行方向角调要求飞机飞行方向角调整的幅度尽量小整的幅度尽量小.飞机编号 横坐标x 纵坐标y 方向角(度)1 150 140 243 2 85 85 236 3 150 155 220.5 4 145 50 159 5 130 150 230 新进入 0 0 52注:方向角指飞行方向与x轴正向的夹角试根据实际应用背景对你的模型进行评价与推广试根据实际应用背景对你的模型进行评价与推广.区域区域4个顶点坐标为个顶点坐标为(0,0),(160,0),(160,160),(0,160).记录数据为记录数据为:两架飞机不碰撞条件两架飞机不碰撞条件 初始位置 时刻t飞机的
31、位置两架飞机的距离222)()()()(tytytxtxrjijiij),(00iiyx),(00jjyx),k(,sin)(,cos)(k0kkk0kkjivtytyvtxtx,)()()sin)(sin()cos)(cos(2)(sin)cos(cos200200000022222jijijijijijijijiijyyxxtyyxxvtvr),sin)(sin()cos)(cos(2,)(sin)cos(cos0000222jijijijiijjijiijyyxxvbva).0()(222ijijijijrtbtatr.64)0()(222ijijijijrtbtatr引入记号 距离表示
32、为 两架飞机不碰撞条件),min(jiijttt 00000000000000000000tan,223tan,23,sin,tan,23tan,2,cos,tan,2tan,20,sin,tan,223tan,20,cosiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiixDyorxyifvyxyorxyDifvxxyDorxDyDifvxDxDyorxDyDifvxDt不必考虑在区域外的碰撞两架飞机都在区域中的时间其中第i架飞机在区域内的时间0|)(02ijtijbtrdtd0ijb.64)0()(222ijijijijrtbtatr64)0(2ijijrc
33、.0)(22ijijijijctbtatr042ijijijcab不碰撞条件的另一表述：初始不碰,距离增加距离减少可能碰撞 发生碰撞必须引入记号 无实根,不可能碰撞,条件为,042ijijijcabijijijijijijacabbT22)(),(),(),(ijjijjijiijiTyTxTyTx,)(0)(,)(0)(,)(0)(,)(0)(DTyorTyDTxorTxDTyorTyDTxorTxijjijjijjijjijiijiijiiji即使有实根,碰撞也不必考虑(不在区域内碰撞),有实根条件 符合条件的实根为 可能碰撞的位置 碰撞不在区域内的条件,04,0)3(04,0)2(;0)
34、1(22ijijijijijijijijijcabbcorcabbcorbc,)(0)(,)(0)(,)(0)(,)(0)(DTyorTyDTxorTxDTyorTyDTxorTxijjijjijjijjijiijiijiiji归结起来不在区域内碰撞的条件为 数学模型数学模型初始位置与方向角000),(iiiyxiii0.),(121NiiNF.,1,6,1,64)(.,),(min2121NijiNjitttrtsFiijijNiiN调整后的方向角 总的调整量 优化数学模型.0)(,0,0,6.,),(min)3(.0,0,6.,),(min)2(.0,6.,),(min)1(2121121
35、21121ijiijijijijiNiiNijijijijiNiiNijiNiiNTxcabbtsFcabbtsFbtsF对应两架飞机不碰撞条件的另一表述的优化模型对应两架飞机不碰撞条件的另一表述的优化模型vDtm/2.,1,6,1,64)(.,),(min2121NijiNjitttrtsFimijNiiN一个简化化数学模型一个简化化数学模型任何一架飞机在区域中停留最长时间 放松到任两架飞机在这段时间不碰撞优化模型简化为.,1,6,1,64)(.,),(min2121NijiNjitrtsFiijNiiN进一步简化化的数学模型进一步简化化的数学模型(永远不碰撞永远不碰撞)可化为线性规划问题0
36、,0,0,6,6,6,3,2,1,8)(.min0,0,0,6,6,6,3,2,1,8)(.min0,0,0,6,6,6,3,2,1,8)(.min321321321321321321321321321jijitrtsjijitrtsjijitrtsijijij目标函数的线性化目标函数的线性化,化为化为8个优化问题之三个优化问题之三个个(N=3)22sin(),22(cos(2sin2)2cos,2sin(2sin2)2cos2sin,2sin2sin(2)sinsin,coscos(jijijijijijijijijijiijijvvvvvvvvv0ij0ijr相对运动的观点相对运动的观点2200jjiiij00102,)0(8sinijijijijijr碰撞临界相对运动方向角以及不碰撞条件碰撞临界相对运动方向角以及不碰撞条件 NitsFiNiiN,1,6.),(min12122j00jiiij最优解最优解 Thats all.Any Questions?Thank you for your attendance!谢金星,清华大学数学科学系,2007.
