1、 2016 年河南科技大学 模拟训练 三 承 诺 书 我们仔细阅读了数学建模 选拔 赛 的 规则 . 我们完全明白，在 做题期间 不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人研究、讨论与 选拔 题有关的问题。 我们知道，抄袭别人的成果是违反 选拔 规则的 , 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文 引用 处 和参考文献中明确列出。 我们郑重承诺，严格遵守 选拔 规则，以保证 选拔 的公正、公平性。如有违反 选拔 规则的行为，我们将 受到严肃处理 。 我们 选择的题号 是 （ 从 A/B/C 中选择一项填写 ） ： B
2、队员签名 ： 1. 段平涛 2. 屈德强 3. 何蕴琪 日期： 2016 年 8 月 27 日 2016 年河南科技大学 数学建模竞赛 选拔 编 号 专 用 页 评阅编号（评阅前进行编号）： 评阅记录（评阅时使用）： 评 阅 人 评 分 备 注 1 农作物用水量预测及智能灌溉方法 摘要 作 为一个农业大国，我国的农业水资源一直处于紧张的状态，大部分地区仍 然处于半干旱甚至干旱的状态，常年的地表灌溉用水不足导致过量抽取地下水，造成了渐渐无水可用的情况。本文旨在通过合理各个农作生产阶段的灌溉量 ，在限制总灌溉水量的条件下给出合理的最优产量灌溉量分配方案，对实际生产有着指导性的作用。 针对问题一，本
3、文先根据实际情况进行分析，选择出了符合实际情况的模型公式，再通过 MATLAB 拟合出二次函数，针对二次函数的精确性再做检验，得出了较为精确的作物需水量和作物单位面积产量 之间的函数关系。 针对问题二，本文根据四个主要阶段对总产量的影响来建立模型，使用了 BP神经网络 算法 和遗传 算法 来求出最优的作物水分生产函数模型。再将这两种算法下的 模型与加法模型和乘法模型作对比，综合评价最优的模型，发现遗传算法 相较 于 BP 神经网络较好 。原因是数据量较小，进而对比 其他两种模型， 发现 乘法模型的整体误差较小，模型优化程度较好。 针对问题三，本文使用动态 规划逆序算法求解 ，在存在阶段变量、决
4、策变量 和 状态变量等多重变量的情况下具有良好的解决效果。通过对比加法模型和乘法模型的结果，结合实际情况，认为乘法模型各个阶段敏感系数较为合理，能有效体现极端情况和一般情况下的非充分灌溉农作物产量的变化。 本文中的模型可以适用与半干旱和半湿润地区的农作物生产模型中，其中使用的遗传模型拟合的精度高，效果好，在非线性、多模型和多目标的优化问题中效果良好 ，用法方便。 关键词： MATLAB 遗传 算法 BP神经网络 动态规划逆序算法 2 一、 问题重述 随着水资源供需矛盾的日益加剧，发展节水型农业势在必行。智能灌溉应用先进的信息技术实施精确灌溉，以农作物实际需水量为依据，提高灌溉精确度，实施合理的
5、灌溉方法，进而提高水的利用率。 灌溉水利用系数是指在一次灌水期间被农作物利用的净水量与水源渠首处总引进水量的比值，它是衡量灌区从水源引水到田间作用吸收利用水的过程中水利用程度的一个重要指标，也是集中反映灌溉工程质量、灌溉技术水平和灌溉用水管理的一项综合指标，是评价 农业水资源 利用，指导节水灌溉和大中型灌区续建配套及节水改造健康发展的重要参考。据有关部门统计分析，我国灌区平均水 利用系数 仅为 0.45，节水仍有较大空间。 作物水分生产函数无论对节水灌溉的区域规划和系统评估，或是非充分灌溉的应用均具有深刻意义。非充分灌溉是指在灌溉水不能完全满足作物的生长发育全过程需水量的情况下，以作物水分生产
6、函数为理论依据，将有限的水科学合理（非足额）安排在对产量影响比较大、并能产生较高经济价值的需水临界期供水，从而建立合理的水量与产量关系及分配模式，在水分利用效率、产量、经济效益三方面寻求有效均衡，实现经济效益最大化。然而由于作物各生育阶段水分对产量影响的机理甚为复杂，目前尚难用严格准确的物理方程来描述。 1. 基于表 1 中作物 1 实测需水量与 产量的对应数据，建立该作物全生育期的水分生产函数的模型，即总产量与需水量之间的解析关系，给出详细过程及拟合效果。 表 1 全生命周期实测需水量与产量 处理号 需水量（ m3/亩） 产量（ kg/亩） 1 187.14 247.63 2 238.38
7、361.18 3 283.12 430.37 4 315.37 462.89 5 337.65 476.88 6 356.21 483.25 7 387.68 483.05 8 414.46 472.01 9 435.57 456.23 3 10 456.82 434.08 2. 作物的全生育期可以分为若干个生育阶段，以水稻为例，可以分为返青、分蘖、拔节孕穗、抽穗开花、乳熟、黄熟 6 个生育阶段。不同阶段灌溉水量不足均会对最终的产量有影响，表 2 为某地晚稻分蘖至乳熟各阶段受旱情况对产量影响的数据。基于表 2的数据，利用附录中材料，选取某类优化算法，寻求最优的作物水分生产函数模型，得到各阶段的
8、蒸发蒸腾量（可以理解为灌水量）与最终产量之间的关系。给出详细过程，并将所得结果与常见的机理模型（见附录）作对比。 表 2 某地晚稻蒸发蒸腾量及产量 处理编号 处理特征 1分蘖 (mm) 2拔节孕穗 (mm) 3抽穗开花 (mm) 4乳熟 (mm) 产量 (kg/hm2) 0 充分灌溉 148.1 111.8 124.7 89.4 7138.5 1 1轻旱 113.2 96.6 92.1 67.2 5757.0 2 1重旱 107.6 88.3 84.9 64.2 4576.5 3 2轻旱 133.9 91.0 106.9 70.3 6111.0 4 2重旱 132.1 77.9 93.9 65
9、.0 4555.5 5 3轻旱 128.2 99.4 85.3 78.7 5520.0 6 3重旱 129.7 92.5 71.9 69.4 5329.5 7 4轻旱 140.5 112.9 108.6 68.6 6345.0 8 4重旱 135.3 108.0 101.7 65.0 6040.5 9 12中旱 110.6 83.3 95.2 72.3 5076.0 10 23中旱 128.4 90.4 83.4 73.6 5442.0 11 34中旱 130.1 102.6 94.7 61.4 6130.5 3. 基于已有材料对于非充分灌溉的描述，结合表 2 晚稻蒸发蒸腾量与产量数据，选择合
10、适的作物水分生产函数（可以选择已得到的模型或者机理模型），建立该地区水稻的非充分灌溉制度的优化模型。在总供水量分别为充分灌溉总水量的 40%， 60%， 80%的情况下，解决供水量在各生育期（从分蘖到乳熟阶段）合理分配的问题，应包含阶段变量、决策变量、状态变量、系统方程、目标函数、初始条件及约束条件等描述（背景材料及提示见附录），最后给出求解方案。 4 二、 问题分析 2.1 问题一分析 根据表 1 的数据我们需要根据数据点的变化特征来拟合函数，函数的选取需要根据实际作物生产的规律和拟合后的检验来分析拟合的效果。 2.2 问题二分析 由表 2 可以发现主要影响产量的阶段是 分蘖 、 拔节孕穗
11、、 抽穗开花 和 乳熟 ，查阅资料可以发现作物的水分生产函数有四种，得到各阶段蒸发量和最后产量的关系后，并需要对每种模型做检验 ， 从而得出精确度较高的函数 。 2.3 问题三分析 作物在全生育期内对灌溉水量的敏感程度是不同的，而且前一个阶段的灌溉水量会对后一个阶段的作物灌溉水量上限产生影响。在影响蒸发蒸腾量的因素中，温度、光照和降水等自然因素 都对模型有不可预测的影响。本文在处理模型的时候通过合理的模型的建立来解决这样的因素变动。在非充分灌溉模型下比较不同的水分生产函数，通过计算得出的结果结合实际情况分析分配方案的合理程度。 三、 模型假设 1、不考虑气温、光照和降水等天气因素 2、不考虑田
12、垄高度 3、总灌溉水量等于作物的需水量 4、可以自由控制总灌溉水量 四、符号说明 变量 变量名 单位 X 作物需水量 2kg/ m Y 作物单位面积产量 2kg/ m R 相关系数 / 检验水平 / ,ij ixy 处理后的作为处理后样本的 输入和输出 / G 最佳的函数表达式 / c 为实常数 / D 并集 / N 群体规模 / 5 F(i) 第 i个 父代 个体的适应度函数 / if 每个个体 所对应 的目标函数值 / s(i) 选择概率 / u 均匀随机数 / sp 选择 操作的概率 / cp 交叉操作的概率 / j ETm 各阶段的潜在腾发量 2kg/ m mY 作物的潜在产量 2kg
13、/ m ij ETa 实际腾发量 2kg/ m aj Y 实际产量 2kg/ m Err 相对绝对误差和 / ()ifX 由遗传程序计算的第 i 组实际输出值 / iy 第 i 组理想 输出值 / i 阶段变量 / im 各个阶段初可供灌溉的水量 2kg/ m iq 各个阶段初可供灌溉的水量 2kg/ m 1F Jensen模型目标函数 / 2F Blank模型目标函数 / 五、模型的建立与求解 5.1.问题一求解 通过查找资料发现，作物需水量 x 和作物单位面积产量 Y 之间存在线性关系和二次抛物线关系 1，只有在一定范围内，作物单位面积产量才会随着总需水量线性增加。当 Y达到一定水平后，再
14、继续增加则要靠其他的农业措施。所以，线性关系一般只适用于灌溉水源不足、管理水平不高、农业资源未能得到充分发挥的中低产地区。随着水源条件的改善和管理水平的提高， Y 与 x 之间的关系会出现一个明显的界限值，当 x 小于这个界限值时， Y 随着 x 的增加而增加，开始增加的幅度较大，然后减少；当达到该界限值时，产量不再增加，随后 Y随 x的增加而减少。因此本文选用二次函数抛物线关系。2 使用 MATLAB处理数据，用最小二乘拟合法可以得到函数式为 2y = - 0 . 0 0 7 x + 5 . 1 8 6 3 x - 4 7 8 . 4 7 2 4 其中 相关系数 0.974R ，检验水平 0
15、.01 ，函数图像如下： 6 图 1：作物需水量 x和作物单位面积产量 Y散点图 图 2：拟合后作物需水量 x和作物单位面积产量 Y的二次函数 7 根据图 2， 可以发现散点图图像呈现如下特征 ： 1、 当 x X0=187.14,238.38,283.12,315.37,337.65,356.21,387.68,414.46,435.57,456.82; Y0=247.63,361.18,430.37,462.89,476.88,483.25,483.05,472.01,456.23,434.08; plot(X0,Y0,*) polyfit(X0,Y0,2) ans = -0.0070 5
16、.1863 -478.4724 X0=187.14,238.38,283.12,315.37,337.65,356.21,387.68,414.46,435.57,456.82; Y0=247.63,361.18,430.37,462.89,476.88,483.25,483.05,472.01,456.23,434.08; plot(X0,Y0,*) hold on x=187.14,238.38,283.12,315.37,337.65,356.21,387.68,414.46,435.57,456.82; y=-0.007*x.2+5.1863*x-478.4724 y = Column
17、s 1 through 8 246.9421 360.0626 428.7743 460.9234 474.6291 480.7406 480.0819 468.6019 Columns 9 through 10 452.4757 429.9416 plot(x,y) 问题二 遗传算法 function bin_gen,bits=encoding(min_var,max_var,scale_var,popsize) bits=ceil(log2(max_var-min_var)./scale_var); bin_gen=radint(popsize,sum(bits); functionvar
18、_gen,fitness=decoding(funname,bin-gen,bits,min_var,max_var) num_var=length(bits); popsize=size(bin_gen,1); scale_dec=(max_var-min_var)./(2.bits-1); bits=cumsum(bits); bits=0 bits; for i=:num_var bin_vari=bin_gen(:,bits(i)+1:bits(i+1); vari=sum(ones(popsize,1)*2.(size(bin_vari,2)-1:-1:0).*bin_vari,2)
19、.*scale_dec(i)+min_var(i); end var_gen=var1,:; for i=1:popsize fitness(i)=eval(funname,(vargen(i,:); end functionevo_gen,best_indiv,maxfitness=selection(old_gen,fitness) 19 popsize=length(fitness); max_fitness,index1=max(fitness); min_fitness,index2=min(fitness); best_indiv=old_gen(index,:); index=1
20、:popsize; index(index1)=0; index(index2)=0; index=nonzeros(index); evo_gen=old_gen(index,:); evo_fitness=fitness(index,:); evo_popsize=popsize-2 ps=evo_fitness/sum(evo_fitness); pscum=cumsum(ps); r=rand(1,evo_popsize); selected=sum(pscum*ones(1,evo_popsize)ones(evo_popsize,1)*r)+1; evo_gen=evo_gen(s
21、electec,:); function new_gen=crossover(old_gen,pc) nouse,mating=sort(rand(size(old_gen,1),1); mat_gen=old_gen(mating,:); pairs=size(mat_gen,1)/2; bits=size(mat_gen,2); cpairs=rand(pairs,1)pc; cpoints=randint(pairs,1,1,bits); cpoints=cpairs.*cpoints; for i=1:pairs new_gen(2*i-1,2*i,:)=mat_gen(2*i-1,2
22、*i,1:cpoints(i) mat_gen(2*i,2*i-1,cpoints(i)+1:bits); end 问题三动态规划程序 function p_opt,fval=dynprog(x,decisfun,subobjfun,transfun,objfun) k=length(x(1,:);x_isnan=isnan(x);t_vub=inf; t_vubm=inf*ones(size(x);f_opt=nan*ones(size(x); d_opt=f_opt; tmp1=find(x_isnan(:,k);tmp2=length(tmp1); for i=1:tmp2 u=feva
23、l(decisfun,k,x(i,k);tmp3=length(u); for j=1:tmp3 tmp=feval(subobjfun,k,x(tmp1(i),k),u(j); if tmp=t_vub f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp; end end end for ii=k-1:-1:1 tmp10=find(x_isnan(:,ii);tmp20=length(tmp10); for i=1:tmp20 u=feval(decisfun,ii,x(i,ii);tmp30=length(u); 20 for j=1:tmp30 tmp00
24、=feval(subobjfun,ii,x(tmp10(i),ii),u(j); tmp40=feval(transfun,ii,x(tmp10(i),ii),u(j); tmp50=x(:,ii+1)-tmp40; tmp60=find(tmp50=0); if isempty(tmp60) if nargin5 tmp00=tmp00+f_opt(tmp60(1),ii+1); else tmp00=feval(objfun,tmp00,f_opt(tmp60(1),ii+1); end if tmp00=t_vubm(i,ii) f_opt(i,ii)=tmp00;d_opt(i,ii)
25、=u(j); t_vubm(i,ii)=tmp00; end end end end end fval=f_opt(tmp1,1);%此语句是否应该为 fval=f_opt(tmp10,1)? fval=fval(find(isnan(fval),1);%此语句是否要去掉？ p_opt=;tmpx=;tmpd=;tmpf=; tmp0=find(x_isnan(:,1);tmp01=length(tmp0); for i=1:tmp01 tmpd(i)=d_opt(tmp0(i),1); tmpx(i)=x(tmp0(i),1); tmpf(i)=feval(subobjfun,1,tmpx(
26、i),tmpd(i); p_opt(k*(i-1)+1,1,2,3,4)=1,tmpx(i),tmpd(i),tmpf(i); for ii=2:k tmpx(i)=feval(transfun,ii-1,tmpx(i),tmpd(i); tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1=0); if isempty(tmp2) tmpd(i)=d_opt(tmp2(1),ii); end tmpf(i)=feval(subobjfun,ii,tmpx(i),tmpd(i); p_opt(k*(i-1)+ii,1 2 3 4)=ii,tmpx(i),tmpd(i),tmpf(i); end end
