动物群落的稳定发展

A题：动物群落的稳定发展

摘要：本文通过对某公园近两年内被运出的某种动物的年龄和性别的数据进行统计分析，并针对题目的四个问题分别建立了符合实际的数学模型，在模型的求解过程中，应用C语言进行编程调试，通过统计学软件SAS，数学软件MATLAB等计算工具，编写相应的程序，对建立的模型进行求解，得出了符合实际的结果。 问题一：我们假设新生幼仔的数量为x0，然后通过对各年龄阶段的存活率?、被运走的动物数量Bj以及该动物的总体数量的分析来建立该群落的动态变化模

d?(k)?t?60(k)60(k)型利用该群落近两年内被运走的各年龄阶段的个体数??xi,1??xi,0，

dti?1i?1量分布，用C语言编程计算，推测出当前该动物的年龄结构（具体结果见 7 页表一）。并利用MATLAB软件对得出的数据用图形表示，利用对比分析法，得到该动物群落的基本分布轨迹，最后用统计软件SAS对模型进行相关性的分析检验，求得相关系数R与P的值，验正了模型的稳定性。

问题二：由于现在采用注射避孕药的方法来维持该种群的稳定，而且已经没有个体被运走或被偷猎的情况，为此我们把该种群的稳定性转化为求目标函数

?1??1x0?[(1??2)c2?c3]（该种群每年的新生幼仔的数量减去该年死亡个体的数

量的差值）；另外从???xi?160(k)i,1k)（即年头的数量与该年年底的数量的差值）??xi(,0i?160当?趋于0时，即认为该群落的个体数量是稳定的，从而把问题的稳定性问题转

化为求单目标的最优化问题建立模型；利用MATLAB对模型进行求得，得出当不考虑不确定性因素影响时要注射药物的雌性动物数量为276头，而当考虑了双胞胎和被重复注射这两个不确定性因素影响后，得到要注射药物的雌性动物数量为352头，其中有110头是被重复注射的。

问题三：其大致模型与问题二相近，不同之处在于要考虑到被运走的动物的数量（b），即目标函数?应考虑上被运走的数量，即只是对问题二的模型进行扩充建立新的目标模型；?1??1x0?[(1??2)c2?c3]－b和???xi?160(k)i,1k)－b；利??xi(,0i?160用MATLAB对不同b值进行求解，从而得出相应的避孕措施。（具体结果见 19 页

表二）

问题四：我们引进了增量加速度的概念，利用c语言进行编程求解，然后用MATLAB软件对得到的数据进行线性回归分析，得到该群落在减少至M时重新壮大该动物群落能力的模型：M＝3.9010+0.0047D。最后应用统计软件SAS对模型进行稳定性分析。

关键字：存活率 年龄结构 新生幼仔数 稳定性 最优目标 增量加速度

一． 问题重述与提出

位于非洲某国的国家公园中栖息着近11000头某种野生动物。管理员要求有一个健康稳定的环境以便维持这个11000头该动物的稳定群落。过去的20年中，整个该动物群是通过一些偷猎枪杀以及转移到外地而稳定下来的。但是近年来，偷猎被禁止，而且每年要转移这些动物也比较困难，因此，要控制现在的数量就使用了一种避孕注射法。用这种方法注射一次可以使得一头成熟雌性动物在两年内不会受孕。

要探讨这种避孕注射法的实用性，我们需要完成以下问题：

1．探讨该动物年龄在2岁到60岁之间的合理的存活率的模型，推测这个动物群落的当前的年龄结构。

2．估计每年在该群落中有多少雌性动物要注射避孕药，可以式群落固定在11000头左右。这里不免有些不确定性，也要估计这种不确定性的影响。 3．假如每年转移50至300头此动物到别处，那么上面的避孕措施将可以有怎样的改变？

4．如果由于某种原因，突然使得注射避孕的方法不得不停止（例如由于一场灾难导致大量该动物的死亡），那时重新壮大该动物群的能力如何？

二． 基本假设与符号说明 （一）模型假设

1．该公园是非开放式的，它与外界不发生关系，从而构成独立的生物群落，

该动物群落不存在与其它动物种群的竞争，或虽有竞争，但其影响只局限于该动物群落的死亡率内。

2．种群是通过雌性个体的繁殖而增长的，所以用雌性个体数量的变化为主

要研究对象。

3．为了讨论的必要，我们把新生的幼儿的存活率定为75％，而其后的存活

率为95％，直到60岁为止。各年龄组的该动物经过一年后即进入高一级的年龄组，而龄超过60即认为全部死亡，退出该系统。

4．由于该公园加强了对该动物群落的保护，我们认为该动物没有再被偷猎

射杀。而该动物群落个体数量的减少只是因为自然死亡以及被运走。 5．假设同一年龄组的动物个体之间是同质的，我们只考虑其平均水平，不

讨论个别差异。

6．题设该动物在10～12岁开始怀孕，我们这里设定为11岁开始，经过22

个月（约两年）的怀孕期后生幼仔，即可认为该雌性动物在13～60岁的时间内可以生幼仔。

7．该群落的自然死亡是在生完幼仔后才发生的，产幼仔只发生在每年的年

初时段，而被运走只发生在年底时段。

（二）符号说明

?1：新生幼儿的存活率，其值为0.75；

?2：1～60岁个体的存活率，其值为0.95； ?3：双胞胎出生的几率，其值为0.0135；

?(k)(t)：该动物第k年时刻的数量；

k)

：该动物第k年初i龄动物的数量； xi(,0

k)

