二次模拟论文2.doc

(一)、摘要高台村用水短缺问题长期以来十分严重，为了解决高台村的用水问题，县政府决定从今年起为该村进行打井和管输建设。本文首先根据高台村现有水井近9年的产水量，通过MATLAB数据拟合分析研究其产水规律，并预测了其在未来五年内的产水情况，再根据该村未来五年的需水量、县政府提供的资金情况等，通过建立多目标规划模型，对该村解决用水问题的方案进行了规划，并给出了具体的解决方案。在建立多目标规划模型时，对定性因素选用0-1变量，以花费资金和每年产水量作为目标，分别设置偏差变量，并给定合适的优先级和权重进行处理，满足目标和约束，得到合理的打井和管输建设计划: 打井计划：第一年打1,3,6,7号井；第二年打4,5号井；第三年打2号井。管输建设：第一年建8.580165公里；第二年建1.430028公里；第三年建10.01019公里。最后将模型进行推广，使其尽可能长的时间满足村民的需水量。关键字 数据拟合 0-1变量 多目标规划 LINGO(二)、问题重述我国许多农村用水短缺问题长期以来十分重，别是一些偏远贫困地区况更加严重。如我国西北地区偏远贫困村高台村，年均降雨不足20mm，是典型的缺水地区，过去该村的用水主要源于两种途经：一是靠每家每户自行建造的蓄水池屯积雨水；二是利用村里现有的四口水井。但由于近年来能源开采,使水源遭到破坏，经常出现数月滴雨不下，使蓄水池的功能完全丧失，同时使原有的四4口水井产水量逐渐减少。表1给出四口水井近9年的产水量统计数字。 为了解决高台村的用水问题，县政府决定从今年起为该村解决用水问题，具体方案包括两个方面：一是在高台村附近找8个可供打井的位置，但由于位置不同其地质构造也不同，因此每个位置打井的费用和年产水量也不同，详见表2，而且预计每口水井的年产水量还会以平均每年10%左右的速率减少。二是从长远考虑，通过建造输水管道，从距该20公里外的柳河将水引入该村，管输的费用为（万元），其中表示每年的可输水量（万吨/年），表示管道长度（公里）。管输工程工期需要三年时间，且每年投资管输的费用为万元的整数倍。要求管输完成之后，每年能够管输至少100万吨水。县政府从2012年开始，连续三年，每年最多可提供60万元用于该村打井和管输工程，为了保证该村从2012至2016年这五年间每年分别能至少获得150、160、170、180、190万吨水。请用数学建模的方法，在不考虑小蓄水池的作用和利息的因素的情况下，为县政府作出一个从2012年起三年的打井和管输建设计划。（三）、问题分析针对高台村的缺水问题，我们可以采用多目标规划的方法解决，即要求花最小的投资获得最大的产水量。首先我们需要设置优先级，假定投资的优先级为p1，产水量的优先级为p2；每年的投资即为每年所建的井乘以相应的费用，再加上每年管输的建设费用；2012年以后高台村的供水（不考虑小蓄水池的作用）主要来自以前的四口井，新建成的井，以及管输建成后的输水。因此高台村的水源可分为三部分考虑。针对第一部分，我们必须预测出以前四口水井从2012年到2016年的产水量。我们可以采用均值算出各水井水量的平均减少率，然后求出相应的水量。针对第二部分，我们用0-1变量先确定前三年所打的井，然后乘以相应的产水量（产水量平均以10%的速率减少）针对第三部分，由于一旦管输修建成功，年产水量将不再改变，故只需要加上管输的年产水量即可。三部分相加即为每年的供水量。其次我们需要设置偏差变量，要使投资最小即要求每年相应正的偏差变量最小，而要使产水量最大即要求每年相应产水量的负的偏差变量最小。而对于同一优先级我们分配以相同的权重因子。最后我们需要统一处理目标与约束条件，整合所建立的模型，利用LINGO软件进行求解。 （四）、模型假设（1）、不考虑小蓄水池的作用和利息的因素；（2）、每年投资的资金当年没有用完则可以累积到下一年；（3）、井和管输建成后当年即可产水；（4）、不考虑其他突发因素导致的费用增加；（5）、不考虑自然因素（降雨）对供水量的影响（五）、符号说明（1）、； （2）、Yi（i=1,2,3,）表示第i年修建的管输长度；（3）、T为管输修建成功后每年的可输水量；（4）、（i=1，2，3）为整型变量，表示第i年投资的管输费用；（5）、（j=1,28）表示修建编号为j的井的费用；（6）、（j=1,28）表示第j号井修建成功当年的产水量；（7）、表示第i年原有第j号井的产水量。（六）、模型建立与求解由表一的数据我们可以借用MATLAB软件拟合出各口井的产水量关于年份的函数方程。为精确起见，对不同的井我们选用不同的函数进行拟合，如一次函数、三次函数。（具体拟合程序见附件） 根据拟合出的函数，进而预测出2012年到2016年各口井的产水量，具体数据如下：年份产水量编号201220132014201520161号井21.420.319.117.816.62号井1.100003号井9.17.14.92.90.84号井11.84.5000总产水量43.431.824.020.717.4由于本题采用多目标规划，故我们需要设置偏差变量，决策变量的含义入下：分别为第一年，前两年，前三年超过投资的金额；分别为第一年，第二年，第三年不足投资的金额；分别为2012年到2016年各年超过需水量的部分；分别为2012年到2016年各年不足需水量的部分。（i=1,2,3；j=1,28）由于每口井最多只能打一次，所以有约束： ；Yi表示第i年修建管输的长度，所以对管道长度有约束： ；由于以及表二的数据，我们可以得到投资、产水量与决策变量的关系：由于假设每年投资的资金当年没有用完可以累积到下一年则第一年的投资应小于60万元，第一年与第二年的投资之和小于120万元，三年投资的综合小于180万元。而对管输的费用必须是万元的整数倍，所以对投资的约束如下：我们已经预测出2012年到2016年原有的四口井的产水量，而新打的井的年产水量平均每年以10%的速率减少，由Xij我们可以算出水井的产水量。由于管输的工期为三年，故只有后两年才有管输的产水量根据2012年到2016年该村的需水量，对产水量的约束如下：；整理该模型如下：目标函数：MIN ；约束条件：；gin（Z1）；gin（Z2）；gin（Z3）；然后根据优先级将目标规划问题分解为一系列的单目标规划问题，再依次求解。我们可将上面的程序分为两部分输入LINGO软件中，即可得到该问题的解决办法。（具体结果见附录）结果分析：由LINGO程序运行结果可设计一个从2012年起三年的打井和管输建设的最优计划：打井计划：第一年打1,3,6,7号井；第二年打4,5号井；第三年打2号井。管输建设：第一年建8.580165公里；第二年建1.430028公里；第三年建10.01019公里。建设费用和以后五年的供水量如下表：20122013201420152016建井费用20107管输费用137750建设总费用338757供水量150.1174.1188.7268.4250.3模型改进由于解决偏远山区用水困难的问题并不是一项以盈利为目的的工程，我们首先要考虑的还是村民的用水问题，于是考虑将所有投资用于打井与管输建设计划，以使供水量更长时间的满足村民的用水需求（即将模型的适用范围延长）。对此，我们算出上述模型在未来10年的供水量，数据如下表：2012201320142015201620172018201920202021原有43.431.92420.717.415.414.213.111.810.6现有107142.3164.1147.7132.9119.6107.696.887.178.4管输000100100100100100100100总量150.4174.2188.1268.4250.3235221.8209.9198.9189根据村民的需水量，我们可以知道其是递增的，但是需水量并不会无限制的增大，于是我们可假设2017年的需水量为200（万吨），2018年的需水量为210（万吨），而以后几年的需水量一直保持210（万吨）。则由上表可知此模型的解决方案2019年的供水量并不能满足村民的需水量。于是我们对此模型进行了改进，使其满足2019年的供水量，得到的结果如下：（七）、模型评价（1）本文把所解决的问题归结为优化问题，建立的数学模型清晰合理。（2）运用MATLAB和LINGO软件处理数据和进行运算，降低运算量，简单易行，有很大的可操作性，且所得数据较为合理可靠。(3）运用01模型解题，全面可靠（4）但在实际运用本方案中还应考虑自然因素对产水量的影响，还有需水量的变化，根据实际情况进行灵活改变。参考文献1姜启源,谢金星,叶俊.数学模型M.北京:高等教育出版社,2007.2江世宏.MATLAB语言与数学实验M.北京:科学出版社，2007.附录MATLAB程序井1：y=32.2,31.3,29.7,28.6,27.5,26.1,25.3,23.7,22.7; x=1:1:9; plot(x,y,r+) p1=polyfit(x,y,1);t=1:1:14 z=polyval(p1,t);plot(x,y,r+,t,z,r) title(一号井产水量拟合图形) xlabel(年份)ylabel(产水量) x=10:1:14; y1=polyval(p1,x)y1 = 21.4472 20.2456 19.0439 17.8422 16.6406井2： y=21.5,15.9,11.8,8.7,6.5,4.8,3.5,2.6,2.0; x=1:1:9; plot(x,y,r*) p2=polyfit(x,y,3); z=polyval(p2,x); t=1:1:14; z=polyval(p2,t); plot(x,y,r*,t,z,r) xlabel(年份) title(二号井产水量拟合图形) ylabel(产水量) x=10:1:14; y2=polyval(p2,x)y2 =1.0690 -0.0595 -1.6448 -3.8732 -6.9310x=1:1:9;井3： y=27.9,25.8,23.8,21.6,19.5,17.4,15.5,13.3,11.2; x=1:1:9; plot(x,y,r+) p3=polyfit(x,y,1); z=polyval(p3,x); t=1:1:14; z=polyval(p3,t); plot(x,y,r*,t,z,r) title(三号井产水量拟合图形)xlabel(年份)ylabel(产水量) x=10:1:14; y3=polyval(p3,x)y3 = 9.1306 7.0456 4.9606 2.8756 0.7906井4：井4; y=46.2 32.6 26.7 23.0 20.0 18.9 17.5 16.3 ; x=2:1:9; plot(x,y,r+) p4=polyfit(x,y,3); z=polyval(p4,x); t=2:1:14; z=polyval(p4,t); plot(x,y,r+,t,z,r) title(四号井产水量拟合图形)xlabel(年份) ylabel(产水量) x=10:1:14; y4=polyval(p1,x)y4= 11.8286 4.5071 -7.0476 -23.8994 -47.1117LINGO程序Min=DPLUS1+ DPLUS2+ DPLUS3;X11+X21+X31=1;X12+X22+X32=1;X13+X23+X33=1;X14+X24+X34=1;X15+X25+X35=1;X16+X26+X36=1;X17+X27+X37=1;X18+X28+X38=20;T=100;Z1=0.667*T0.51*Y1;Z2=0.667*T0.51*Y2;Z3=0.667*T0.51*Y3;gin(Z1);gin(Z2);gin(Z3);5*X11+7*X12+5*X13+4*X14+6*X15+5*X16+5*X17+3*X18+0.667*T0.51*Y1-DPLUS1+DMINUS1=60;5*X11+7*X12+5*X13+4*X14+6*X15+5*X16+5*X17+3*X18+0.667*T0.51*Y1+5*X21+7*X22+5*X23+4*X24+6*X25+5*X26+5*X27+3*X28+0.667*T0.51*Y2-DPLUS2+DMINUS2=120;5*X11+7*X12+5*X13+4*X14+6*X15+5*X16+5*X17+3*X18+0.667*T0.51*Y1+5*X21+7*X22+5*X23+4*X24+6*X25+5*X26+5*X27+3*X28+0.667*T0.51*Y2+5*X31+7*X32+5*X33+4*X34+6*X35+5*X36+5*X37+3*X38+0.667*T0.51*Y3-DPLUS3+DMINUS3=180;43.1+25*X11+36*X12+32*X13+15*X14+31*X15+28*X16+22*X17+12*X18-DPLUS4+DMINUS4=150;31.8+25*X21+36*X22+32*X23+15*X24+31*X25+28*X26+22*X27+12*X28+25*X11*0.9+36*X12*0.9+32*X13*0.9+15*X14*0.9+31*X15*0.9+28*X16*0.9+22*X17*0.9+12*X18*0.9-DPLUS5+DMINUS5=160;24.0+25*X31+36*X32+32*X33+15*X34+31*X35+28*X36+22*X37+12*X38+25*X21*0.9+36*X22*0.9+32*X23*0.9+15*X24*0.9+31*X25*0.9+28*X26*0.9+22*X27*0.9+12*X28*0.9+25*X11*0.81+36*X12*0.81+32*X13*0.81+15*X14*0.81+31*X15*0.81+28*X16*0.81+22*X17*0.81+12*X18*0.81-DPLUS6+DMINUS6=170;20.7+25*X31*0.9+36*X32*0.9+32*X33*0.9+15*X34*0.9+31*X35*0.9+28*X36*0.9+22*X37*0.9+12*X38*0.9+25*X21*0.81+36*X22*0.81+32*X23*0.81+15*X24*0.81+31*X25*0.81+28*X26*0.81+22*X27*0.81+12*X28*0.81+25*X11*0.729+36*X12*0.729+32*X13*0.729+15*X14*0.729+31*X15*0.729+28*X16*0.729+22*X17*0.729+12*X18*0.729+T-DPLUS7+DMINUS7=180;17.4+25*X31*0.81+36*X32*0.81+32*X33*0.81+15*X34*0.81+31*X35*0.81+28*X36*0.81+22*X37*0.81+12*X38*0.81+25*X21*0.729+36*X22*0.729+32*X23*0.729+15*X24*0.729+31*X25*0.729+28*X26*0.729+22*X27*0.729+12*X28*0.729+25*X11*0.6561+36*X12*0.6561+32*X13*0.6561+15*X14*0.6561+31*X15*0.6561+28*X16*0.6561+22*X17*0.6561+12*X18*0.6561+T-DPLUS8+DMINUS8=190;结果 Local optimal solution found. Objective value: 0.000000 Objective bound: 0.000000 Infeasibilities: 0.000000 Extended solver steps: 0 Total solver iterations: 26 Model Class: MINLP Total variables: 47 Nonlinear variables: 4 Integer variables: 27 Total constraints: 22 Nonlinear constraints: 6 Total nonzeros: 211 Nonlinear nonzeros: 15 Variable Value DPLUS1 0.000000 DPLUS2 0.000000 DPLUS3 0.000000 X11 0.000000 X21 1.000000 X31 0.000000 X12 0.000000 X22 1.000000 X32 0.000000 X13 0.000000 X23 1.000000 X33 0.000000 X14 0.000000 X24 1.000000 X34 0.000000 X15 0.000000 X25 1.000000 X35 0.000000 X16 0.000000 X26 1.000000 X36 0.000000 X17 0.000000 X27 1.000000 X37 0.000000 X18 0.000000 X28 1.000000 X38 0.000000 Y1 8.580165 Y2 1.430028 Y3 10.01019 T 100.2395 Z1 60.00000 Z2 10.00000 Z3 70.00000 DMINUS1 0.000000 DMINUS2 10.00000 DMINUS3 0.000000 DPLUS4 0.000000 DMINUS4 106.9000 DPLUS5 72.80000 DMINUS5 0.000000 DPLUS6 34.90000 DMINUS6 0.000000 DPLUS7 103.7495 DMINUS7 0.000000 DPLUS8 74.16847 DMINUS8 0.000000 Row Slack or Surplus Dual Price 1 0.000000 -1.000000 2 0.000000 0.000000 3 0.000000 0.000000 4 0.000000 0.000000 5 0.000000 0.000000 6 0.000000 0.000000 7 0.000000 0.000000 8 0.000000 0.000000 9 0.000000 0.000000 10 0.2038527E-01 0.000000 11 0.2394711 0.000000 12 0.000000 0.000000 13 0.000000 0.000000 14 0.000000 0.000000 15 0.000000 0.000000 16 0.000000 0.000000 17 0.000000 0.000000 18 0.000000 0.000000 19 0.000000 0.000000 20 0.000000 0.000000 21 0.000000 0.000000 22 0.000000 0.000000第二次程序Min=DMINUS4+DMINUS5+DMINUS6+DMINUS7+DMINUS8;X11+X21+X31=1;X12+X22+X32=1;X13+X23+X33=1;X14+X24+X34=1;X15+X25+X35=1;X16+X26+X36=1;X17+X27+X37=1;X18+X28+X38=20;T=100;Z1=0.667*T0.51*Y1;Z2=0.667*T0.51*Y2;Z3=0.667*T0.51*Y3;gin(Z1);gin(Z2);gin(Z3);5*X11+7*X12+5*X13+4*X14+6*X15+5*X16+5*X17+3*X18+0.667*T0.51*Y1-DPLUS1+DMINUS1=60;5*X11+7*X12+5*X13+4*X14+6*X15+5*X16+5*X17+3*X18+0.667*T0.51*Y1+5*X21+7*X22+5*X23+4*X24+6*X25+5*X26+5*X27+3*X28+0.667*T0.51*Y2-DPLUS2+DMINUS2=120;5*X11+7*X12+5*X13+4*X14+6*X15+5*X16+5*X17+3*X18+0.667*T0.51*Y1+5*X21+7*X22+5*X23+4*X24+6*X25+5*X26+5*X27+3*X28+0.667*T0.51*Y2+5*X31+7*X32+5*X33+4*X34+6*X35+5*X36+5*X37+3*X38*0.667*T0.51*Y3-DPLUS3+DMINUS3=180;43.1+25*X11+36*X12+32*X13+15*X14+31*X15+28*X16+22*X17+12*X18-DPLUS4+DMINUS4=150;31.8+25*X21+36*X22+32*X23+15*X24+31*X25+28*X26+22*X27+12*X28+25*X11*0.9+36*X12*0.9+32*X13*0.9+15*X14*0.9+31*X15*0.9+28*X16*0.9+22*X17*0.9+12*X18*0.9-DPLUS5+DMINUS5=160;24.0+25*X31+36*X32+32*X33+15*X34+31*X35+28*X36+22*X37+12*X38+25*X21*0.9+36*X22*0.9+32*X23*0.9+15*X24*0.9+31*X25*0.9+28*X26*0.9+22*X27*0.9+12*X28*0.9+25*X11*0.81+36*X12*0.81+32*X13*0.81+15*X14*0.81+31*X15*0.81+28*X16*0.81+22*X17*0.81+12*X18*0.81-DPLUS6+DMINUS6=170;20.7+25*X31*0.9+36*X32*0.9+32*X33*0.9+15*X34*0.9+31*X35*0.9+28*X36*0.9+22*X37*0.9+12*X38*0.9+25*X21*0.81+36*X22*0.81+32*X23*0.81+15*X24*0.81+31*X25*0.81+28*X26*0.81+22*X27*0.81+12*X28*0.81+25*X11*0.729+36*X12*0.729+32*X13*0.729+15*X14*0.729+31*X15*0.729+28*X16*0.729+22*X17*0.729+12*X18*0.729+T-DPLUS7+DMINUS7=180;17.4+25*X31*0.81+36*X32*0.81+32*X33*0.81+15*X34*0.81+31*X35*0.81+28*X36*0.81+22*X37*0.81+12*X38*0.81+25*X21*0.729+36*X22*0.729+32*X23*0.729+15*X24*0.729+31*X25*0.729+28*X26*0.729+22*X27*0