工厂选址问题最终稿.doc

承 诺 书我仔细阅读了中国大学生数学建模竞赛的竞赛规则。我完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我将受到严肃处理。参赛选择的题号是（从A/B/C中选择一项填写）： A 报名号： 所属学院、专业、学号、报名号： 承诺人姓名 ： 郝培良 冀航 李贤伟 日期： 2014 年 5月 8 日A题 公司新厂选址问题摘要在资源稀缺的市场竞争时代,如何优化资源配置是每个生产公司在日益激烈的市场竞争中求生存、促发展的有效途径和理智选择。本案例是根据其上一年11月向各城市供应的需求量，来预测下一年对各城市的供应量。并且新厂还要根据对各城市的供应量来扩大生产规模。其生产规模则可选择一个生产量最大的城市来算出其生产规模。模型I 根据所给供应量对未来一年的供应量的做出预测。先用SPSS对所给数据的散点图进行曲线回归，通过比较不同回归曲线的sig值，得到最适合各个城市需求量的拟合曲线为：（三次多项式）。再用MATLAB进行数据处理，建立三次多项式拟合的模型，画出拟合曲线图，并求出月份相对于城市需求量的回归方程，来预测其走势。其拟合效果，大体上符合短期预测所给出的未来一年各个城市的需求量。模型II 根据18个城市每个月的需求量总和，我们选出11个月中最大生产量的月来对生产规模进行规划。根据所给的数据求出其成本最小化。我们建立了动态规划模型，并且将动态规划模型转化为0-1规划模型。利用LINGO进行LP优化求解得到成本最小值，并用以确定各个厂的生产规模。为了方便建模，我们忽略了交叉工资弹性和运输竞争。生产成本=月总工资+货物量的运输成本。求出各厂的生产量，与各厂大概人数。模型III 根据第二问的0-1变量，确定出各个工厂所供应的城市。在需求点（各个城市的坐标）、运输量及线性运输费率不变的前提下，我们将运输成本作为唯一决策因素。于是对该静态选址问题建立模型。通过MATLAB，利用重心法为各个工厂确定最优选址，对运输路径进行优化。得到其运输成本最优解。关键词 三次多项式拟合 0-1规划模型 LP优化 重心法一、 问题的重述公司现场选址问题公司要根据各城市的需求量、各工厂到各城市的距离和工资标准，对工厂进行规划和重新选址，使公司的成本最低，获利最大。现给出18个城市某年一到十一月的月需求量，需建立模型预测未来一年中各城市每月的产品需求量；通过给定的六个加工厂距各城市的距离和六个工厂所在的工资标准及运输价格等信息，在满足约束条件的情况下，确定各加工厂的生产规模，将各城市合理分配给六个工厂，确定各工厂的人数，使运输费用和工资最小；进一步优化，根据相关条件，重新设定新厂址位置，让运输费用更低，并向公司提出有用的建议。二、 问题的分析问题一：我们将18个城市11个月的月需求量数据,通过SPSS绘制散点图并进行曲线回归,通过对不同曲线的sig值的比较，最终决定用三次多项式进行拟合。在MATLAB中输入所给的数据，采用多项式拟合的方法进行预测，对其进行多项式拟合，预测下年的各个地区每月的产品需求量，建立模型求解。问题二：根据所给工资标准及运输价格等条件确定各加工厂的生产规模，生产规模应该是指工厂的人数，根据所给的运输路程和基本工资，结合其约束条件的规划模型，可以运用规划优化模型进行求解。由于案例提供了各城市的月需求量、各工厂到各城市的距离和工资标准，我们可以建立以生产量和加班时间为决策变量；通过决策变量来表示工厂的人数，构建工资和运输成本最低的目标函数；再根据附件数据可以找出一个约束条件：工厂对每个城市的供应量不能超过它的生产量，再通过国家对每月加班时间的要求和我们的合理假设（每个城市只由一个工厂提供产品）再得到两个约束，在目标函数和约束条件都合乎逻辑的情况下，建立规划模型。问题三：在问题二中我们得到了一个优化后的结果，即：每个工厂所供应的城市。在需求点（各个城市的坐标）、运输量及线性运输费率不变的前提下，我们将运输成本作为唯一决策因素。于是对该静态选址问题建立模型。通过MATLAB，利用重心法为各个工厂确定最优选址，对运输路径进行优化。得到最短运输距离，从而使使运输成本降低。三、 基本假设模型一：设做短期预测，波动大合理模型二：1.设工人每天正常工作时间为8小时，每月只工作24天； 2.设加班时间每月不超过36小时； 3.设每个城市由一个厂提供； 4.设新厂的生产量不受最大容量影响； 5.工厂对每个城市的供应量不能超过它的生产量； 6.不考虑工资弹性和交叉工资弹性。模型三：1.设运输成本只受运输距离影响，不受需求量影响 2.货物运输正点交付，不存在货损和货差。3.不考虑运输竞争。 四、 符号说明：总工资；：总得运输成本；：第i个工厂当地的标准工资；: 第i个工厂的生产量；：第i个工厂的加班时间；：第i个工厂的员工人数；：第j个城市五月份的需求量；：第i个工厂到第j个城市的距离；：第j个城市是否由第i个工厂提供产品。五、 模型的建立与求解问题一：模型I 通过对实际数据的分析，发现1月到11月各城市需求量的起伏波动无规律，并且又是对短时间内的需求量做预测，所以找不到一种可以使其更精确的预测方法。于是我们用SPSS描绘了数据的散点图，并且采用回归曲线的方式建立了多个拟合模型。 注：源程序及相关说明见附件。模型汇总和参数估计值因变量:天津需求量方程模型汇总参数估计值R 方Fdf1df2Sig.常数b1b2b3线性.048.45019.519360876.3645660.000对数.2402.83519.127304934.63656502.206倒数.4948.78119.016456300.983-223886.847二次.6597.74128.013172040.00092815.245-7262.937三次.88718.28237.001-8222.727241713.578-36976.5731650.758复合.1011.01219.341331885.2521.025幂.3404.64119.060279571.968.200S.62014.65819.00413.064-.746增长.1011.01219.34112.713.025指数.1011.01219.341331885.252.025Logistic.1011.01219.3413.013E-6.976自变量为 月份。通过比较不同回归曲线的sig值，我们发现三次多项式曲线和S型曲线的拟合程度明显高于其它曲线，考虑到S型曲线拟合建模的难度较大，最终我们建立了三次多项式拟合模型，即。根据模型，由MATLAB求得18个城市未来一年需求量的图（散点：实际值 红色曲线：预测值）一 天津后一年各个月份的预测值为：419985，511504，657841，868902，1154593，1524820，1989489 2558506，3241777，4049208，4990705，6076174，7315521二 太原后一年各个月份的预测值为：496318，687132，969214，1357300，1866126，2510428，3304942 4264404，5403550，6737116，8279838，10046452三 石家庄后一年各个月份的预测值为：410790，442640，503430，599460，737030，922440，1161990，1461980 1828710，2268480，2787590，3392340，4089030四 济南后一年各个月份的预测值为（1.0e+006）*：0.3807，0.3732，0.3798，0.4048，0.4524，0.5268，0.6323，0.7730，0.9533，1.1773，1.4492，1.7734，2.1540 五 郑州后一年各个月份的预测值为（1.0e+006 *）：0.4663，0.5369，0.6492，0.8110，1.0304，1.3151，1.6730，2.1121，2.6402，3.2652，3.9950，4.8375，5.8006六 西安后一年各个月份的预测值为（1.0e+006 *）：0.5798，0.7277，0.9395，1.2258，1.5972，2.0642，2.6374，3.3273，4.1446，5.0997，6.2033，7.4659，8.8981七 上海后一年各个月份的预测值为（1.0e+006 *）：0.4782，0.5935，0.7640，0.9993，1.3085，1.7011，2.1865，2.7740，3.4730，4.2929，5.2430，6.3327，7.5713八 南京后一年各个月份的预测值为（1.0e+007*）：0.0716，0.0947，0.1265，0.1684，0.2216，0.2874，0.3672，0.4623，0.5738，0.7032，0.8518，1.0208，1.2115九 合肥后一年各个月份的预测值为（1.0e+007*）：0.0691，0.0913，0.1224，0.1637，0.2165，0.2822 0.3623，0.4579，0.5705，0.7015，0.8522，1.0239，1.2180十 武汉后一年各个月份的预测值为（1.0e+006*）：0.3817，0.4013，0.4416，0.5071，0.6025，0.7325，0.9017，1.1146，1.3760，1.6904，2.0626，2.4970，2.9984十一 重庆后一年各个月份的预测值为（1.0e+006*）：0.5211，0.6576，0.8569，1.1294，1.4852，1.9345，2.4876，3.1547，3.9460，4.8717，5.9421，7.1673，8.5576十二 杭州后一年各个月份的预测值为（1.0e+007*）：0.0630，0.0859，0.1183，0.1617，0.2175，0.2872，0.3723，0.4742，0.5944，0.7344，0.8955，1.0794，1.2874十三 长沙后一年各个月份的预测值为（1.0e+007*）：0.0732，0.1006，0.1378，0.1864，0.2478，0.3233 ，0.4144，0.5226，0.6491，0.7955，0.9632，1.1535，1.3679十四 南昌后一年各个月份的预测值为（1.0e+006*）：0.4905，0.6190，0.8012，1.0454，1.3598，1.7528，2.2327，2.8078，3.4865，4.2770，5.1877，6.2269，7.4028十五 贵阳后一年各个月份的预测值为（1.0e+006*）：0.4905，0.6190，0.8012，1.0454，1.3598，1.7528，2.2327，2.8078，3.4865，4.2770，5.1877，6.2269，7.4028十六 福州后一年各个月份的预测值为（1.0e+007*）：0.0627，0.0828，0.1103，0.1465，0.1924，0.2491，0.3177，0.3993，0.4949，0.6057，0.7327，0.8771，1.0398十七 广州后一年各个月份的预测值为（1.0e+006*）：0.4540，0.5578，0.7097，0.9175，1.1887，1.5310，1.9520，2.4594，3.0607，3.7636，4.5757，5.5047，6.5582十八 南宁后一年各个月份的预测值为（1.0e+007*）：0.0663，0.0944，0.1335，0.1852，0.2512，0.3331，0.4323，0.5506，0.6896，0.8508，1.0358，1.2463，1.4838 问题二：模型II本题要解决的是各加工厂的员工人数、加班时间和生产量以及工厂合理配置的问题。由于案例提供了各城市的月需求量、各工厂到各城市的距离和工资标准，在忽略交叉工资弹性和运输竞争的情况下，我们可以建立以生产量和加班时间为决策变量。通过决策变量来表示厂的人数，构建工资和运输成本最低的目标函数；再通过国家对每月加班时间的要求和我们的合理假设（每个城市只由一个工厂提供产品）再得到两个约束，在目标函数和约束条件都合乎逻辑的情况下，建立线性规划模型，并用LINGO进行优化。源程序见附录。通过计算，可以轻松地得到5月份各个城市的需求量总和最大。设：每个厂的五月份的产量为：；每个厂五月份的加班时间为： 各个工厂的人数： ；各个工厂的总工资：；各个工厂的总运输成本：； 我们建立模型如下：目标函数： 约束条件：工厂对每个城市的供应量不能超过它的生产量i工厂提供j城市产品为1，否则为0运用lingo得到以下结果：各个厂五月份的生产量：x1=190640000;x2=187830000;x3=89980000;x4=93260000;x5=128370000;x6=142960000;需要加班的工厂：t2=36;t4=36;各个工厂的人数规模：a1 =1.2411e+005；a2 =1.02977e+005；a3 =5.8581e+004；a4 =5，1130e+004；a5 =8.3574e+004；a6 =9.3073e+0041号工厂提供天津、太原、石家庄、济南的产品；2号工厂提供郑州、上海、南京、合肥的产品；3号工厂提供西安、武汉的产品；4号工厂提供重庆、贵阳的产品；5号工厂提供长沙、广州、南宁的产品；6号工厂提供南昌、杭州、福州的产品。问题三：模型III首先根据网上资料得到18个城市的经纬度坐标如下：天津 （39.095963,117.180176） 太原 （37.872685,112.549438）石家庄（38.056742,114.510498） 济南 （36.668419,116.982422）郑州 （34.759666,113.620605） 西安 （34.252676,108.94043）上海 （31.231592,121.464844） 南京 （32.045333,118.806152）合肥 （31.821565,117.224121） 武汉 （30.581179,114.32373）重庆 （29.573457,106.54541） 杭州 （30.278044,120.168457）长沙 （28.207609,112.939453） 南昌 （28.690588,115.839844）贵阳 （26.667096,106.633301） 福州 （26.086388,119.300537）广州 （23.14036,113.269043） 南宁 （22.816694,108.363647）；根据重心法写程序得到以下结果（源程序见附录）：即可得到6个新厂址可设在距离重心81.4095处、 77.5951处 、 75.1697处、 81.6040处、 75.3840处、 77.4427处。从这些客观实际问题讨论使成本最小的基础上获利最大，这样更科学、更具有代表性，更有利于公司选择正确的厂址。六、 模型的分析与检验模型I我们对未来一年需求量的预测，因为短期预测随机性太大，不能做出精确预测，由拟合图可看出，三次多项拟合出的预测数据基本符合预测要求。模型II根据建立的模型，我们得出各厂生产量，与各厂人数。即：各个厂五月份的生产量：x1=190640000；x2=187830000；x3=89980000；x4=93260000；x5=128370000；x6=142960000；符合工厂产量的约束条件。需要加班的工厂：t2=36；t4=36；各个工厂的人数规模：a1 =1.2411e+005；a2 =1.02977e+005；a3 =5.8581e+004；a4 =5.1130e+004；a5 =8.3574e+004；a6 =9.3073e+004由于对最大生产容量不考虑的原因，所以各厂生产x1至x6,结果合理可行。各个工厂的人数最大的有12万多人，其余均为几万人，而二厂和四厂即使正常工作时间的生产量供应不足，也可以通过加班增产。模型III 通过第二问各厂对指定城市的供应，我们采用0-1变量法确定制定城市，则使得重新建厂时，在需求点（各个城市的坐标）、运输量及线性运输费率不变的前提下，我们将运输成本作为唯一决策因素。利用重心法可得到最优解，将求解得到的值带入计算，运输成本有明显降低。七、 模型的推广与评价模型评价：利用数据拟合可以很好处理一些比较繁杂的数据，并且起能很好反映出结果，因此问题一采用此解法比较合理。优化模型便是求最优解，而第二题便是要求最小成本的基础下来设定工厂规模，而0-1变量也能使问题简便化，使问题不会那么繁琐。模型推广：本题当中所涉及的采购模型考虑到的影响因素、约束条件还不够多。因此结果可能有些偏差。故对大型工厂相对误差较大，但适合小型的工厂预测和小型工厂规模优化。参考文献：1姜启源，谢金鑫，叶俊，数学模型第四版，高等教育出版社。2胡本超，数学建模仓库选址问题.3百度地图。附录一：必做题一一 模型假设(1) 由托里拆利定律：从小孔流出的液体的流速与水面高度的平方根成正比。由题目可知，水塔的最低水位为：270.3024=8.1648m最高水位为：270.3024=10.7352m又因为sqrt(10.7352/8.1648)1.1467，所以可忽略水位对流速的影响。（2） 将流量定义成单位时间流出水的高度,则可以将其看成时间的连续函数，又因为水塔截面积为S=(57*0.3048)2*pi/4=237.8(m2).因此得到流量后乘以s即为水位高度。二 建立模型（1）求任意时间内的水流量做出水位随时间的散点图如下：【1】拟合出时间水位函数从题目所给的表一数据，结合散点图可知，一天有两次供水时段和三次未供水时段，选取第一，二次未供水时段的数据作多项式拟合，便可得到水位函数。(这里拟合多项式的次数选为3，因为第三段供水时间的数据太少，拟合多项式的次数过高便不能得到很好的拟合)设t,h分别为已输入的时刻和水位测量记录，则根据第一段未供水时段，利用MATLAB拟合出水位函数：t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66；h=9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c1=polyfit(t(1:10),h(1:10),3)；%采用三次拟合的水位函数;tp1=0:0.1:8.9;x1=polyval(c1,tp1)；%求出函数值;plot(tp1,x1)；拟合出的图像如下：同理，拟合出第二段未供水时段水位函数： c2=polyfit(t(11:21),h(11:21),3); tp2=10.9:0.1:20.9; x2=-polyval(c2,tp2);plot(tp2,x2)；拟合出图像如下：【2】确定流量时间函数因为流量函数是水位函数的导函数，因此可以对对之前的水位函数求导：c1=polyfit(t(1:10),h(1:10),3);c2=polyfit(t(11:21),h(11:21),3);a1=polyder(c1);%对第一段时间内的水位函数求导a2=polyder(c2);%对第二段时间内的水位函数求导tp1=0:0.01:8.97;tp2=10.95:0.01:20.84;x11=-polyval(a1,tp1);x111=-polyval(a1,0:0.01:8.97);1=100*trapz(tp1,x111);%用梯形法计算第一段未供水时间内的总用水量x12=-polyval(a1,7.93,8.97); x23=-polyval(a2,tp2);x112=-polyval(a2,10.95:0.01:20.84);2=100*trapz(tp2,x112);%用梯形法计算第二段未供水时间内的总用水量x24=-polyval(a2,10.95,12.03); x25=-polyval(a2,19.96,20.84); subplot(1,2,1)plot(tp1,x11*100) subplot(1,2,2)plot(tp2,x23*100)运行程序，得到图像如下： 这样就得到了第1,2未供水时段的-时间流量图。接下来，对于第1供水时段的流量用供水前后的流量进行拟合。为使流量函数在t=9和t=11连续，我们取t=7.93,8.97,10.95,12.03这四个点，用3次多项式进行拟合 D1=7.93,8.97,10.95,12.03;Ds=x12, x24; N1=polyfit(D1, Ds,3); N2=7.93:0.1:12.03; N=polyval(N1 ,N2); G1=8.97:0.01:10.95;Gs1=polyval(N1,8.97:0.01:10.95);Gsl1=100*trapz(G1,Gs1);%用梯形法积分出第一段供水时的总用水量plot(N2, 100*N)由此得到第一供水时段的时间-流量图。接下来，在第2供水时段之前取t=19.96,20.84两点的流量，用第三未供水时段的3个记录做差分得到两个流量数据21.62,18.48，然后用这4个数据做3次多项式拟合得到第2供水时段与第3供水时段的时间-流量图， 程序如下t3=19.96,20.84,t(22),t(23); ls3=x25*100,21.62,18.48; Ns=polyfit(t3,ls3,3); tp3=19.96:0.01:25.91; xx3=polyval(Ns, tp3);g2=20.84:0.01:24;Gs2=polyval(Ns, 20.84:0.01:24);Gsl2=trapz(g2,Gs2);%该语句计算第2供水时段的总用水量 plot(tp3,xx3)做出图像如下：这样，通过分时间段进行拟合，便描述了任何时刻（包括水泵正供水时）从水塔流出的水流量。【3】对总用水量的估计分别对供水的两个时段和不供水的两个时段使用trapz语句进行梯形yom积分(流量对时间)，最后将各时段所得结果求和便得到一天的总用水高度。其值约为526.8435表2 时段第1未供水段第2未供水段第1供水段第2供水段全天用水用水高度145.65260.6346.6073.9635526.8435三 得出结论结合之前的模型分析，一天内的总用水量为：526.8435*237.8*10=1.2528e+006（升）必做题二一 模型假设（1）销售量的变化虽然是离散的，但对于大量的销售而言，可设销售量的变化随售价的增加而线性递减。（2）销售增长因子虽然也是离散的，但当广告费逐渐增加时，可设销售增长因子也是连续变化的。二 模型建立【1】售价与预期销售量的模型根据条件（表1）描出散点图， 程序如下： x=2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0；y=41,38,34,32,29,28,25,22,20；scatter(x,y,k*);画出的散点图如下：因此可以假设售价与预期销售量为线性关系，那么可得基本模型： 接下来可以直接用MATLAB里自带的CFT工具箱（曲线拟合工具箱）进行求解。其中，date x和date y分别取x,y。在fitting下的polynomial中选择 linear polynomial，运行结果如下：Linear model Poly1: f(x) = p1*x + p2Coefficients (with 95% confidence bounds): p1 = -5.133 (-5.573, -4.694) p2 = 50.42 (48.58, 52.27)Goodness of fit: SSE: 3.622 R-square: 0.9909 Adjusted R-square: 0.9896 RMSE: 0.7193则解得售价与预期销售量的线性回归方程为：=50.42-5.133x【2】广告费与销售增长因子的模型根据条件（表2）描出散点图，程序如下：m=0,1,2,3,4,5,6,7; k=1.00,1.40,1.70,1.85,1.95,2.00,1.95,1.80; scatter(m,k,k*)得到散点图：因此可以假设广告费与销售因子为非线性关系，则设其基本模型为：再次利用CFT工具箱，date x，date y分别改为m,k.,在fitting下的polynomial中选择madratic polynomial，则得到如下结果：Linear model Poly2: f(x) = p1*x2 + p2*x + p3Coefficients (with 95% confidence bounds): p1 = -0.04256 (-0.04701, -0.03811) p2 = 0.4092 (0.3768, 0.4416) p3 = 1.019 (0.9702, 1.067)Goodness of fit: SSE: 0.002515 R-square: 0.997 Adjusted R-square: 0.9957 RMSE: 0.02243因此，得到广告费与销售因子的回归方程为：k=-0.04256*m2+0.4092*m+1.019【3】预期利润的模型 为了最大预期利润，建立预期利润的目标函数为其中，限制条件为：要解目标函数max Z，等价于求min (-Z),接下来，利用MATLAB求解min (-Z)程序如下：f=104*x(2)-(x(1)-2)*(50.42-5.133*x(1)*103*(-0.0426*x(2).2+0.4092*x(2)+1.0187)f =104*x(2)-(x(1)-2)*(50.42-5.133*x(1)*103*(-0.0426*x(2).2+0.4092*x(2)+1.0187) x=fminsearch(f,4.5,3.5),f=eval(f)x = 5.9114 3.3082f = -1.1661e+005因此，得到的结果为：;x=5.9114;m=3.3082;三 得出结论：利润在预期的条件下获得最大利润为116610元此时最佳广告费用33082元售价应定为5.9114元。必做题三一 问题分析题目的要求是希望得到利润最大的方案。根据题目条件，只有在一个深度上挖了四块，才能在下一层再挖一块。反过来考虑，一旦挖取了第四层，则必须挖取之前所有的长方形块。根据题目中给出的条件，计算出其总收入为(1.5*5+0.75*5+2.0+1.0*2+0.5*2+0.25+8.0+9.0+6.0+0.5+27+6)*0.2=14.6000（万元），挖取费用为（万元），因为14.6-14.4=0.2（万元）依然是盈利的。因此，还必须考虑第四块挖取的情况。二 模型分析根据题目条件，挖取深层的块的前提是必须先挖取前一层相对应的块。例如，若挖取第二层的第一块长方形块，则必须先挖取第一层能够包含这一块的四块长方形。所以，应先确定底层长方形块的挖取是否能盈利，之后由底层的盈利情况反过来确定上层长方形块能否挖取。三 建立模型以各块含金属的百分数建立矩阵AN（N=1,2,3,4）则有：A1=1.5,1.5,1.5,0.75;1.5,2.0,1.5,0.75;1.0,1.0,0.75,0.5;0.75,0.75,0.5,0.25A1 = 1.5000 1.5000 1.5000 0.7500 1.5000 2.0000 1.5000 0.7500 1.0000 1.0000 0.7500 0.5000 0.7500 0.7500 0.5000 0.2500A2=4.0,4.0,2.0;3.0,3.0,3.0;2.0,2.0,0.5A2 = 4.0000 4.0000 2.0000 3.0000 3.0000 3.0000 2.0000 2.0000 0.5000A3=12.0,6.0;5.0,4.0A3 = 12 6 5 4A4=6.0A4 = 6接下来，可以将每一个AN矩阵的值转移为获得的利润情况，建立矩阵PN（N=1,2,3,4）表示每一个块的利润值（正值为盈利，负值为亏损）。在MATLAB中有：P1=A1*0.2-0.3P1 = 0.0000 0.0000 0.0000 -0.1500 0.0000 0.1000 0.0000 -0.1500 -0.1000 -0.1000 -0.1500 -0.2000 -0.1500 -0.1500 -0.2000 -0.2500P2=A2*0.2-0.6P2 = 0.2000 0.2000 -0.2000 0.0000 0.0000 0.0000 -0.2000 -0.2000 -0.5000P3=A3*0.2-0.8P3 = 1.6000 0.4000 0.2000 0P4=A4*0.2-1.0P4 = 0.2000从矩阵P1中可以观察到，在第一层中只有第6个元素为正值，因此在第一层只有挖第六块才够盈利。我们不难得到以下推论：1、在挖取多个层的长方形块的时候若要盈利，必须在底层块里选择利润为正值的长方形块（选0也可，只要保证不亏损），否则必将亏损。2、在考虑同一层的长方形块时，如要盈利，应选择挖取利润为正值的长方形块。若挖取亏损的长方形块，则总利润必将减少。接下来考虑挖取方案。（1）先从第四层开始考虑，如果挖了第四层（只能挖一块），那么就必须要将第一，二，三层全部挖开，从P4矩阵直接看出，总盈利额为0.2万元（2）接着从第三层开始挖取。先考虑挖取一块的情况。例如挖取第三层第1块，则必须挖取第二层的第1、2、4、5块以及第一层的第1、2、3、5、6、7、9、10和11块，由PN矩阵，则可以计算出每种情况的盈亏值之和为：第四块：0-0.7-0.45-0.40-0.25=-1.8第三块：0.2+0-0.4-0.45=-0.65第二块：0.4+0-0.65=-0.25第一块：1.6+0.4-0.25=1.75将所得数据填入表格：表1 仅挖取一块时第三层各块的盈亏值长方形块1234盈亏值（万元）1.75-0.25-0.65-1.8由表1容易得到，只有在挖取第1块的时候是盈利的，而其他情况下都是亏损的，因此，在第三层挖取2、3、4块时，根据之前的推论2，其总盈利值必然是减少的，因此不用考虑挖取更多块数的情况。（3）接下来，从第二层开始挖取，同样只挖取一块，采用如上述同样的方法，则可得到第二层各块长方形块的盈亏值如表2所示。表2 仅挖取一块时第二层各块长方形块的盈亏值长方形块123456789盈亏值（万元）0.30.3-0.5-0.1-0.15-0.9-0.5-0.8-1.3由表2得到，第二层只挖取一块的方案能获得的最大收入仅为0.3万元。那么如果从第二层挖取两块，从P2矩阵可以看到，只有第1、2块的盈利值大于0，其他均为负，因此，其总利润为0.2+0.2+0.1=0.5万元，所以第二层挖一块和两块的方案都不是最优；挖取三块或更多块的时候，由于其他均为负，那么可由推论2知其盈利值只能越来越小，因此无需考虑。（4）若第二层也都不挖取，只挖取第一层的长方形块，显然只有挖取第6块会盈利，其值为0.1万元，不是最优方案。综上所述，挖取第三层第1块长方形块的时候，收入是最大的，其值为1.75万元。四 得出结论为使该公司盈利额最大，挖取第一层的第1、2、3、5、6、7、9、10、11块，第二层的第1、2、4、5块和第三层的第1块长方形块时，收入与费用（即纯收入）之差最大，其最大值为1.75万元。附录二：Matlab程序及相关说明模型汇总和参数估计值因变量:天津需求量方程模型汇总参数估计值R 方Fdf1df2Sig.常数b1b2b3线性.048.45019.519360876.3645660.000对数.2402.83519.127304934.63656502.206倒数.4948.78119.016456300.983-223886.847二次.6597.74128.013172040.00092815.245-7262.937三次.88718.28237.001-8222.727241713.578-36976.5731650.758复合.1011.01219.341331885.2521.025幂.3404.64119.060279571.968.200S.62014.65819.00413.064-.746增长.1011.01219.34112.713.025指数.1011.01219.341331885.252.025Logistic.1011.01219.3413.013E-6.976自变量为 月份。确定使用三次多项式拟合后，使用MATLAB的CFT（曲线拟合工具箱），在“polynomial”选项里选择“cubic polynomial”（三次多项式拟合）拟合出的图像如下：（红线画出）结果显示：Linear model Poly3: f(x) = p1*x3 + p2*x2 + p3*x + p4Coefficients (with 95% confidence bounds): p1 = 1651 (610.1, 2691) p2 = -3.698e+004 (-5.591e+004, -1.804e+004) p3 = 2.417e+005 (1.417e+005, 3.417e+005) p4 = -8223 (-1.531e+005, 1.366e+005)Goodness of fit: SSE: 8.375e+009 R-square: 0.8868 Adjusted R-square: 0.8383 RMSE: 3.459e+004因此，我们可以得出回归方程为：Y1=1651*x3 -3.698e+004 *x2 + 2.417e+005*x-8223依据MATLAB设计程序，得到之后一年的预测值：a=12:1:24;p1=1651,-3.698e+004 ,2.417e+005, -8223;z1=polyval(p1,a)z1 = Columns 1 through 7 419985 511504 657841 868902 1154593 1524820 1989489Columns 8 through 13 2558506 3241777 4049208 4990705 6076174 7315521同理可得其他城市未来一年的预测值。程序验证：除了直接的出回归方程，我们还可以直接使用三次拟合函数进行验证：程序设计如下：x=1,2,3,4,5,6,7,8,9,10,11;y=206900290600487900480900476100440900415300424600382200 379000358800;p3=polyfit(x,y,3);xcurve=12:1:24;p3curve=polyval(p3,xcurve);plot(xcurve,p3curve,-)做出图像如下：求具体的预测值：table=p3curvetable = 1.0e+006 * Columns 1 through 9 0.4202 0.5117 0.6580 0.8691 1.1547 1.5249 1.9894 2.5583 3.2415 Columns 10 through 13 4.0488 4.9901 6.0753 7.3145可以看出，未来一年里的预测值与我们回归方程得到的预测结果大体一致，但从CFT能获得更加精确的解。模型缺陷：结合现实中的实际需求量考虑，使用三次多项式拟合的到的预测结果，随着月份的增长，其误差会显著增大。所以该模型可能只适用于在短期内需求量的预测。附录三：LINGO程序min=(x1*100/(192+t1)*8)*(1700+1.3*1700*t1/192)+(x2*100/(192+t2)*8)*(1540+1.3*1540*t1/192)+(x3*100/(192+t3)*8)*(1510+1.3*1510*t3/192)+(x4*100/(192+t4)*8)*(1600+1.3*1600*41/192)+(x5*100/(192+t5)*8)*(1640+1.3*1640*t5/192)+(x6*100/(192+t6)*8)*(1450+1.3*1450*t6/192)+47610000*(2.97*y11+5.59*y21+9.3*y31+15.2*y41+15.62*y51+14*y61)+47730000*(2.55*y12+5.5*y22+5.91*y32+12.38*y42+14*y52+13.37*y62)+48120000*(1.16*y13+4.7*y23+6.95*y33+13.13*y43+14.