：该动物第k年初底i龄动物的数量 xi(,1

Bj：第j年被运走的动物的数量；

?(k)?0?：表示该动物第k年初时的总数量；

?1:表示每年没有注射避孕药的雌性动物生幼仔的几率，其值为

1； 3.5?2：表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔的几率，其

值为

1； 5.5表示被注射过避孕药但在两年内被重复注射的雌性动物生幼仔的几率，其值?3：为

1； 6.5c1：表示从13～60岁该动物的雌性个体的总数；

c2：表示从1～59岁该动物的个数总和；

c3：表示60岁该动物的个体总和；

y1：表示13～60岁雌性动物没有被注射避孕药部分的数量；

y2： 表示13～60岁雌性动物被注射过避孕药但在两年内不再被注射部分的数

量；

表示13～60岁雌性动物被注射过避孕药但在两年内被重复注射部分的数量； y3：

?1：表示每年出生幼仔的数量与该年个体死亡的数量的差值；

表示该种群每年的新生幼仔的数量减去该年死亡个体的数量与运走个体数量?2：

的和的差值；

?3：表示该动物群落在年底时的总数量与年初的数量加上被运走的个体数量b

的差值。

三． 问题分析与模型建立 问题一：

1．我们要研究该动物群落的稳定性问题，首先要根据存活率确定其当前的年龄结构。该动物的新生幼仔存活率较低，题设是70％到80％之间，为了讨论的需要，我们这里设定为75％。在1岁后的存活率比较高，在这里设为95％，直到60岁，而超过60岁则认为退出该系统。因此，我们先建立出该动物群落中年龄在2岁到60岁之间的合理的存活率的模型。

模型一：

?d?(k)?t?60(k)60(k)??xi,1??xi,0????????????(1)?dti?1i?1??60(k)(k)x??(0)??????????????????(2)??i,0?i?1?d?(k)?t??Bj???????????????????(3)dt??x(k?1)?x(k)......(i?0,1,..59)???????????(4)i,1?i,0?x(k)?0............(i?60)??????????????(5)?i,1?x(k)??x(k)......(i?0)??????????????(6)

1i,0?i,1k)k)?xi(,1??2xi(,0......(i?1,2,...,59)??????????(7)?60?(k)x?i,0?(k)?????????????????(8)?x0,0??1i?13?2式（1）表示该动物第k年增长的数量；

式（2）表示该动物第k年初时的总数量，可由已有的数据计算出?(k)(0)来； 式（3）表示该动物被运走的数量；

式（4）和（5）表示该动物第i龄到了年底全部转化为（i＋1）龄； 式（6）和（7）表示该动物各年龄段的变化； 式（8）表示该动物新生的幼仔数量。

2．通过对该公园近两年内从这个地区运出的该动物的年龄和性别的数据进行统计分析，并利用编程工具Turbo C 2.0对该模型进行编程计算（源程序及计算过程见附录1），可得到当前该动物群落的年龄结构，如下表所示：

表一 该动物的年龄结构统计表

前一年数量前一年运走前一年剩下前兩年数前兩年运走前两年剩下假设无运走今年数年龄（岁） （头） 数量（头） 数量（头） 量（头） 数量（头） 数量（头） 数量（头） 量（头） 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

792 594 0 564 0 535 0 508 3 482 4 457 7 434 20 412 9 391 15 371 9 352 22 334 3 317 23 301 5 285 13 270 21 256 0 243 22 230 14 218 5 207 13 196 10 186 0 786 594 564 535 505 478 450 414 403 376 362 330 331 294 296 272 249 256 221 216 213 194 186 186 806 604 0 573 20 544 21 516 13 490 12 465 13 441 22 418 14 397 40 377 14 358 26 340 13 322 14 305 27 289 3 274 14 260 12 246 20 233 25 221 17 209 14 198 10 188 0 806 604 553 523 503 478 452 419 404 357 363 332 327 308 278 286 250 248 226 208 204 195 188 188 746 559 531 504 478 454 431 409 388 368 349 331 314 298 283 268 254 241 228 216 205 194 184 174 800 600 569 540 512 486 461 437 415 394 374 355 337 320 303 287 272 258 245 232 220 208 197 187 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

176 13 167 30 158 14 150 12 142 0 134 20 127 6 120 3 113 5 107 8 101 12 95 10 90 3 85 7 80 14 75 10 71 16 67 21 63 13 59 10 56 12 53 6 50 3 47 6 44 9 41 13 38 10 36

3

163 137 144 138 142 114 121 117 108 99 89 85 87 78 66 65 55 46 50 49 44 47 47 41 35 28 28 33

178 2 169 3 160 4 151 4 143 3 135 2 128 3 121 13 114 16 108 13 102 10 96 10 91 12 86 16 81 12 76 10 72 12 68 19 64 13 60 24 56 17 53 16 50 25 47 12 44 45 41 23 38 34 36

13

176 166 156 147 140 133 125 108 98 95 92 86 79 70 69 66 60 49 51 36 39 37 25 35 -1 18 4 23

165 156 148 140 132 125 118 112 106 100 94 89 84 79 75 71 67 63 59 56 53 50 47 44 41 38 36 34

177 168 159 151 143 135 128 121 114 108 102 96 91 86 81 76 72 68 64 60 56 53 50 47 44 41 38 36

52 53 54 55 56 57 58 59 60 生幼仔的雌性数量 13—60岁雌性数量 总数量

34 32 30 28 26 24 22 20 18

6 21 15 4 13 10 32 14 0

28 11 15 24 13 14 -10 6 18

34 32 30 28 26 24 22 20 18

16 10 17 13 13 12 3 22 20

18 22 13 15 13 12 19 -2 -2

32 30 28 26 24 22 20 18 17

34 32 30 28 26 24 22 20 18

784

792

734 789

2735 11714

232 622

2503 11092

2772 11876

302 876

2470 11000

2569 11006

11808

注1：0岁表示新生幼仔。

注2：由于每个年龄段的数据均为推测值，而实际上运走的各年龄段的数量不一定全部与预测值相符，故表中“剩下数量”两组数据中出现负数可认为是独异点，不影响模型整体的准确性。

每年新生幼仔的数量（x0）减去生幼仔的雌性的数量（

?x/2i13603.5，由于雄性

与雌性的数量比接近1：1，我们可近似地认为13——60岁个体的雌雄数量相等），其差值即为双胞胎的数量，这个差值与生幼仔的雌性数量之比即为双胞胎的几率（?3）1.35％。由表中数据可得，

792?784806?792?0.0102，?0.0177，784792746?734800?789?0.0163，?0.0139，这些比例都基本上接近题设的双胞胎的735789几率0.0135，说明以上推测得出的数据是准确的。

利用Matlab软件对以上四组数据用图形表示，并进行比较，得到该动物群落的基本分布情况图（源程序见附录2），如下图所示

图1

分析该图，可以看出，这四组曲线的轨迹、分布情况基本相同。由于“预测当前的年龄结构情况（无运走）”一组数据没有减去被运走的个体数量，故其每个年龄层的数量都略多于前三组的数量，因此其曲线比前三组的曲线略高一点，

利用SAS软件对模型进行相关性的分析检验（源程序见附录3），得到如下结果：

图2

程序的分析及统计结论：

程序中的x1是前一年的该动物群落的年龄结构，x2是前两年该动物群落的年龄结构，x3是该动物群落没有被运出是的年龄结构，x4是预测的当前的该动物群落的年龄结构。过程中的PROC CORR是分析变量中两两变量之间的PEAROS简单相关的。

输出结果中的结果1是一些基本的描述统计量，结果2是两两变量之间的相关矩阵，其中包括相关系数和显著性检验的概率。由结果可知前一年的该动物群落的年龄结构（x1）与前两年的该动物群落的年龄结构（x2）的相关系数R=0.99744,P=0.0001<0.01，所以前一年的该动物群落的年龄结构（x1）与前两年的该动物群落的年龄结构（x2）之间存在着极显著的正相关；前一年的该动物群落的年龄结构（x1）与没有运走是的该动物群落的年龄结构（x3）的相关系数R=0.99896,P=0.0001<0.01，所以前一年的该动物群落的年龄结构（x1）与没有运走时的该动物群落的年龄结构（x3）之间存在着极显著的正相关；同理可知x1与x4的相关系数R=0.99896,P=0.0001<0.01;x2与x3的相关系数R=0.99846,P=0.0001<0.01；x2与x4的相关系数R=0.99854,P=0.0001<0.01；x3

与x4的相关系数R=0.99999,P=0.0001<0.01 。

由以上的分析可知，x1,x2,x3,x4之间的相关系数接近1，可见模型一的稳定行很强，而且由公式推出的前一两年的数据与该公园已有的数据基本相符合，可见模型是很优的。

问题二：

由于目前该动物已经很少被移出或移入，而且偷猎枪杀的情况微乎其微，所以暂时不列入考虑范围内。因此对该动物群落若不采用人工手段控制，则其在一定时间范围内会大幅度增加，从而破坏该种群的动态平衡。为了保持该种群的平衡，而又不必每年运走一定数量动物，现在使用一种避孕注射法，可使该动物群落的数量固定在一定范围内，用这种方法注射一次可以使得一头成熟雌性动物在两年内不会受孕，但不会引起其它附加的反应。我们所要做的就是估计出每年在该群落中要注射避孕药的雌性动物的数量，并且要考虑到各种不确定性因素的影响。

为了分析的方便，我们先建立初步模型，该模型暂时不考虑注射避孕药所产生的不确定性因素的影响，即不考虑两年内被重复注射的雌性数量及双胞胎的几率。在这里我们只认为新生幼仔的数量由两部分组成，一部分为没注射过避孕药的雌性个体所生，另一部分为被注射过避孕药的雌性个体所生。另外，由于已经没有个体被运走或被偷猎的情况，为此我们把该种群的稳定性转化为求目标函数，当?趋于0?（该种群每年的新生幼仔的数量减去该年死亡个体的数量的差值）

时，即认为该群落的个体数量是稳定的，从而把问题的稳定性问题转化为求单目标的最优化问题。从而建立模型如下:

模型二：

minz???1y1+?2y2?x0=0?? ? y1?y2?c1????x?[(1??)c?c]10223?

其中：?1＝

1，表示每年没有注射避孕药的雌性动物生幼仔的几率。 3.51?2?，表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔

5.5的几率；

这里设c1＝2596，表示从13～60岁该动物的雌性个体的总数为2596； 设c2＝10243，表示从1～59岁该动物的个数总和为10243； 设c3＝17，表示60岁该动物的个体总和为17；

?：表示该种群每年的新生幼仔的数量与该年死亡个体的数量总和的差

值。

把已知的数据代入上述模型，从而得到以下模型：

minz???0.2857y1+0.1818y2?x0=0?y1?y2?2569 ??0.75x0?529???利用Matlab软件进行编程（源程序见附录4），求得

y1?2293,y2?276,x0?705,??2.0286?10?16.

?16??2.0286?10?0，说明该模型是稳定的，即该动物群落的数量被其中

控制在一定的范围内。由y2?276可知，每年大约有276头雌性动物要注射避孕药，才能使该群落的数量保持在11000头左右。

模型三（模型二的改进）：

由于每年被注射的雌性动物数量一定，所以被注射过后两年有可能又被注射，则将其归入到“被注射过避孕药但在两年内不再被注射的雌性动物”内。因此，注射避孕药所产生的不确定性因素之一即为“被注射过避孕药但在两年避孕期内被重复注射的雌性动物”。另外不确定因素之二即双胞胎的几率问题，所以这里我们加入一个参数?3。

当第一年采用注射避孕药的方法时，是不会发生有雌性个体被重复注射的情

况的，故有以下模型：

minz????1y1+?2y2+?3(y1+y2)?x0=0?y1?y2?c1 ?????x?[(1??)c?c]10223?其中：

1，表示每年没有注射避孕药的雌性动物生幼仔的几率； 3.51?2?，表示被注射过避孕药但在两年内不再被注射的雌性动物生幼仔的

5.5?1＝

几率；

?3?几率；

1，表示被注射过避孕药但在两年内被重复注射的雌性动物生幼仔的6.5?1?0.75，表示新生幼仔的存活率； ?2＝0.95, 表示1～60岁动物的存活率；

?3?0.0135, 表示新出生的幼仔中双胞胎的概率。

minz???0.2857y1+0.1818y2+0.0135(0.2857y1+0.1818y2)?x0=0?y1?y2?2569 ????0.75x0?[(1?0.95)10243?17]?

利用Matlab软件求解（源程序见附录5），求得

y1?2202,y2?367,x0?705,??2.4791?10?16

?16??2.4791?10?0，说明该模型是稳定的，即该动物群落的数量被其中

控制在一定的范围内。由y2?367可知，每年大约有367头雌性动物要注射避孕

药，才能使该群落的数量保持在11000头左右。由于考虑双胞胎的机率，故必须增加注射的数量。

当注射避孕药一年后再次注射时，就会有某些数量的雌性个体被重复注射的情况出现，但这部分一定比前一年注射的雌性个体的数量少，故建立以下模型

模型四：

minz????1y1+?2y2+?3y3+?3(?1y1+?2y2+?3y3)?x0=0?y3???? y?y?y?c1231?????1x0?[(1??2)c2?c3]?代入已知数据，得到以下目标函数模型：

minz???0.2896y1+0.1843y2+0.1698y3?x0=0?y3????

y?y?y?2569123????0.75x0?529?其中?是前一年被注射避孕药的雌性个体的数量。利用Matlab软件编写程序求解，通过对?的不同取值进行调试，（源程序见附录6）求得当?＝242时，模型稳定，此时

y1?2217,y2?242,y3?110,x0?705,??1.6409?10?12

其中x0?705与模型三的解相同，说明新生幼仔的数量比较稳定，而每年注

?12射避孕药的雌性动物则减少到352头，其中有110头是在注射后一年又被重复注射的数量，而??1.6409?10?0则说明了该动物群落的个体数量是稳

定的，被控制在一定范围内。由此可知当每年注射352头时，就能把该动物群落的数量控制在11000头范围内。

问题三:

假如每年转移50至300头此动物到别处，其大致模型与问题二相近，不同之处在于要考虑到被运走的动物的数量（b），即?应表示该种群每年的新生幼仔的数量减去该年死亡个体的数量与运走个体数量的和的差值。相应的避孕措施将改变如下： 模型五：

minz????1y1+?2y2+?3y3+?3(?1y1+?2y2+?3y3)?x0=0?y3????

y?y?y?c1231?????1x0?[(1??2)c2?c3?b]?其中b表示该年被运走的动物的数量（50～300）； 代入已有数据，并把方程标准化，得到以下模型：

minz???0.2896y1+0.1843y2+0.1698y3?x0=0?y3????y1?y2?y3?c1

????0.75x0?(0.05c2?c3)?b?其中b的值在50～300之间变化，我们以50头为间距（即b的值分为50，100，150，200，250，300六种情况），利用软件MATLAB6.0计算得到每年转移从50至300头到别处时所应采取的避孕措施（源程序见附录7），当取b＝50时，得到以下输出结果： ：

由输出结果的提示可知，在100000次迭代后仍然无法找到最优值；由exitflag=-1可知，以上模型得到的数据是发散的，从而说明上述模型本身并不存在问题，只是计算机无法在有限的迭代次数内找到最优值。为此解决该问题，我们对以上模型进行修改，把目标函数?3定义为该动物群落在年底时的总数量与年初的数量加上被运走的数量b的差值，从而得到以下模型： 模型六：

minz??3?60??(k)60(k)3??xi,1???xi,0?b－－－－－－－－（1）i?1i?1?60??x(k)i,0?11000－－－－－－－－－－－（2）?i?1?x(k?1)(k)?i,0?xi,1......(i?0,1,..59)－－－－－－（3）??(?xk)i,1?0............(i?60)－－－－－－－－（4） ?x(k)(k)i,1??1xi,0......(i?0)－－－－－－－－（5）??x(k)x(k)i,1??2i,0......(i?1,2,...,59)－－－－－（6）??x(k)0,0??1y(k)1??2y(k)2－－－－－－－－－（7）?60?x(k)?i,0y(k)(k)?i?13??1?y2?2－－－－－－－－－（8）

其中

式（1）表示目标函数，其值趋于0；

式（2）表示该动物在第k年初的数量，其值取11000；

式（3）和（4）表示该动物第i龄到了年底全部转化为（i＋1）龄； 式（5）和（6）表示该动物各年龄段的变化； 式（7）表示该动物新生的幼仔数量；

式（8）表示该动物在13～60岁的雌性个体数量之和；

对b以10为间距，分别从50取到300头，利用MATLAB软件进行求解（源程序见附录8），得到以下表二：

被运走的数量（头） 被注射的数量（头） 50 60 70 80 90 100 110 120 130 140 370 361 351 342 332 322 313 304 294 285

被运走的数量（头） 被注射的数量（头） 150 160 170 180 190 200 210 220 230 240 275 266 256 247 237 228 218 209 199 190

被运走的数量（头） 被注射的数量（头） 250 260 270 280 290 300 180 171 161 152 142 133

问题四：

用人工手段可以使得该动物群落的数量减少，使其数量控制在一定范围内，但是，我们并不知道用了这种人工手段后，对于一次突发事件后要增加其数量时，该种群的复壮能力如何。因此，我们建立数学模型，分析研究其复壮能力。我们

假设由于某种原因（例如由于一场大灾难导致该动物的死亡），使得该动物减少至M头，此时不得不停止注射避孕的方法。在初步模型中，我们先假设灾难发生后该种群虽然数量减少，但其年龄结构没有改变。建立该动物群落重新壮大的能力的模型如下： 模型七：

?d?(k)?t?60(k)60(k)??xi,1??xi,0－－－－－－－(1)?i?1i?1?dt?60(1)??xi,0?M－－－－－－－－－－－－－(2)??i?1k?1)k)?xi(,0?xi(,1.....(i?0,1,..59)－－－－－－－(3)?k)?xi(,1?0...........(i?60)－－－－－－－－－(4) ?(k)(k)?xi,1??1xi,0......(i?0)－－－－－－－－－(5)?(k)(k)x??x?2i,0.....(i?1,2,...,59)－－－－－－(6)?i,1对M取不同的值，利用c语言进行编程求解，（源程序见附录9）可以得到以下表： M值 第t年 t=1 增量d1 当前数量M1 t=2 增量d2 当前数量M2 t=3 10000 9000 8000 7000 6000 5000 4000 3000 2000 665 594 526 456 7456 489 7945 521 8466 557 388 6388 414 6082 444 7246 472 7718 507 322 5322 342 5664 365 6029 390 6419 418 250 4250 267 4517 286 4803 305 5108 325 179 3179 191 3370 205 3575 219 3794 235 113 2113 122 2235 130 2365 140 2505 145 10665 9594 8526 707 638 562 11372 10232 9088 757 恢复 699 602 增量d3 当前数量M3 10932 9690 729 640 t=4 增量d4 当前数量M4 恢复 10330 9023 689 596 t=5 增量d5

当前数量M5 t=6 11019 9619 741 恢复 638 8225 538 6837 445 7272 472 7744 507 8251 543 5433 348 5781 371 6152 401 6553 426 6979 456 7435 22.8061 4029 251 4280 268 4548 286 4548 307 4835 327 5142 16.4364 2650 156 2806 166 2972 178 3150 191 3341 204 3545 10.7636 增量d6 当前数量M6 10257 8792 682 581 t=7 增量d7 当前数量M7 10929 9373 729 恢复 623 9996 665 t=8 增量d8 当前数量M8 t=9 增量d9 当前数量M9 10661 8794 707 581 t=10 增量d10 当前数量M10 11368 9375 46.6042.6800 57 38.8810 35.6424 28.6121 d?2?t?46.00增长加速度：D? dt200

该表对M进行不同的取值，分别表示某一年（t＝0）灾难发生后该动物群落剩下的数量。我们用第一题的C程序即可算出该群落在下一年（t＝1）的增量，把这个增量加上上一年群落的总数即为第一年（t＝1）群落的总数M1,接着我们又用程序算出这一年的增量，加上M1即为第二年（t＝2）的个体总数M2。以这样的算法推算下去，一直到群落数量恢复到11000左右或者前十年（t＝10）为止。

将计算出来的各年增量按时间t（1，2??10）列在座标上，用MATLAB进行曲线拟合（具体程序见附录10），得到的图形如下所示：

从图中可以看出，每年的增量都是遵从一定的线性关系的，也就是说，该群

d?2?t?落的增长加速度是一定的，通过MATLAB的曲线拟合可以求出不同M

dt2值所对应的增长加速度(程序见附录 )，数据如下： M值 10000 9000 46.6000 8000 42.6857 7000 38.8810 6000 35.6424 5000 28.6121 4000 22.8061 3000 16.4364 2000 10.7636 D（只/46.00年^2） 00

通过这些数据，我们可以知道，在灾难后群落数量减少至M时的重新壮大能力。

再用MATLAB求出M与D的关系，进行曲线拟合求出斜率k后画出图形（程序见附录11），如下所示

可知M与D是服从一定的线性关系的，由MATLAB算出的的数据可知

M＝3.9010+0.0047D。

该模型的稳定性分析见本论文第四部分。

四． 模型的稳定性分析

1．对于问题一：我们利用C语言进行编程，求得了第k年时各年龄结构的分布情况，并利用MATLAB软件对得到的数据进行画图分析，从输出的图可看出各年该动物各个年龄的分布是围绕一定的轨迹进行分布的，从而说明了该模型是稳定的，其次，用统计软件SAS对该模型得出的数据进行相关系数的分析，得到相关系数R都接近于1，而P＝0.0001<0.01则表示其相关性显著。故建立的模型一是很稳定的，其各年的该动物的总数基本上都在11000左右。

2．对于问题二、问题三：我们从两个不同方面对问题进行分析，分别从?＝新生的幼仔数量－自然死亡的数量～0，以及?＝该动物群落在年底时的总数量与年初的数量加上被运走的数量b的差值～0两个不同方面进行模型的建立，从而把该动物群落的稳定性问题转化为求目标函数最小化问题。（这是本论文的一

个亮点）。利用MATLAB对模型进行求解，得出的数据最后都能使得该动物群落维持稳定，可见对于问题二、问题三所建立的模型是很稳定的，结果也是很优的。

3．对于问题四，我们用C语言的编程求得了不同M值在10年内各年的增量，求出该群落在不同M值所对应的增长加速度D。为了求出不同的M与对应的D的关系，我们在MATLAB中输入M与D的数组，进行稳定性分析（程序见附录11）。在回归分析及检验中，我们得到，?0＝3.9010，?1＝0.0047；?0的置信区间为[-1.5746 , 9.3766], ?1的置信区间为[0.0039, 0.0055];

r2=0.9616,F=175.1021,p=0.0000,p<0.05,可知回归模型

M＝3.9010+0.0047D 成立 对模型进行残差分析,作残差图,得图如下:

从残差图可以看出,除最后一个数据外,其余数据的残差离零点均较近,且残差的置信区间均包括零点,这说明回归模型M＝3.9010+0.0047D能较好的符合原始数据,而最后一个数据视为异常点.(造成异常点的原因是M=10000时群落只要两年就可恢复到原来的数量11000,数据太少无法作出较接近实际的曲线拟合。

在预测及作图后得各数据点及回归方程得图形如图：

可以看出,只有最后一个数据点离回归直线距离较远是异常点。 由以上分析可知M与D的关系是比较稳定的。

五． 模型优缺点及改进方向 我们的模型有以下优点：

1．模型的稳定性好，模型给出的方案使得该动物群落的各年龄段的数量保持平衡，从而满足了维持该动物群落的稳定。

2．模型的适用范围广，易于推广，模型对于其它生态经济现象（如最优捕鱼策略问题）同样适用。

3．基本模型对问题的描述准确、合理、推导严谨，理论性强； 4．模型结合实际，具有很高的实用价值。 我们的模型有以下缺点：

1. 由于题目给的数据不够多，所以使得该模型无法更加接近实际的情况。 2. 由于题目没有以前各年偷猎枪杀该动物的数据，故模型没有考虑偷猎枪杀该

动物所产生的影响。

3. 对于一个动物群落，其影响的不确定性因素很多，而模型只是局限在题目中提到的，故与实际情况有一定的偏差。 模型的改进方向：

1．由生活常识可知，对于该动物群落超过60岁后的死亡情况应基本上服从正态分布。故可在模型一中的死亡数中加上这个因素。

2．由于动物的产幼仔是随机分布的，故可利用模拟仿真的方法建立更具有真实性的模型。

参考文献

【1】 樊欣,邵谦谦，SAS 8.X经济统计，北京希望电子出版社，2003.3. 【2】 苏金明，张莲花，MATLAB工具箱应用，北京电子工业出版社，2004.1. 【3】 姜启源，谢金星，数学模型（第三版），高等教育出版社，2003.8. 【4】 余世孝，数学生态学导论，科学科技文献出版社，1995.8.

附录一:（模型一用C语言实现的过程与运行结果） main()

{int x[70],n; /* 定义一个长度为70的数组,每个单位的值表示各年

龄的数量,定义变量n */

int sum,sum1,sum2; /* 三个变量分别代表该年群落个体数量,13—59岁

雌性数量,1—59岁个体数量 */

clrscr(); sum=0; sum1=0; sum2=0;

printf(\

scanf(\/* 输入当年出生的幼仔数量,设幼仔年龄为0岁 */ printf(\

x[1]=0.75*x[0]; /* 幼仔从0岁向1岁过渡时存活率为75% */ for(n=2;n<=60;n++)

{x[n]=0.95*x[n-1]; /* 向下一年龄段过渡时的存活率为95% */ }

for(n=0;n<=60;n++)

{sum=sum+x[n]; /* 对各年龄段数量累加得到种群总体数量 */ printf(\ }

for(n=13;n<=60;n++)

{sum1=sum1+x[n]; /* 对13—60岁个体数量累加 */ }

for(n=1;n<=59;n++)

{sum2=sum2+x[n]; /* 对1--59岁个体数量累加 */ }

printf(\

printf(\/* 13—60岁个体总数的一 printf(\半即为该年龄段雌性数量

再除以3.5即为当年生幼仔的雌性的数量 */

通过不断输入新生幼仔的初值,直到算出当年群落的个体总数与题设数据吻合,则该组数据即为当年群落年龄结构的分布(未算入当年运走的数量).

运行结果:

1.前兩年运走的个体数量为876头,因此前一年的个体总数为11000+876=11876(头),通过调试程序使总数sum达到11876左右. 第一次调试结果:

2．前两年年末即前一年年头所剩下的13—60岁的雌性个体数量决定了前一年年末个体的数量即2470/3.5＝706头； 前两年剩下的数量加上前一年新生幼仔的数量 11000+706=11706，其值为前一年年末的个体数量

用程序算得11714，误差为8头，其误差由于变量是整型变量，数字取整数所致。 第二次调试结果:

在前一年的年末运走622头，即该群落的个体数量剩下11714-622=11092 根据表得到2503/3.5=715，所以11092+715=11807为该群落今年年尾的个体总数。

程序计算得11808，误差为1，其误差由于变量是整型变量，数字取整数所致。

3.对于以上数据是否符合假设,我们还要用一组初步数据來进行对比.因此我们设每年运出的个体数量为0, 个体总数为11000 头,通过调试程序使总数sum达到11000左右.

第三次调试结果:

该结果为当每年没有个体被运走时数量停留在11000左右的数据，用于与初步数据进行对比。

附录2（用MATLAB画出模型一轨迹）： clear all;clc

y1=[792 594 564 535 508 482 457 434 412 391 371 352 334 317 301 285 270 256 243 230 215 207 196 186 176 167 158 150 142 134 127 120 113 107 101 95 90 85 80 75 71 67 63 59 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18]';

y2=[806 604 573 544 516 490 465 441 418 397 377 358 340 322 305 289 274 260 246 233 221 209 198 188 178 169 160 151 143 135 128 121 114 108 102 96 91 86 81 76 72 68 64 60 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18]';

y3=[746 559 531 504 478 454 431 409 388 368 349 331 314 298 283 268 254 241 228 216 205 194 184 174 165 156 148 140 132 125 118 112 106 100 94 89 84 79 75 71 67 63 59 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18 17]';

y4=[800 600 569 540 512 486 461 437 415 394 374 355 337 320 303 287 272 258 245 232 220 208 197 187 177 168 159 151 143 135 128 121 114 108 102 96 91 86 81 76 72 68 64 60 56 53 50 47 44 41 38 36 34 32 30 28 26 24 22 20 18]';

x1=[0 0 0 0 3 4 7 20 9 15 9 22 3 23 5 13 21 0 22 14 5 13 10 0 13 30 14 12 0 20 6 3 5 8 12 10 3 7 14 10 16 21 13 10 12 6 3 6 9 13 10 3 6 21 15 4 13 10 32 14 0]';

x2=[0 0 20 21 13 12 13 22 14 40 14 26 13 14 27 3 14 12 20 25 17 14 10 0 2 3 4 4 3 2 3 13 16 13 10 10 12 16 12 10 12 19 13 24 17 16 25 12 45 23 34 13 16 10 17 13 13 12 3 22 20]'; m1=y1-x1; m2=y2-x2; m3=y3; m4=y4; t=0:60

plot(t,m1,'b.'),hold on,plot(t,m2,'y*'),hold on,plot(t,m3,'g--'),hold on,plot(t,m4,'m-')

legend('前一年剩余数量的年龄结构情况','前两年剩余数量的年龄结构情况','没有被运走时的年龄结构情况','预测当前的年龄结构情况（无运走）') title('今年时该动物群落的年龄结构情况')

附录3（用SAS对模型一进行检验的程序及调试结果）：

data A;

input x1 x2 x3 x4 @@; cards;

594 604 559 600 564 553 531 569 535 523 504 540 505 503 478 512 478 478 454 486 450 452 431 461 414 419 409 437 403 404 388 415 376 357 368 394 362 363 349 374 330 332 331 355 331 327 314 337 294 308 298 320 296 278 283 303

272 286 268 287 249 260 254 272 256 248 241 258 221 226 228 245 216 208 216 232 213 204 205 220 194 195 194 208 186 188 184 197 186 188 174 187 163 176 165 177 137 166 156 168 144 156 148 159 138 147 140 151 142 140 132 143 114 133 125 135 121 125 118 128 117 108 112 121 108 98 106 114 99 95 100 108 89 92 94 102 85 86 89 96 87 79 84 91 78 70 79 86 66 69 75 81 65 66 71 76 55 60 67 72 46 49 63 68 50 51 59 64 49 36 56 60 44 39 53 56 47 37 50 53 47 25 47 50 41 35 44 47 35 -1 41 44 28 18 38 41 28 4 36 38 33 23 34 36 28 18 32 34 11 22 30 32 15 13 28 30 24 15 26 28 13 13 24 26 14 12 22 24 -10 19 20 22

6 -2 18 20 18 -2 17 18 ;

proc corr; var x1 x2 x3 x4; run;输出结果如下：

附录4（用MATLAB求解目标最优化问题） clear all;clc f=[0;0;0;1];

Aeq=[0.2857 0.1818 -1 0;1 1 0 0;0 0 0.75 -1]; beq=[0 2569 529 ]; lb=zeros(4,1);

[x,fval,exitflag,output,lambda]=linprog(f,[],[],Aeq,beq,lb) 输出结果如下：

>> Optimization terminated successfully. x =

1.0e+003 * 2.2934 0.2756 0.7053 0.0000

fval = ％ 返回目标函数值 ％ 2.0286e-016

exitflag = ％ 表示求得的解是收敛的 ％ 1 output =

iterations: 6 cgiterations: 0

algorithm: 'lipsol' ％ 求目标最优化的方法 ％ lambda = ％ 所占的空间 ％ ineqlin: [0x1 double] eqlin: [3x1 double] upper: [4x1 double] lower: [4x1 double]

附录5（用MATLAB求解模型三目标最优化问题） clear all;clc f=[0;0;0;1];

Aeq=[0.2896 0.1843 -1 0;1 1 0 0;0 0 0.75 -1]; beq=[0 2569 529 ]; lb=zeros(4,1);

[x,fval,exitflag,output,lambda]=linprog(f,[],[],Aeq,beq,lb) 输出结果如下：

> Optimization terminated successfully. x =

1.0e+003 * 2.2020 0.3670 0.7053 0.0000 fval = 2.4791e-016 exitflag = 1 output =

iterations: 6 cgiterations: 0 algorithm: 'lipsol'

附录6: （用MATLAB求解模型四目标最优化问题） clear all;clc f=[0;0;0;0;1]; A=[0 0 1 0 0]; b=[242];

Aeq=[0.2896 0.1843 0.1698 -1 0;1 1 1 0 0;0 0 0 0.75 -1]; beq=[0 2569 529]; lb=zeros(5,1);

[x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb) 输出结果如下：

>> Optimization terminated successfully. x =

1.0e+003 * 2.2171

0.2421 ％由此可知当?＝242时模型稳定％ 0.1098 0.7053 0.0000 fval = 1.6409e-012 exitflag = 1 output =

iterations: 6 cgiterations: 0 algorithm: 'lipsol' lambda =

ineqlin: 1.9145e-014 eqlin: [3x1 double] upper: [5x1 double] lower: [5x1 double]

附录7（用MATLAB求解模型五目标最优化问题） clear all;clc f=[0;0;0;1];

Aeq=[0.2896 0.1843 -1 0;1 1 0 0;0 0 0.75 -1]; beq=[0 2569 579]; lb=zeros(4,1);

[x,fval,exitflag,output,lambda]=linprog(f,[],[],Aeq,beq,lb) 输出结果如下：

>> Exiting: One or more of the residuals, duality gap, or total relative error

has grown 100000 times greater than its minimum value so far: the primal appears to be infeasible (and the dual unbounded). (The dual residual < TolFun=1.00e-008.) x =

1.0e+003 * 2.5685 0.0005 0.7439 0.0000 fval = 4.8164e-006

exitflag = -1 output =

iterations: 3 cgiterations: 0 algorithm: 'lipsol' lambda =

ineqlin: [0x1 double] eqlin: [3x1 double] upper: [4x1 double] lower: [4x1 double]

附录8（用MATLAB求模型六目标最优化问题）

clear all;clc for i=50:10:300

％表明得到的数据是发散的％

f=[]; A=[]; b=[]; Aeq=[]; beq=[]; lb=zeros(1,1);

[x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb) end

输出结果如下： x ＝

370 361 351 342 332 323 313 304 294 285 275 266 256 247 228 218 209 199 190 180 171 161 152 142 133

附录9（用C语言求不同M值时各年的增量） main() {int x[70],n,d; int sum,sum1,sum2; clrscr(); sum=0; sum1=0; sum2=0;

printf(\ scanf(\

printf(\ x[1]=0.75*x[0]; for(n=2;n<=60;n++) {x[n]=0.95*x[n-1]; }

for(n=0;n<=60;n++)

{sum=sum+x[n];

printf(\ \ /*打印出各年龄段的数量*/ }

for(n=13;n<=60;n++) {sum1=sum1+x[n]; }

printf(\ /*用总和来控制，使之等于M在第t年

数量*/

printf(\/*每年的增量约等于13—16岁雌性个体

的数量*/

}

附录10（求不同M值增量d的关系） clear all;clc

t=[1 2 3 4 5 6 7 8 9 10]; t1=[1 2]; t2=[1 2 3]; t5=[1 2 3 4] t3=[1 2 3 4 5 6]; t4=[1 2 3 4 5 6 7 8]; M1=[665 707 757]; M2=[594 638 699 729]; M3=[526 562 602 640 689 741];

M4=[456 489 521 557 596 638 682 729]; M5=[388 414 444 472 507 538 581 623 665 707]; M6=[322 342 365 390 418 445 472 507 543 581]; M7=[250 267 286 305 325 348 371 401 426 456]; M8=[179 191 205 219 235 251 268 286 307 327]; M9=[113 124 134 141 152 161 173 185 199 212];

f=inline('k(1)*t+k(2)','k','t'); k1=[0.5 30];

m1=lsqcurvefit(f,k1,t2,M1) z1=f(m1,t2) k2=[0.5 30];

m2=lsqcurvefit(f,k2,t5,M2) z2=f(m2,t5) k3=[0.5 30];

m3=lsqcurvefit(f,k3,t3,M3) z3=f(m3,t3) k4=[0.5 30];

m4=lsqcurvefit(f,k4,t4,M4) z4=f(m4,t4) k5=[0.5 30];

m5=lsqcurvefit(f,k5,t,M5) z5=f(m5,t) k6=[0.5 30];

m6=lsqcurvefit(f,k6,t,M6) z6=f(m6,t) k7=[0.5 30];

m7=lsqcurvefit(f,k7,t,M7) z7=f(m7,t) k8=[0.5 30];

m8=lsqcurvefit(f,k8,t,M8) z8=f(m8,t) k9=[0.5 30];

m9=lsqcurvefit(f,k9,t,M9) z9=f(m9,t)

plot(t2,M1,'b.'),hold

on,plot(t2,z1,'r-'),plot(t5,M2,'b.'),hold

on,plot(t5,z2,'r--'),plot(t3,M3,'b.'),hold on,plot(t3,z3,'y-'),plot(t4,M4,'b.'),hold on,plot(t4,z4,'y--'),plot(t,M5,'b.'),hold on,plot(t,M6,'b.'),hold

on,plot(t,z5,'g-'),hold

on,plot(t,z6,'g--'),plot(t,M7,'b.'),hold

on,plot(t,z8,'m--'),hold

on,plot(t,z7,'m-'),plot(t,M8,'b.'),hold on,plot(t,M9,'b.'),hold on,plot(t,z9,'k-') title('M取不同的值时每年增量的分布情况')

其运行结果如下：

Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m1 =

46.0000 617.6667 z1 =

663.6667 709.6667 755.6667 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m2 =

46.6000 548.5000 z2 =

595.1000 641.7000 688.3000 734.9000 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m3 =

42.6857 477.2667 z3 =

519.9524 562.6381 605.3238 648.0095 690.6952 733.3810 Optimization terminated successfully:

Norm of the current step is less than OPTIONS.TolX

m4 =

38.8810 408.5357 z4 =

Columns 1 through 7

447.4167 486.2976 525.1786 564.0595 602.9405 641.8214 680.7024 Column 8 719.5833

Optimization terminated successfully:

Norm of the current step is less than OPTIONS.TolX m5 =

35.6424 337.8667 z5 =

Columns 1 through 7

373.5091 409.1515 444.7939 480.4364 516.0788 551.7212 587.3636

Columns 8 through 10 623.0061 658.6485 694.2909 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m6 =

28.6121 281.1333 z6 =

Columns 1 through 7

309.7455 338.3576 366.9697 395.5818 424.1939 452.8061 481.4182

Columns 8 through 10 510.0303 538.6424 567.2545 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m7 =

22.8061 218.0667 z7 =

Columns 1 through 7

240.8727 263.6788 286.4848 309.2909 332.0970 354.9030 377.7091

Columns 8 through 10 400.5152 423.3212 446.1273 Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m8 =

16.4364 156.4000 z8 =

Columns 1 through 7

172.8364 189.2727 205.7091 222.1455 238.5818 255.0182 271.4545

Columns 8 through 10 287.8909 304.3273 320.7636 Optimization terminated successfully:

First-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected m9 =

10.7636 100.2000 z9 =

Columns 1 through 7

110.9636 121.7273 132.4909 143.2545 154.0182 164.7818 175.5455

Columns 8 through 10 186.3091 197.0727 207.8364

附录11（M、D关系的稳定性检验）

x=[2000 3000 4000 5000 6000 7000 8000 9000 10000];

y=[10.7636 16.4364 22.8061 28.6121 35.6424 37.5357 40.4000 46.6000 46.0000];

f=inline('k(1)*x+k(2)','k','x'); k=[0.5 30];

m=lsqcurvefit(f,k,x,y) z=f(m,x)

plot(x,y,'.'),hold on,plot(x,z,'-') 其输出结果如下：

>> Optimization terminated successfully:

Relative function value changing by less than OPTIONS.TolFun m =

0.0046 4.0893 z =

Columns 1 through 7

13.2742 17.8666 22.4591 27.0516 31.6440 36.2365 40.8290

Columns 8 through 9 45.4214 50.0139

附录12（问题四稳定性分析）

>> x=[2000 3000 4000 5000 6000 7000 8000 9000 10000]'; >> X=[ones(9,1)x];

Warning: Future versions of MATLAB will require that whitespace, a comma, or a semicolon separate elements of a matrix. Please type

\matrix_element_separators\at the MATLAB prompt for more information.

>> Y=[10.7636 16.4364 22.8061 28.6121 35.6424 38.8810 42.6857 46.6000 46.0000]';

>> [b,bint,r,rint,stats]=regress(Y,X); >> b,bint,stats b = 3.9010 0.0047 bint =

-1.5746 9.3766 0.0039 0.0055 stats =

0.9616 175.1021 0.0000 >> rcoplot(r,rint) >> z=b(1)+b(2)*x z = 13.2832 17.9743 22.6653 27.3564 32.0475 36.7386 41.4296 46.1207 50.8118

>> plot(x,Y,'k+',x,z,'r')

B题：人员安排问题

问题综述：该问题主要是为了求解在客户的要求下公司每天收益的最大化，属于优化问题；

我们在对这个问题建模时，主要是基于客户的两个要求来建立的： （1）客户对员工的人数要求； （这个要求是本来题目有的） （2）客户对工期的要求； （这个要求是我们进一步假设的）

对于第一个要求我们建立了基本模型，而对于第二个要求，我们在第一个要求的基础上，进一步改进了基本模型，从而建立了某个项目先完工的模型。具体的解题思路如下图所示：

Max：公司一天的直接收益

只考虑客户对员工的人数要求 进一步考虑客户对工期的要求

改进 基本模型 某个项目先完工的模型

一．模型基本假设：

1．假设客户对项目的工期没有限制，项目的工期由公司决定，且四个项目同时开工，同时完工，中间也不停工。 2． 假设所有人员总能在岗位上工作，不考虑由于生病或是其他意外事件而造成人员的缺席。 3．假设四个项目同时需要的最多人数不超过现有公司工作人员的人数，即使超过，也只分配公司现有的工作人员。

4．假设C、D两个项目的管理费由公司支付；

5．假设所有工作人员都安排完毕，即每个人都有工作。

6．假设同等级别的工作人员的技术水平是相同的，即他们可以接受任意等同的任务。

二．符号说明：

i：用i＝1,2,3,4分别表示高级工程师，工程师，助理工程师和技术员。 j：用j＝1,2,3,4分别表示项目A,B,C和D。

Xij：公司分配第i级别工作人员到第j个项目上的人数。例如X23表示公司分配工程师到

项目C上的人数。

aij：第i级别工作人员分配到第j个项目上的收费。

。 bij: 第i级别工作人员分配到第j个项目上时公司的开支（包括工资和管理费）

Aij: 表示到项目j工作的第i级别工作人员为公司贡献的纯利润收入。

。 ?j： 表示第j个项目的总工时（即项目j的总工作量）

。 Tj： 表示第j个项目客户所要求的工期（即项目j所需要的完工时间）

Mj： 表示客户要求第j个项目一天所必须完成的工作量。

mj： 表示公司分配给第j个项目的所有工作人员一天能够完成的工作量。 三．模型的建立：

1．基本模型：

公司每天的直接收益等于公司每天对各个项目的总收费减去公司每天的总支出，公司每天的总收费为：

??a44ij Xij，每天的总开支为：??bijXij，于是公司每天的直接收益为：

44i?1j?1i?1j?1444444f(Xij)???aijXij－??bijXij=i?1j?1??(aij?bij)Xij i?1j?1i?1j?1在客户的人数要求下，可得基本模型为： 44Max：f(Xij)???(aij?bij)Xij

i?1j?14s.t.

?Xi1?10； i?1?4Xi2?16； i?1?4Xi3?11； i?1?4Xi4?18； i?1?4X1j?9； j?1?4X2j?17； j?1?4X3j?10； j?1?4X4j?5； j?11?X11?3; 2?X21?17;

1）

（2）

3）

4）

5）

6）

7）

8）

（9）

（10） （11）

（（（（（（（

2?X31?10; （12） 1?X41?5; （13） 2?X12?5; （14） 2?X22?17; （15）

2?X32?10; （16） 3?X42?5; （17）

X13?2; （18）

2?X23?17; （19） 2?X33?10; （20） 1?X43?5; （21） 1?X14?2; （22） 2?X24?8; （23）

1?X34?10; （24） X44?0; （25）

（i,j=1,2,3,4） (26) Xij为整数；

此模型即为在这些约束条件下对公司一天的收益求解最大值，属于整数线性规划模型。

2．基本模型的求解：

为了便于数据的对照和模型的求解，我们先把题目所提供的表2的结构改成与表3相同的结构，如表1所示： 表1

收费（元/天） 高级工程师 工程师 助理工程师 技术员 A 1000 800 600 500 B 1500 800 700 600 C 1300 900 700 400 D 1000 800 700 500 利用Matlab和Lingo软件求解基本模型得出解如表2所示：（具体解法与源程序请见附录1和附录2，至于我们为什么两个软件都用，请务必看附录2中的注3的解释，谢谢！）

表2 分配（人） 高级工程师 工程师 助理工程师 技术员 A 1 6 2 1 B 5 3 5 3 C 2 6 2 1 D 1 2 1 0 并解得公司一天直接收益最大值为：27150(元) 3．结果分析：

在符合模型基本假设的条件下，该模型的解既满足了客户的人员要求，又实现了公司每天收益的最大化。另外，从解中可以知道，公司的人员分配满足了项目A,B,C的最大人员要求，然而项目D的人员却只有4人，显然这个解是符合公司利益的，因为从收费表中知道，项目D的收费相对较低，而且公司还要付出一定的管理费，从公司的收益来考虑，当然希望把最少的人分配到项目D上，而把更多的人分配到其他项目上去，但是这种情况也恰恰从另一方面反映了这个模型的简单与不足，那就是这个模型更多的是基于公司的利益来考虑问题的，虽然满足了客户对人员的要求，但并没有考虑到现实中客户对项目的工期等其他方面的要求，所以需要进一步改进。

四：模型的改进：

1．某个项目j先完工的模型：

在基本模型中，我们主要是基于模型基本假设1来建立模型的，为了更加符合现实情况，我们将基本假设1修改为：

（1）假设各项目同时开工，但客户要求某项目j先完工，而其他项目可以后完工，而且同时完工。

并进一步假设如下：

（2）假设不同级别的工作人员的工作能力与其工资水平成正比。 （3）假设客户对项目的工期要求是合理的。 在这三个假设下，我们令一个技术员一天能够完成的工作量为一个标准工时，由于不同等级的工作人员的工作能力（即一天能完成的工作量）与其工资水平成正比，所以一个高级工程师，工程师和助理工程师一天能够完成的工作量分别相当于：250/110,200/110,170/110 个工时，假设项目j完工总共需要的工时为?j，而客户要求项目必须在Tj天内完成，则项目

j 每天的工时为：

Mj??jTj （27）

而该项目所分配人员一天所能完成的工作量记为：

mj?250200170X1j?X2j?X3j?X4j （28） 110110110对于客户给定的?j和Tj，我们可以确定Mj，于是得到模型为：

Max：f(Xij)???(ai?1j?144ij?bij)Xij

s.t. （2）～（26）式，

250200170X1j?X2j?X3j?X4j?Mj (29) 110110110因此某项目j先完工的模型就是在基本模型上加多一条约束式（29）。 mj?2．Mj的确定：

由于公司所分配给每个项目的工作人员是有限制的，所以他们在一天内所能完成的工作量也是有限制的，这样，在项目j的总工作量?j一定的情况下，如果客户要求的工期Tj过短，那么每天所要求完成的工作量Mj有可能超过公司分配给该项目的工作人员一天能够完成的工作量mj，这就涉及到客户要求的工期的合理性的问题，在总的工作量?j一定时，如果客户所要求的工期Tj使得每天的工作量Mj在mj的范围[mjmin mjmax]内，那么就是合理的，如果超出这个范围那么就是不合理的。为此，我们先求出在客户所要求的人数限制条件下，公司所分配给该项目的工作人员一天所能完成的工作量mj的范围[mjmin

mjmax]，具体如下：

注1：mj的范围[mjmin mjmax]表示mjmin?mj?mjmax，以下同。

我们所要求的就是在约束条件（2）～（26）下mj的范围，例如对于项目D（m4），现在求

m4的范围，具体求解时可依照基本模型的解法，只需将目标函数

Max：f(Xij)?Max：m4???(ai?1j?144ij?bij)Xij 改为

250200170X14?X24?X34?X44 和 110110110250200170X14?X24?X34?X44 Min：m4?110110110即可以求出m4的范围在最小值m4min和最大值m4max之间，解得：

7.45?m4?25.27 （30）

依照该方法可求出其他三个项目所分配人员一天所完成的工作量的范围，将得到的解表示在

表3中：

表3

项目（mj） 范围 A（m1） [10.00 18.18] B（m2） [14.27 28.36] C（m3） [12.27 19.54] D（m4） [7.45 25.27] 从表3知道，要确定Mj，只需要给定的值Mj在这些范围内即可。

3．模型的求解：

在这个模型中，我们具体求解的是项目D先完成的模型，并在满足(30)式的范围内，给定项目D客户每天要求的工作量M4。

但由于如果单取M4为一个值时，该模型可能没有整数解，所以这里具体取M4在满足（30）式下的更小的范围12?M4?14,那么依照基本模型的解法，可解得新的人员安排如表4所示： 表4 分配（人） 高级工程师 工程师 助理工程师 技术员 总计（各项目） A 1 3 2 1 7 B 5 3 5 3 16 C 2 6 2 1 11 D 1 5 1 0 7 总计（各级别人员） 9 17 10 5 41 并解得公司一天直接收益最大值为：27000(元) 而取22?M4?24时的解如表5： 表5 分配（人） 高级工程师 工程师 助理工程师 技术员 总计（各项目） A 1 2 2 1 6 B 5 2 2 3 12 C 2 5 2 1 10 D 1 8 4 0 13 总计（各级别人员） 9 17 10 5 41

并解得公司一天直接收益最大值为：26650（元)

取M4在满足（30）式下的不同范围内，将得到公司每天收益不同的最优解，具体如表6所示：（这里只列出收益最优解，不列出新的人员安排解） 表6 M4的范围 公司每天收益最优解 [12 14] [14 16] 27000 26950 [16 18] 26900 [18 20] 26850 [20 22] 26750 [22 24] 26650

在实际中M4取值的变化则表现为如果客户对工期的要求不同，将导致公司对各个项目人员的分配发生改变，来维持公司每天收益的最大化。 4．项目D完成后公司一天的收益求解：

项目D完成后，根据基本假设6，公司将重新对A,B,C三个项目进行人员分配,当然，在分配的时候，原先在各自的项目工作的人员，只要该项目人员不减少，那么在实际中仍然保留这些人员。此时公司每天的直接收益可依照基本模型的方法求解，具体只要将对D 项目的约束式（1）和（22）～（25）解除，并注意此时公司的人员并不一定全部能够分配到位，即此时基本假设5不作为该改进模型的假设，因此将约束式(6),(7),(8),(9)做如下修改：

?Xj?1441j?9；

?Xj?1441j?9 （6）

?Xj?142j?17； ?X2j?17 （7）

j?1?Xj?143j?10； ?X3j?10 （8）

j?144?Xj?14j?5； ?X4j?5 （9）

j?1即将约束式中的等号改为小于或等于号，求得其解如表7所示： 表7 分配（人） 高级工程师 工程师 助理工程师 技术员 总计（各项目） A 2 5 2 1 10 B 5 6 2 3 16 C 2 6 2 1 11 D 0 0 0 0 0 总计（各级别人员） 9 17 6(未满人) 5 37（未满人）

并解得此时公司一天直接收益最大值为：24550(元)

且由这个表知道公司的员工并没有全部分配完，只有10+16+11=37(人)。 5．结果分析：

在这个改进模型中，由于假设了D项目最先完成，而其他项目没有限制时间，所以我们可以看到结果公司分配给D项目的人员增多了，由原来的1+2+1=4(人)增加到1+5+1=7(人)再增加到1+8+4=13（人），并且，更具体的分析，如果M4的取值越大，也就是每天的工作量越大，那么客户所要求的工期就是越短的，所以分配的员工就越多，我们可以看到当M4的取值由范围[12 14]变到[22 24]时，员工由原来的1+5+1＝7（人）增加到1+8+4=13（人）。这个结果是符合实际情况的。另一方面，我们可以从表6看到随着M4取值的增大，公司的收益是逐渐降低的，这是由于：

（1）分配到项目D 的人数增加了，结果导致管理费的增加；

本文来源：https://www.bwwdw.com/article/5q65.html