农作物施肥分析1 问题重述某研究所为了研究氮肥(N)、磷肥(P)、钾肥(K)三种肥料对于土豆和生菜的作用，分别对土豆和生菜两种作物做了三组实验，实验中将每种肥料的施肥量分为10个水平，在考察其中一种肥料的施肥量对土豆或生菜的产量的影响时，其他两种肥料的施肥量均在第七个水平上，实验数据如下面表格所示，其中ha 表示公顷,t 表示吨,kg 表示千克，试建立反映施肥量与产量关系的模型，并从应用价值和如何改进等方面作出评价. 施肥量与产量关系的实验数据： 土豆: 氮肥 磷肥钾肥生菜： 氮肥 磷肥钾肥施肥量 (/)kg ha 产量 (/)t ha 0 15．18 34 21．36 67 25．72 101 32．29 135 34．03 202 39．45 259 43．15 336 43．46 404 40．83 471 30．75 施肥量 (/)kg ha 产量 (/)t ha0 33．46 24 32．47 49 36．06 73 37．96 98 41．04 147 40．09196 41．26245 42．17 294 40．36 342 42．73施肥量 (/)kg ha 产量(/)t ha 0 18．98 47 27．35 93 34．86 140 38．52 186 38．44 279 37．73 372 38．43 465 43．87 558 42．77 65146．22施肥量 (/)kg ha 产量 (/)t ha 0 6．39 49 9．48 98 12．46 147 14．33 196 17．10 294 21．94 391 22．64 489 21．34 587 22．07 685 24．53 施肥量 (/)kg ha 产量 (/)t ha0 11．02 28 12．70 56 14．56 84 16．27 112 17．75 168 22．59224 21．63 280 19．34 336 16．12392 14．11施肥量 (/)kg ha 产量 (/)t ha 0 15．75 47 16．76 93 16．89 140 16．24 186 17．56 279 19．20 372 17．97 465 15．84 558 20．11 65119．402模型假设：1、忽略土壤现有肥力对农作物生长的影响.2、假设在所有作物的生长环境的阳光，空气、水、温度变化都相同.3、假设施肥时肥料的质量准确.3问题分析：我们研究的是不同肥料的施肥量对农作物产量的影响，就是要求出施肥量与作物产量之间的关系，目前我们仅有单一施肥量与作物产量的数据，我们需要用线性回归的方法拟合出每一种肥料对作物产量影响的方程，然后根据第七个水平的共同施肥量建立多种肥料对一种作物产量影响的方程，从而找到最佳施肥量使得作物的产量最大.4模型的建立和求解4.1氮肥施肥量改变时对土豆产量的影响由现有数据建立氮肥施肥量对土豆产量的影响散点图如下:图 1 氮肥施肥量与土豆产量关系散点图由图 1 氮肥施肥量与土豆产量关系散点图我们发现氮肥的施肥量与土豆产量大致成二次关系，故我们建立如下模型：111212122212111p p n n np x x x xx x X x x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦，12n y y Y y ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦ ，01p b b B b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦ . 因1111112121121121222111111121221111111111nni ip i i p n nnn p i i i ip T i i i p pnp n n np nnnipip i ip i i i n x x x x x x x x x x x x xx x X X x x x x x x xxxx ========⎡⎤⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦∑∑∑∑∑∑∑∑1111211211121111n i i nn i i T i p p np n n ip ii y y x x x y x y X Y x x x y x y ===⎡⎤⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦∑∑∑ 于是（4.5）式可写成T T X XB X Y = （4.5）’ 这就是正规方程组的矩阵形式. 在()4.5’式两边左乘T X X 的逆矩阵()1TX X -（设()1TX X - 存在）得到()4.5’ 的解.011ˆˆˆ()ˆT T p b b B X X X Y b -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦， 这就是我们需要求的01(,,)T p b b b 的最大似然估计. 我们取011ˆˆˆˆp pb b x b x y ++= 作为12011(,,,)p p p x x x b b x b x μ=+++的估计. 则可以得到p 元线性回归方程方程：011ˆˆˆˆp py b b x b x =++经计算得正规方程组的解为0112ˆ14.74164ˆˆ()0.00034ˆ0.19715T T b B b X X X Y b -⎡⎤⎡⎤⎢⎥⎢⎥===-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦则求的原假设模型为20.000340.1971514.74164N NN y x x =-++ 用2R 来检验模型的拟合度，将所得模型与现有数据的散点绘制在同一张图上得氮肥施肥量与土豆产量拟合图形如下：计算得20.986R = 拟合度良好，可以较好的反映氮肥施肥量与土豆产量的关系.4.2 磷肥施肥量改变时对土豆产量的影响由现有数据画出磷肥施肥量与土豆产量关系的散点图：图 2 磷肥施肥量与土豆产量关系散点图y = -0.0003x 2+ 0.1971x + 14.742R² = 0.9863010203040500100200300400500土豆产量氮肥施肥量由图 2 磷肥施肥量与土豆产量关系散点图可以看出磷肥施肥量与土豆产量大致成二次关系，故可建立磷肥施肥量与土豆产量关系的模型如下：2012P P P y b b x b x =++为简化计算，计算出2p x 令之为2P x ，则上述假设模型可化为：2120p p p y b x b x b =++经计算得正规方程组的解为0112ˆ32.916ˆˆ()0.000138ˆ0.071859T T b B b X X X Y b -⎡⎤⎡⎤⎢⎥⎢⎥===-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦ 则假设模型可解得：20.0001380.07185932.916P p p y x x =-++用2R 来检验模型的拟合度，将所得模型与现有数据的散点绘制在同一张图上得磷肥施肥量与土豆产量拟合图形如下：经计算得20.864R = ，模型拟合程度良好，可以较好的反映磷肥施肥量与土豆产量之间的关系.4.3 钾肥施肥量改变时对土豆产量的影响由现有钾肥施肥量与土豆产量的数据画出如下钾肥施肥量与土豆产量关系散点图y = -0.0001x 2+ 0.0719x + 32.916R² = 0.8645010203040500100200300400土豆产量磷肥施肥量图 3 钾肥施肥量与土豆产量关系散点图由图 3 钾肥施肥量与土豆产量关系散点图可以初步看出钾肥的施肥量与土豆的产量大致成对数关系，我们因此假设土豆的产量与钾肥的施肥量模型如下：13ln K k y c x ε=+为简化计算令ln ln k k x x = 先将全部的ln k x 计算出来用lnk x 来代换，则假设模型可化为：ln k k y a bx ε=++求一元线性回归的解法如下： （1）模型假设()20,y a bx N εεσ=++⎧⎪⎨⎪⎩（2）参数估计i i i y a bx ε=++()i i i y a bx ε-+= ，()20,iN εσ()()2211,nni i i i i y a bx Q a b ε===-+=⎡⎤⎣⎦∑∑()(),0,0Q a b aQ a b b∂⎧=⎪⎪∂⎨∂⎪=⎪∂⎩112111n ni i i i n n ni i i i i i i na b x y x a x b x y=====⎧⎛⎫+=⎪ ⎪⎪⎝⎭⎨⎛⎫⎛⎫⎪+= ⎪ ⎪⎪⎝⎭⎝⎭⎩∑∑∑∑∑()221111211ninn n i ni i i ni i i ii i i xnn x x n x x x x ======⎛⎫=-=-≠ ⎪⎝⎭∑∑∑∑∑∑故方程有唯一解：()()()121niii ni i x x y y b x xa y bx==⎧--⎪⎪=⎪⎨-⎪⎪=-⎪⎩∑∑（3）得到回归方程为：y a bx =+经计算得ˆ 5.811ˆ7.021b a ⎧=⎪⎨=⎪⎩ 则假设模型为：ln 5.8117.021k k y x =+因此可以得出原假设模型为：5.811ln 7.021k k y x =+用2R 检验模型拟合度，将所得模型与现有数据的散点绘制在同一张图上得钾肥施肥量与土豆产量拟合图形如下：经计算得模型20.865R = ，因此模型拟合度良好，可以较好的反映钾肥施肥量与土豆产量之间的关系.y = -0.3137x 2+ 5.8579x + 16.576R² = 0.892101020304050024681012土豆产量钾肥施肥量4.4 氮、磷、钾肥施肥量共同改变时对土豆产量的影响氮肥、磷肥、钾肥三种肥料对土豆的产量的影响模型组如下：220.000340.1971514.741640.0001380.07185932.9165.811ln 7.021N N N P P P KK y x x y x x y x ⎧=-++⎪=-++⎨⎪=+⎩ 由以上模型组可以假设土豆产量（y ）与氮肥、磷肥、钾肥三种肥料的施肥量的关系模型如下：2222ln N N N P P p k N P y a x a x b x b x c x ε=+++++将30组数据汇总在一起，在没有其他肥料的施肥量数据的地方填上对应的肥料的第七水平的施肥量，分别计算出2N x 、2p x 、ln k x ，然后做多元线性回归分析，经计算得：220.0003130.18380.0001470.07536.066729.7174N NP Pa ab bc ε=-⎧⎪=⎪⎪=-⎨=⎪⎪=⎪=-⎩ 因此解出假设模型为：220.0003130.18380.0001470.0753 6.0667ln 29.7174N N P P K y x x x x x =-+-++-计算2R 来检验模型拟合度，得20.953R =，接近1拟合度很好，可以很好的反映氮肥、磷肥、钾肥的施肥量对土豆产量的影响.计算出最大产量时的氮磷肥三种肥料的施肥量如下表（MATLAB 代码见附录）：表 1 氮磷钾肥施肥量对应土豆产量表氮肥施肥量（/kg ha ） 磷肥施肥量 （/kg ha ） 钾肥施肥量 （/kg ha ）土豆产量（t/ha ） 293.6104256.125065146.2116由表 1 氮磷钾肥施肥量对应土豆产量表得当氮肥施肥量为293.6104/kg ha ，磷肥施肥量256.1250/kg ha ，钾肥施肥量651/kg ha 时土豆的产量最大值为46.2116t/ha .4.5 氮肥施肥量改变时对生菜产量的影响由现有氮肥施肥量的改变对生菜产量的影响数据，建立氮肥施肥量与生菜产量关系的散点图如下：图 4 氮肥施肥量与生菜产量关系散点图由图 4 氮肥施肥量与生菜产量关系散点图可以看出生菜的产量与氮肥的施肥量大致成二次关系，故假设生菜的产量与氮肥施肥量之间的模型为：2120=b n n y x b x b ε+++生为了简化计算先计算出2n x 用2n x 代换，则原模型可化为：2120=b n n y x b x b ε+++生经计算得正规方程组的解为：011210.229ˆ()0.0002380.10132T T b B b X X X Y b -⎡⎤⎡⎤⎢⎥⎢⎥===-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦则原假设模型解得：2=-0.0002380.1013210.229nn y x x ++生 对模型进行2R 检验，将所得模型与现有数据的散点绘制在同一张图上得氮肥施肥量与生菜产量拟合图形如下：经计算得20.925R = 拟合度较好，可以较好的反映氮肥施肥量与生菜产量之间的关系.4.6 磷肥施肥量改变时对生菜产量的影响为了初步分析生菜产量在仅有磷肥改变时所受的影响，在现有数据的基础上建立磷肥施肥量与生菜产量关系的散点图如下：图 5 磷肥施肥量与生菜产量关系散点图由图 5 磷肥施肥量与生菜产量关系散点图初步可以看出磷肥的施肥量与生菜产量之间大致呈现二次关系，故建立如下模型：2012=b +b p p y x b x ε++生计算出2p x 用2p x 代换则原假设模型可以化为如下多元线性回归方程：y = -0.0002x 2+ 0.1013x + 10.229R² = 0.924905101520250100200300400500生菜产量氮肥施肥量2012=b +b p p y x b x ε++生经计算得正规方程组的解为0112 6.87566ˆ()0.0000550.060597T T b B b X X X Y b -⎡⎤⎡⎤⎢⎥⎢⎥===-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦ 则模型为：2=0.0000550.060597 6.87566p p y x x -++生 用2R 检验模型拟合度，将所得模型与现有数据的散点绘制在同一张图上得磷肥施肥量与生菜产量拟合图形如下：经计算的20.959R =，模型拟合度较好，故可以较好的反映磷肥的施肥量与生菜产量之间的关系.4.7 钾肥施肥量改变时对生菜产量的影响由已知钾肥施肥量与生菜产量关系数据，画出钾肥施肥量与生菜产量关系散点图如下：y = -5E-05x 2+ 0.0606x + 6.8757R² = 0.95860510152025300200400600800生菜产量磷肥施肥量图 6 钾肥施肥量与生菜产量关系散点图由图 6 钾肥施肥量与生菜产量关系散点图可以看出第八组数据与其他数据存在较大差异，为了保证数据的准确性，对本次实验数据进行残差检验，残差分析图如下：图7 钾肥施肥量对生菜产量影响的残差分析图从残差图可以看出数据的残差离零点的远近，除了第八组数据的95%残差置信区间不包含零点，其他组数据的残差95%置信区间均包含零点，则第八组数据可视为异常点剔除.在剔除异常点后画出钾肥施肥量与生菜产量关系的修正散点图如下：图 8 钾肥施肥量与生菜产量关系修止散点图由图 8 钾肥施肥量与生菜产量关系修止散点图可以看出钾肥的施肥量与生菜的产量之间大致呈现二次关系，故建立模型如下：2012=b +b k k y x b x ε++生同理将2k x 计算出来，用2k x 代换，进行多元线性回归分析.经计算得正规方程组的解为011215.836ˆ()0.0000070.11T T b B b X X X Y b -⎡⎤⎡⎤⎢⎥⎢⎥===-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦ 则原假设模型为：2=-0.000007+0.01115.836k k y x x +生用2R 检验模型的拟合度，将所得模型与现有数据的散点绘制在同一张图上得钾肥施肥量与生菜产量拟合图形如下：计算得20.830R =模型的拟合度良好，可以较好的反映钾肥施肥量与生菜产y = -7E-06x 2+ 0.0106x + 15.836R² = 0.830105101520250200400600800生菜产量钾肥施肥量量之间的关系.4.8 氮、磷、钾肥施肥量共同改变时对生菜产量的影响由4.5，4.6，4.7得氮肥、磷肥、钾肥三种肥料对生菜的产量的影响模型组如下：22=-0.0002380.1013210.229=0.0000550.060597 6.87566=0.00594916.227n n N p p p k k y x x y x x y x ⎧++⎪-++⎨⎪+⎩生生生 根据上述模型组，我们可以建立生菜产量（y ）与氮肥（n x ）、磷肥（p x ）、钾肥（k x ）之间的模型：2220123456n n p p k k y b b x b x b x b x b x b x ε=+++++++经计算得正规方程组的解为： 4.8154250.0002320.098706ˆ0.0000380.0487480.0000300.027461B -⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎣⎦则可以解出模型为:2220.0002320.0987060.0000380.0487480.0000300.027461 4.815425n n p p k k y x x x x x x =-+-+-+-用2R 检验模型拟合度，计算得20.866R =接近于1，拟合度较高，可以较好的反映氮、磷、钾三种肥料同时施肥时对生菜产量的影响.经计算得出上述模型的最大产量时的氮磷钾的施肥量如下表（代码见附录）：表 2 氮磷钾肥施肥量对应生菜产量表氮肥施肥量（/kg ha ） 磷肥施肥量 （/kg ha ） 钾肥施肥量 （/kg ha ） 生菜产量 （t/ha ） 212.8283641.4227457.685523.5586由表 2 氮磷钾肥施肥量对应生菜产量表可知当氮肥施肥量为213/kg ha ，磷肥施肥量641/kg ha ，钾肥施肥量458/kg ha 时生菜的产量最大值为23.56t/ha（为了方便施肥，对小数点数据后进行四舍五入）.5 模型的评价与改进我们种植农作物时，不能单一的只考虑高产量，同时还需考虑作物的售价和肥料价格，以达到收益的最大化.目前各种作物价格以及肥料价格如下表表 3 作物及肥料价目表种类 价格（元/kg ）土豆 3.6 生菜 6.0 氮肥 1.6 磷肥 0.7 钾肥3.5用1w 、2w 分别表示土豆、生菜的利润，t p 、s p 分别表示土豆、生菜的价格，n z 、p z 、k z 分别表示氮、磷、钾三种肥料的单价.5.1 土豆最佳施肥根据上述价格及方程则可以建立土豆的利润模型为:221(0.0003130.18380.0001470.0753 6.0667ln 29.7174)N N P P K tn n p p k kw x x x x x p x z x z x z =-+-++----带入单价得上述模型为（MATLAB 代码见附录）：2213600(0.0003130.18380.0001470.0753 6.0667ln 29.7174)1.60.7 3.5N N P P K n p kw x x x x x x x x =-+-++----经计算得：表 4 土豆最佳施肥计算结果表(kg/ha)(kg/ha)(kg/ha)利润/ha)(元 292.9003 255.4611 651.0000163430由表 4 土豆最佳施肥计算结果表可知道当土豆的施肥量为氮肥293kg/ha ，磷肥255kg/ha ，钾肥651kg/ha ，可获得最大利润163430/ha 元.5.2 生菜最佳施肥根据上述价格及方程则可以建立生菜的利润模型为:22226000(0.0002320.0987060.0000380.0487480.0000300.027461 4.815425) 1.60.7 3.5n n p p kk n p kw x x x x x x x x x =-+-+-+----经计算得（MATLAB 代码见附录）：附录土豆计算最大产量的MATLAB 代码 M 文件：function f=fun3(x)f=-(-0.000313*x(1)*x(1)+0.1838*x(1)-0.000147*x(2)*x （2）+0.0753*x(2)+6.0667*log(x(3))-29.7174);主程序：x0=[1;1;1];A=[ 1 0 0;0 1 0;0 0 1]; b=[471;342;651]; Aeq=[];beq=[];VLB=[0;0;0]; VUB=[];x0=[0;0;0];[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,VLB,VUB)结果： x =293.6104 256.1250 651.0000 fval = -46.2116即293.6104256.125065146.2116n p k x x x y =⎧⎪=⎪⎨=⎪⎪=⎩生菜计算最大产量的MATLAB 代码： M 文件：function f=fun3(x)f=-(-8.858450+0.098706*x(1)+0.048748*x(2)+0.027461*x(3)-0.000232*x(1)*x(1)-0.000038*x(2)*x(2)-0.000030*x(3)*x(3))主程序：x0=[1;1;1];A=[ 1 0 0;0 1 0;0 0 1]; b=[392;685;651];Aeq=[];beq=[];VLB=[0;0;400]; VUB=[];x0=[0;0;0];[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,VLB,VUB)结果:x =212.7283641.4227457.6855fval =-23.5586土豆最佳施肥计算代码：M文件function f=fun3(x)f=-(3600*(-0.000313*x(1)*x(1)+0.1838*x(1)-0.000147*x(2)*x(2)+0.0753*x (2)+6.0667*log(x(3))-29.7174)-1.6*x(1)-0.7*x(2)-3.5*x(3));主程序x0=[1;1;1];A=[ 1 0 0;0 1 0;0 0 1]; b=[471;342;651];Aeq=[];beq=[];VLB=[0;0;0]; VUB=[];x0=[0;0;0];[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,VLB,VUB)结果x =292.9003255.4611651.0000fval =-1.6343e+05。