第七章 系统模拟 第一节 系统模拟概述 第二节 系统 模拟示例 第三节 蒙 特卡罗模拟 第四节 系统 动力学

系统模拟(Simulation)又称系统仿真。 它作为研究、分析和设计系统的一种有效技术 正被广泛应用。农业系统的研究也不例外。模拟技 术的不断实践也扩展了自身的概念。可以说凡是利 用计算机在计算机上模拟而不是在真实系统上进行 实验、运行的研究方法都可认为是模拟。

第一节 系统模拟概述 一、概念 模拟（Simulation） 是对真实过程或系统在整个时间内运行的模仿。 两种情况: （1）在某些情况下，所研究的模型足够简单，可以用数 学方法表示并能求解，从而获得所关心的系统性能参数； （2）许多真实系统是非常复杂的，无法用数学关系或数 学方法来求解。这时利用仿真就可以像观察测试真实系统那 样，在模拟模型上得到系统性能随时间而变化的情况，得到 有关系统的性能参数，从模拟中获得想要的东西。

系统模拟技术主要包括三个过程： （1）建立模型 是实现模拟的基础。即先对被模拟的系统进行分析、假 设、简化，从而得出能反映真实系统本质及变化的模型，再 将模型转变成计算机能运行的模拟模型。 （2）对模型进行实验、运行 （3）对模拟的有效性分析 如：模型有效性分析、输入输出数据分析，输出结果与 现有真实数据的比较验证等。

模拟一般可分为两大类： （1）物理模拟或实物模拟 构造一个物理模型进行实验。 （2）计算机模拟 通过在计算机上模拟系统的变化，通过模拟输 出的数据得到各种结果。 计算机模拟又可分为离散系统模拟和连续系统 模拟两大类。

二、离散事件系统模拟 离散系统 是所研究系统的状态变量在一些离散时间 点上发生变化的系统。 特定时刻 这些离散的时间点。 在这些特定时刻系统状态发生变化，而在相邻的两个特 定时间之间的其它任何时刻系统状态保持不变。而在这些特 定时刻是由于有事件发生从而引起了系统状态变化。 常见的离散系统有排队系统、库存系统等。 主要特征 随机性——一个或多个输入为随机变量，而不是确定变 量。所以它的输出也往往是随机变量。

蒙特卡罗统计模拟的思想 用模拟随机变量的方法(通常用计算机获取一个随机变 量的一个现实，即所谓抽签，来不断获得随机输入，用所谓 系统模型来产生系统的输出。）实现对系统的模拟。 在这类系统的模拟中对随机型输入、输出进行分析是模 拟有效性分析的一个重要内容。

（1）微分方程的数值解法 三、连续系统模拟 连续系统 是指系统中的状态变量随时间连续变化的系统。 由于连续系统的模型要描述系统实体变量的速率R，所 以连续系统模型通常都是由微分方程组成的。当系统比较复 杂尤其是引进非线性因素后，此微分方程经常不可求解，至 少非常困难。所以通常也就采用模拟方法进而获得答案或者 近似结果。 系统两类： （1）微分方程的数值解法

计算数学的一个分支。要是算法 。 常用的方法： 欧拉（Euler）法，龙格—库塔（Runge-kutta） 法及阿达姆斯（Adams）法等。 对连续系统进行模拟，首先要保证这一数值解的稳定性, 即在初始值有误差，计算机在舍入误差影响下，误差不会积 累而导致计算失败，确保其解在满足一定精度下的正确性。 （2）离散相似法 是将连续系统进行离散化处理，包括按照时间以及其它 变量离散化。用离散化的模型直接代替连续系统的模型。一 般是将微分方程用其相应的差分方程来代替或等效。

比如:常系数微分方程用其常系数差分方程来等效。通过 选取适当的离散间隔就可以方便地用迭代方法在数字计算机 上直接求解差分方程。只要采取了合理的算法，就可以将连 续模型变成模拟模型，然后编程、运行得到系统的模拟结果。 系统动力学（System Dynamics） 由美国麻省理工学院J.W福雷斯特(Jay. w. Forrester)教授 所创。 以差分方程代替微分方程为重要思想基础，并借助系统 各要素之间的因果关系，以及将系统作为信息反馈系统，从 而建立连续系统模拟模型的一种方法 。 系统动力学软件 即所谓DYNAMO计算机语言及软件系统。

四、系统模拟的一般步骤 首先要确定问题，并明确系统模拟的目的和任 务。确定问题包括确定系统的边界和组成部分以及 评价模拟结果的指标等。 1、拟定问题 首先要确定问题，并明确系统模拟的目的和任 务。确定问题包括确定系统的边界和组成部分以及 评价模拟结果的指标等。 2、收集和整理数据 在系统模拟中需要大量的数据。数据的可靠性对模拟输 出结果的正确性影响很大。 3、建立模拟模型 把所研究的系统进行适当的抽象化，从而建立数学模型 或结构模型。 4、转换模型 应用通用的计算机语言或专用的模拟程序语言，将系统 模型或结构模型转变为由计算机程序组成的模拟模型，以便 在计算机上进行模拟运行。

5、验证模型 进行调试运行、模拟并分析、验证其模型是否能正确反 映现实系统的本质以及模拟模型是否正确实现原型。在调试 中应及时修改补充模型（包括模拟模型、模拟程序）甚至调 整补充数据，直到满意为止。 6、模拟运行方案的制定 根据目的及系统的真实条件（包括约束条件，环境条件） 建立模型运行的试验条件，包括设定系统运行的初始条件及 种类，模拟运行时间或次数等。 7、模拟运行 对所研究的系统进行大量的模拟运行，获取大量的输出 资料。

五、模拟语言 8、模拟结果分析 对模拟结果进行分析，包括对结果的各种统计特性分析, 与已知部分有关信息及预期目标等的比较分析，进一步说明 其模拟结果的可靠性与有效性，并做出模拟结论。 五、模拟语言 编程可采用的两类语言： （1）仿真专用语言 离散系统模拟语言— GPSS语言； 连续系统模拟语言—面向微分方程的CSS1、面向结构图 的CSS2、 ACSL 、DYNAMO等。 （2）通用语言 FORTRON、BASIC、C++等。

六、模拟的应用 1、模拟技术在农业中的应用 如：农产品销售库存系统动态模拟，作物生长 过程模拟,以及农业生产过程、病虫害群体的发生、 发展和防治以及水利灌溉模拟等。 农业系统是一个随机因素更多，不可控因素更多的复杂 系统。应用模拟技术对提高农业效益，减少风险，实现可持 续发展有着重要意义。 2、其它方面的应用 排队系统的模拟已被应用于各类服务系统中，如加油站 服务系统、交通通道系统、通信交换服务系统等。 计算机集成制造系统、专门操作人员的培训 。

第二节 系统模拟示例 例1： 一个大型超市每日都从农村采购新鲜农产品 出售，正常情况下每公斤可获例1元。如果采购数量 例1： 一个大型超市每日都从农村采购新鲜农产品 出售，正常情况下每公斤可获例1元。如果采购数量 过多，次日只能减价出售，每公斤将亏损0.4元，现 在该市采用以下采购策略：以前一天的市场需求量 作为当天的采购量。据统计分析，每天平均需求量 为100 kg,标准差为30 kg。在这种情况下，该超市经 营一个月能获多少利润？

第一步 建立数学模型 设Qt为第t日的采购量(kg); Dt 为第t日的需求量（kg); P 为正常出售每千克利润，P = 0.1元/kg；L为减价出售每千克 赔本额，L=0.04元/kg；Rt为第t日的总利润（元）；Tt为第t 日的累积利润（元）根据变量之间的关系，可得该问题的 数学模型为：

第二步 作出模拟框图, 将数学模型转变为计算 机模拟模型 开 始 t=0, D0=100, T0=0 t=t+1 Qt=Dt-1 产生Dt 第二步 作出模拟框图, 将数学模型转变为计算 机模拟模型 产生Dt t=0, D0=100, T0=0 t=t+1 Qt=Dt-1 开 始 Dt≤Qt Rt=0.1Dt-0.04(Qt-Dt) Tt=Tt-1+Rt T<30 停 止 Rt=0.1Qt 否 是

由题知,我们可以取需求量Dt服从于正态分布N(100,302)。 产生30个正态分布随机数，于是可将它们作为随机需求量进 第三步 上机运行，得出模拟结果 由题知,我们可以取需求量Dt服从于正态分布N(100,302)。 产生30个正态分布随机数，于是可将它们作为随机需求量进 行模拟计算。现选定t0=0，D0=100,T0=0,得到超市30天总利润。 T Dt Qt Rt Tt 备注 100 1 104.02 10 Qt<Dt 2 118.37 10.4 20.4 3 68.6 4.87 25.27 Qt>Dt 4 124.69 6.86 32.13 5 69.44 4.73 36.86 6 99.64 6.94 43.81 7 111.69 9.96 53.77 8 125.91 11.17 64.94 9 46.66 1.5 66.44 148.31 4.67 71.1 11 67.32 3.49 74.6

12 112.59 67.32 6.73 81.33 Qt<Dt 13 127.36 11.26 92.59 14 74.81 5.38 97.97 Qt>Dt 15 98.68 7.48 105.45 16 45.14 2.37 107.82 17 121.82 4.51 112.33 18 74.29 5.53 117.86 19 77.41 7.43 125.29 20 85.76 7.74 133.03 21 46.07 8.58 141.61 22 105.14 8.88 150.49 23 112.49 10.51 161 24 80.62 6.79 167.79 25 82.03 8.06 175.85 26 135.65 8.2 184.05 27 87.25 190.84 28 73.59 6.81 197.65 29 62.85 5.85 203.51 30 118.39 6.28 209.79

例2: 这是一个关于空调机生产发展问题。某区随着家庭户 数的增加，要求成套住房数也随之增加。房数的增加，则要 例2: 这是一个关于空调机生产发展问题。某区随着家庭户 数的增加，要求成套住房数也随之增加。房数的增加，则要 求供应空调机数量增加。随着已售出空调机数量增加，市场 就越来越饱和，售出空调机数量反而下降。对住房也是如此。 而家庭户数增加，人口增加，成年人增加，因而结婚数增 加，使得家庭户数又增加。某管理机构，想了解家庭户数， 售出的房数，空调机售出数量的变化趋势。为了分析方便， 现令R1为家庭户数增长率；R2为售出住房增长率；R3为空 调机出售增长率；H为家庭户数；Y为已售出住房数；X为空 调机已售出数量。并假定： R1=K1H 式中： H-Y—潜在的住房市场需求量 R2=K2(H-Y) Y-X—为潜在的空调机市场需求量 R3=K3(Y-X)-K4X K1,K2,K3—常数 K4—损坏率，也是一个常数

空调机发展关系图 家庭户数增 长率（R1) 售出住房增 长率（R2) 空调机出售增 长率（R3) 家庭户数（H) 已售出住房 数（Y) 已售出空调 机数量（X) H-Y Y-X k1 k2 k3 k4 空调机发展关系图

第一步 建立数学模型 第二步 将数学模型转化为模拟模型 具有递推关系的差分方程组：

只要给定初始条件，确定出有关常数，则可编程得到模 拟模型。现假定所获的有关数据为： （1）100个家庭平均10个月增加一个家庭； （2）一幢住房从开始建筑到销售使用要经过5个月的时间; （3）空调机从订货到安装要经历10个月时间； （4）售出的空调机平均25个月损坏一台（损坏后不能修 复）。 于是可知：

k1 = 1/（100×10）= 0.001户/(月·户)； k2 =1/T2 =1/5=0.2 住房数/月； k3=1/T3=1/10=0.1空调机数/月； k4=1/T4=1/25=0.04台数/月（损坏）。 设初始条件为： H0=1000 户；Y0=0；X0=0；模拟时间步长△t = 1月。于是 有差分方程为（代入常数后）： Ht=Ht-1+0.001Ht-1 Yt=Yt-1+0.2(Ht-Yt-1) Xt=Xt-1+0.1(Yt-Xt-1)-0.04Xt-1

空调机趋势模拟计算表 模拟次数 月 Ht Yt Xt 1000 1 1001 200.2 20.02 2 1002.001 360.56 1000 1 1001 200.2 20.02 2 1002.001 360.56 53.273 3 1003.002 489.086 98.98 4 1004.002 592.039 152.25 …

第三节 蒙特卡罗模拟 蒙特卡罗（Monte Carlo）模拟，又称模拟抽样法或统计 试验法，是一种以概率与数理统计理论为指导的模拟技术。 它的实质是用按一定的概率分布产生随机数的方法来模拟可 能出现的随机现象。 蒙特卡罗法通过抓住事物运动过程的数量和物理特征,运 用数学方法来进行模拟。每一次模拟描述系统可能出现的一 次情况，经过成百上千次模拟后，便可得到一些很有价值的 结果。 如：顾客需求行为分析、交叉影响分析、经济计量模型模拟 求解，预测方法的评价和误差分析等等，都可以运用蒙特卡 罗模拟技术来求解。

（1）在某些情况下，用一个适当的理论分布（诸如正 蒙特卡罗模拟法求解问题的三个主要步骤： 1、构造模型，确定描述研究对象行为的概率分布。 （1）在某些情况下，用一个适当的理论分布（诸如正 态分布或泊松分布）来描述研究对象的经验概率分布，既是 可能的，也是可取的（例如，在一定时间内，服务台到达的 顾客量即服从泊松分布，随机因素所产生的干扰常服从正态 分布）。理论分布往往更能代表所研究的那个“总体”。用 于确定经验分布的抽样信息就是从这个“总体”中得出来的。 （2）某些经济问题常常没有可以直接引用的分布律。通 常的作法是根据历史记录或主观的分析判断求得研究对象的 一个初始概率分布。例如：在需求预测中，可以根据过去的 实际需求量分布状况估计预测目标的初始分布，或运用主观 概率法、专家调查法给出一个事件出现的概率分布。

2、运行模拟 即根据确定的模型结构（概率分布及其结构 关系）进行随机抽样，故又称为数字模拟。 具体作法： 将研究目标的概率分布映射到一实数区间(常取0 ~1或1 ~ 100), 然后利用随机数表或随机数发生器产生随机数，若该随 机数落在某一事件发生概率所对应的区间内，则认为该事件 发生了。 例如:某一研究目标出现E1、E2、E3的概率P(E1)、P(E2)、P(E3) 分别为0.2、0.5、0.3，将其分别映射到0~19，20~69，70~99区 间，然后任意地从随机表中抽取一数，例如抽取的数为46，此 数在20~69之间，故我们就认为E2事件发生了一次。

随 机 模 拟 示 意 图

发生的频数，并运用数理统计知识求取各种统计量。 实际模拟中，往往要规定一组不同的作业条件（特性), 完 3、根据模型约定随机模拟结果和预测要求，统计各事件 发生的频数，并运用数理统计知识求取各种统计量。 实际模拟中，往往要规定一组不同的作业条件（特性), 完 成一种条件下的模拟后，在返回到第二步骤进行另一种条件 下的模拟。如此经过多次模拟试验，再对不同条件下的模拟 结果进行比较，选择一个最接近实际系统的作业条件和一个 最佳模拟方案。

蒙特卡罗模拟工作流程图 确定初始概率分布 数字模拟（随机抽样） （多方案分析） 修正初始概率分布 计算模拟统计量 不合理 进行评价与预测分析 输出模拟结果 （多方案分析） 修正初始概率分布 合理 不合理 蒙特卡罗模拟工作流程图

例3 : 某企业计划投资建立某产品生产线，为此必须对该产 品未来能实现多少利润进行预测和分析。建立该生产线需投 例3 : 某企业计划投资建立某产品生产线，为此必须对该产 品未来能实现多少利润进行预测和分析。建立该生产线需投 资5万元。产品能实现多少利润主要受以下三个不确定因素 的影响：售价、成本与年销售量。经过有关生产、计划、销 售人员分析，考虑到原材料供应、市场竞争和价格浮动等因 素的作用，初步估计售价、成本与年销售量可能出现的情况 及其发生概率如下表。 售价 （元） 发生概率 成本 年销售量 （万件） 5 0.3 2 0.1 3.5 0.2 6 0.5 0.6 4 0.4 6.5 4.5

解析法求解甚为不便（实际工作中的问题，不确定影响因素 更多，往往存在成百上千种可能状态，用解析法求解几乎是 这个问题虽已简化，但可能出现的状态仍有33=27种,用 解析法求解甚为不便（实际工作中的问题，不确定影响因素 更多，往往存在成百上千种可能状态，用解析法求解几乎是 不可能的），试用蒙特卡罗模拟法分析建立此产品生产线的 未来盈亏状况。 解：为了进行随机模拟，首先必须将上述三因素各种状态的 发生概率映射到一数字区间，如下表所示： 发生概率 （售价） 对应区间 （成本） （年销售量） 0.3 00~29 0.1 00~09 0.2 00~19 0.5 30~79 0.6 10~69 0.4 20~59 80~99 70~99 60~99

在均匀随机数表中，从任何一行、任何一列 中任意抽取00~99间的一个数，其出现的概率都是 均等的，都是1%。所以，从均匀随机数表中任意 抽取一个数据，该数据在00~29之间的概率为30%，再30~79 之间的概率为50%，在80~99之间的概率为20%。每抽取一次 数据称为一次模拟。如果抽取的数据在00~29之间，则认为 在各种不确定因素的影响下，未来售价应定为5元；若在 30~79之间，则售价应定为6元；若在80~99之间，则售价应 定为6.5元。另两个因素（成本、年销售量）可能发生的状态 可以用类似的方法进行模拟。按此方法进行20次数字模拟的 结果如下表：

样本号 随机数 售价 成本 年销售量 实现利润 1 16 5 75 4.5 29 4 -3 2 86 6.5 26 3.5 42 7 3 34 6 45 71 6.25 96 64 95 8.5 38 41 51 33 23 56 28 85 76 -2.75 8 10 17 48 9 07 98 08 -3.25 80 11 61 12 27 05 13 54 93 00 0.25

利润（万元）=（单价-成本）×销售量（万件）-5（万元） 14 07 5 01 2 74 4.5 8.5 15 54 6 25 3.5 56 4 16 35 27 17 68 83 1.75 18 37 98 87 19 36 89 20 86 6.5 70 3 平均利润=3.575

（1）产生（0，1）上的均匀分布的随机数，公式为： 乘同余法步骤： （1）产生（0，1）上的均匀分布的随机数，公式为： 式中：A, M为常数。这是一个递推公式，即原来的随机数Xn乘以常数A，再取以M为模的同余数，就得到新的随机数Xn+1。因此，只要给定了初值X0和常数A与M，就可以产生一个随机数列{Xn}。 （2）再用公式 产生一个（a,b）上的均 匀分布随机数列{Zn}。

模拟精度和模拟实现次数有以下一些关系： 1、如果n次模拟得到事件A发生的频率为P*， 由二项分布可知，实际概率P的样本标准差为： 在P*±2S范围内。 2、用蒙特卡罗法计算事件A的概率时，为了使最大可能 误差不大于给定值△，取95%的置信度，则由上式可解出模 拟次数不应少于n，n按下式计算： 式中：P — 事件A的概率

3、如果n次模拟得到随机变量X的算术 平均值 则随机变量X的数学期望将处在 区间的可能为95%。 为n次模拟样本均值

4、用模拟方法求随机变量X的数学期望时，如果要求模拟误差不大于给定值A，则模拟次数n应不 少于 例如：根据第一批次模拟结果取 然后再随着模拟次数的增加不断加以校正。

第四节 系统动力学 一、系统动力学概述 （一） 系统动力学及其发展 系统动力学（System Dynamics, 简称SD）是综合应用控 制论、信息论和决策论等有关理论和方法，以计算机模拟技 术为手段，研究系统动态行为的一种计算机模拟技术。系统 动力学通过建立系统动力学模型，并编制程序，进行模拟试 验，用所获得的模拟信息来分析和研究真实系统的结构和行 为，为正确的决策提供科学依据。 第二次世界大战以后，随着科学技术和工业化的进展,一

些国家存在的诸如城市人口过多，环境污染， 资源短缺等问题日趋严重。如何研究处理并解 决这些问题成为研究的焦点。 福雷斯特于1961年发表了“工业动力学”又于1969年1971 年发表了“城市动力学”和“世界动力学”等著作。并于1972年 正式提出“系统动力学的名称”。在20世纪70年代初期，西方 学术界把系统动力学赞誉为人类文化和人类进化的第三个里 程碑。在20世纪70年代以后福雷斯特教授又致力于全美国有 关社会经济问题系统动力学模型的研究，已取得了令人瞩目 的成就。

（1）应用系统动力学研究社会系统，能够容纳大量变量, （二）系统动力学的特点 系统动力学主要特点: （1）应用系统动力学研究社会系统，能够容纳大量变量, 一般可达数千以上，而这正好符合社会系统的需求； （2）系统动力学模型，既有描述系统要素之间因果关系 的结构模型，以此来认识和把握系统结构；又有专门形式表 现的数学模型，用它来实现模拟计算，从而模拟系统的动态 行为。因此，系统动力学是一种定性分析与定量分析相结合 的产物。它结合定性分析与定量分析的各自优势，是符合人 类思维发展进程的模拟技术。 （3）系统动力学的模拟试验能起到实际系统实验室的作 用。通过人和计算机的结合，既能发挥人对社会系统的了解

分析、推理、评价、创造等能力的优势， 又能利用计算机高速计算和迅速跟踪的功 能，以此来试验和剖析实际系统，从而获 得丰富而深化的信息，为选择最优决策提供有力的 依据; （4）系统动力学通过反馈回路耦合，可以处理 社会系统的复杂而高度非线性的问题; （5）系统动力学相对简单易学，易于使用; （6）虽然它是研究社会系统的得力工具，但并 不排除它也可用来解决一般工程技术中的有关问题。

(三） 系统动力学的研究对象及其特征 研究对象 主要是社会系统。社会系统包括十分广泛，凡是涉及到 人类和经济活动的系统都属于社会系统。企业，团体是属于 社会系统，经济管理系统，环境系统，资源系统，人口系统 等都也属于社会系统。社会系统一般总包含物理系统。 社会系统总是经过采集信息，并对采集的信息进行加工 处理，最后进行决策，从而导致系统行为的产生。 决策是社会系统的一大特征。 社会系统的基本特征: (1)自律性 自律性就是系统具有自主决策，自我管理，控制和约束 自身行为的能力。

二、因果关系及因果关系图 (2) 非线性性 非线性性是指社会现象中原因和结果之间呈现的极端非 线性关系。 (2) 非线性性 非线性性是指社会现象中原因和结果之间呈现的极端非 线性关系。 比如原因和结果在时间上或空间上的分离性、滞后性、 出现事件的随机性等等。这些都与信息反馈有着紧密的联系。 二、因果关系及因果关系图 系统动力学通过分析所研究系统中相互联系的元素（变 量）之间的因果关系来描述系统内部结构，构成系统的反馈 结构，从而建立模型，实现对系统的模拟。 （一）因果箭 对两个有因果关系的元素可用连接因果要素的有向边来 表示。这种有向边称为“因果箭”，箭尾始于原因要素，箭 头终于结果要素。箭头有引起、影响、导致之意。

因果关系按其影响作用的性质可以分为正因果关系和负 因果关系，正因果关系是指当由原因引起结果时，原因和结 A是原因, B是结果。可以说A导致B的发生。 因果箭的极性 因果关系按其影响作用的性质可以分为正因果关系和负 因果关系，正因果关系是指当由原因引起结果时，原因和结 果的变化方向是一致的。负的因果关系其原因和结果的变化 方向是相反的。它们分别用箭上的“+”,“—”号表示。 例如: 出生人口（原因）和总人口（结果）农业投资(原因)与 粮食产量（结果）等关系体现了正因果关系。而死亡人口与 总人口，产品出库（原因）与产品库存量（结果）之间关系 则体现了负因果关系。

（二）因果链 两个或两个以上的因果箭首尾相连（不闭合） 串联而成的因果关系称为因果链。 因果链是诸多要素因果关系具有递推性质的一种表示。 如: 要素A是要素B的原因，要素B又是要素C的原因，而要素 C又是要素D的原因。最终要素A也成为要素D的原因。 国民收入(A)增加，使得食物和营养水平(B)也提高，这样人的期望寿命(C)也相应增长，最终导致人口总数D增加。

因果链极性确定的一般方法： 1、如果因果链中所有因果箭都呈正极性， 则因果链就是正极性； 2、如果在因果链中所有的因果箭中含有偶数个 负因果箭，则因果链仍为正极性； 3、如果链中含奇数个负的因果箭，则因果链的极性就 为负极性。 结论： 因果链的极性符号与全部因果箭的极性符号的乘积相同。 因果链的极性的意义： 正极性说明起始因果箭的原因与终止箭的结果呈正因果 关系； 负极性说明它们的关系是负因果关系。

它的极性为正（偶数个负箭）。可解释成为由于人口数量（A）增加,导致资源总量（B）减少,于是又导致污染量（C）增加,最终使非健康人数（D）随之增加。人口数量（A）与非健康人数（D）呈正因果关系。这说明我们一定要注意消除污染，改善环境。

它的极性为负（奇数个负箭），可解释为，由于某种商 品销售量（A）增加，使其库存量（B）减少，这样向工厂的 订货量（C）增加，则工厂生产该商品的产量（D）也随之 增加；但工厂的生产能力有限，该种产品产量增加，必然导 致其它产品产量（E）的下降。

（三）反馈回路 若因果链首尾相接形成闭合回路则称为反馈回路。 反馈回路是自然界中和社会系统中作用与反作用普遍现 象的反映与表示。常常一些原因和结果总是相互作用的。原 因引起结果，而结果又作用于形成该原因的有关元素，从而 使原因又产生变化。这样就形成了反馈回路。 正反馈回路和负反馈回路 判别方法与因果链的判别方法相同。 性质： （1）正反馈回路的要素之间递推互促,使其要素沿着原 变化方向不断变化。即正反馈回路具有自我强化(或削弱)的 作用，是系统中促进系统发展(或衰退),进步(或退步)的因素。

例如: 由于国民收入增加使购买力增强，致使商品数量减少, 从而促使生产量增加；更进一步，生产量的增加又使国民收 例如: 由于国民收入增加使购买力增强，致使商品数量减少, 从而促使生产量增加；更进一步，生产量的增加又使国民收 入增加。国民收入增加使国民收入又进一步增加，具有自我 强化作用，是一个正反馈回路。 生产 增加 商品 减少 收入 - +

（2）负反馈回路的要素之间渐进递推，使其要素的属性 沿着原来变化相反的方向变化。因此，它具有内部调节器 （平衡器）的效果。即负反馈回路可以控制系统的发展速度 或衰减速度。是系统具有自我调节功能必不可少的环节。 粮食 亩产 土壤 肥力 负荷 + - 例如：土壤肥力增加,这样就使得粮食亩产量增加,单产增加使得土壤负荷增大,从而使得土壤肥力下降,因而起到调节，平衡的作用,它是一个负反馈回路。 -

（四）多重反馈回路 若系统结构中存在着两个或两个以上的反 馈回路，就称为多重反馈回路。 人口系统中人口的总数与年出生人数之间存在着正反馈回路，而死亡人数与人口总数之间存在着负反馈回路。年出生人数与年死亡人数共同影响着人口总数的变化。然而年出生人数和年死亡人数变化的影响因素十分复杂，它受到社会、政治、经济、环境等因素的影响，如果在继续深究下去，就会发现更多的多重反馈回路。 人口系统两重反馈回路

工业投资两重反馈回路 工业 资本 投资 折旧 + - 经济过程中，也存在着正的反馈回路与负的反馈回路。当投入一定量的工业资本，就会有一定的产出，产品的利润一部分作为投资扩大再生产。从而形成新的工业资本。所以工业资本与投资形成了互相促进的正反馈回路。而工业资本与折旧费用形成了负反馈回路，这是因为工业资本的增加，使每年的折旧费增加，折旧费增加使得工业资本要相对减少。

(1) 由人口、工业化和食物三因素形成一个正反馈回路。 (2) 工业化和自然资源，工业化和污染形成的两个负反馈回路。 世界模型基本要素因果关系图 － + 工业化 自然资源 污染 自然空间 人口 食物 占用土地 (1) 由人口、工业化和食物三因素形成一个正反馈回路。 (2) 工业化和自然资源，工业化和污染形成的两个负反馈回路。 (3) 工业化、污染和人口形成一个负反馈回路。 (4) 工业化、食物、占用土地、自然空间和人口形成一个负反馈回路。

三、系统动力学模型 (一) 实际系统的动力学描述 例：水流由塔1通过阀门2流入水箱3，在通过阀门4 流出。现假定阀门4调整到某个流量不变，那么控制 者是如何调整控制水箱3中的水位呢？ 控制过程是：控制者通过观察水箱液面而获得液面状态的信息,并与所期望的液面状态相比较,然后做出调节阀门2的决策,并通过手动调节改变阀门2的流速,使水箱液面状态发生变化,变化的液面状态信息再被观察者获取,不断循环这一过程可使水箱水位控制在一个较恒定的位置。 1.水塔 2.阀门 3.水箱 4.阀门

这里信息的传递使其成为反馈回路。其中虚线表示信息传递。因果关系图描述了调节控制过程中相关要素之间关系，从而也描述了控制的过程。 + － 水位距 目标之 差 水箱 水位 调整 流率 这里信息的传递使其成为反馈回路。其中虚线表示信息传递。因果关系图描述了调节控制过程中相关要素之间关系，从而也描述了控制的过程。 水位调整反馈回路

水塔1称为源,阀门称为”流率”(flow rate),水箱中的液面叫做“平”(level,也称为“流位”),带箭头的实线表示“流”(物质流),带箭头的虚线表示信息流。 在系统动力学中,水平、流率、流和信息是4个基本要素,它们在反馈回路中形成一个整体而发挥作用,使得真实系统的变化过程生动地加以描述,这正是系统动力学的关键。 反馈 回路 水 流 源 汇 阀4 流率 水位目标 (决策函数) 水位(状 态变量) 差值 系统动力学流程图

（二）动力学流图 常用的流图符号： （1）流（flow） （2）水平（level） （3）变化率（rate） 流表示系统中物质、能量、信息的流动。流可以是物流、资金流、信息流等,用带有各种符号的有向边描述,通常只分实物流(实线表示),信息流(虚线表示)。 （2）水平（level） 水平表示系统中要素的状态,通常是一状态变量。它是实物流的累积变量,如库存量,人口总数,水位等用矩形框表示。对它有流入流出之分,流入流出使水平变量朝相反的方向变化。 （3）变化率（rate） 变化率用来描述系统中流随时间变化的比率(速率)。例如物资的入库速率,出库速率,人口的出生率,死亡率等。速率表示决策函数,对系统状态变量即水平起着控制作用。 （1） 流 物 ￥￥ 货币 人 信息 （2）水平 L （3）速率 R

（5）辅助变量（auxiliary variable） （4）常数 （4）常数（parameter） 常数是表示系统在一次运行过程中保持不变的量, 例如:调整生产所用时间, 库 存量最大值, 期望液面高度等。常数一旦确定, 在同一次模拟中它是一个不变的常量。 （5）辅助变量（auxiliary variable） 辅助变量是为了人们易于理解以及简化变化率、变量等,人为引进的一种变量。 （6）源（source）与汇（sink） 源是指物流等的来源，汇表示物流的归宿。源取之不尽，汇永填不满。它们都是从现实中抽象出来的概念。 （7）信息（information）源 信息可来自水平，变化率等处,用带小圆型尾端的虚线箭头表示。小圆表示信息的出处，箭头表示信息的接收端。 （5）辅助变量 A （6）源与汇 （7）信息的取出

它是在输一个参数时，而这个参数具有在不同条件下取不同值的性质。这时，可以方便地用所谓的表函数来表示（即表格表示的函数）。 （8）延时 out in 1 2 3 4 1－水平变量名称 2－延时函数名称 4－延时常数 3－输出变量名 （8）延时（delay） 延时指物流或信息流有时会出现时间上的滞后。这是因为信息和物质的传递需要一定的时间才能完成。于是就会导致原因和结果，输入和输出，发送和接收等之间的滞后。滞后是造成社会系统包括物理系统非线性的根本原因。 （9）表函数 它是在输一个参数时，而这个参数具有在不同条件下取不同值的性质。这时，可以方便地用所谓的表函数来表示（即表格表示的函数）。 (10）外生变量 当系统的边界以外有变量影响时用外生变量表示。它用双线圆表示。 （9）表函数 表函数 （10）外生变量

例： 现研究一个经营单一商品的零售店的订货策略问题，要 求应用系统动力学模型进行模拟，以选择最优订货策略。 零售商店应考虑：商店的销售量,这是问题的因;商店的库存量,这是一个累积变量;还有商店的订货量,它要引起库存增加,并促使工厂生产,是速率变量。 工厂应考虑:工厂未供订货量,即商店向工厂订货,但工厂未能供应的数量,它是一个累积变量,可促使工厂扩大生产;工厂的生产能力,工厂的计划产量等。

商店 订货 工厂未 供货量 计划 产量 工厂生 产能力 工厂 库存 销售 购货 顾客数 + - 因果关系图及反馈回路 系统边界

例如：商店订货量增加，会使工厂未供货量增加， 计划产量增加，这就要求工厂扩大生产能力，从 而使产量增加。工厂产量增加会使商店库存增加， 从而又使订货减少。整个是一个负极性反馈回路。 有两种实体流。即商品流与订货流。商店库房里积累而 形成的库存量（L2），工厂积累而形成的未供货量（L1）， 它们是水平变量。另外还可看出，影响水平变量的速率变量 有三个，它们是商店订货率（R1），工厂生产率（R2），以 及商店的销售速率（R3）。而工厂生产能力和计划产量属于 辅助变量，分别用A1 及A2 表示。最后顾客购买量属于外生 变量，用 D来表示。

工厂子系统的流图 商店流图

系 统 整 体 流 图 S1—平均销售量,它是一辅助变量, 顾客购货数D是外生变量, D1(周)—是调整生产的时间，D2 —完成未供订货所需的时间，D3 —商店平均订货时间。

(三）流图的DYNAMO表达 所谓DYNAMO是dynamic model的缩写，是 动力学模型之意。 1、DYNAMO 的时间规定 态变量是连续的而且是对时间一阶可导的。 图中： J—过去时刻 K—现在时刻 L—未来时刻 JK—由过去时刻到现在时刻的时间间隔 KL—由现在时刻到未来时刻的时间间隔 时间标记示意图 时间 JK KL J K L

系统动力学采用逐步（step by step） 模拟的方法。模拟时间间隔(或称为步长), 就是JK, KL 的时间间隔，用单位时间DT (delta T) 来 表示。DT的单位可以取年，月，日，周，时，分， 秒，视模拟精度而定。必要时可以取得很小，以便 很好地逼近连续时间系统。由此可见，时间步长DT 要选择适当。一般DT的选择可以根据经验或在计算 机上进行运行试验，然后比较结果从而确定出合适 的 DT值。

DYNAMO中的数采用十进制表示,最大可用8位。 “十”、“一”号和小数点也各算一位。超过8位用浮点数表示。 DYNAMO中的所有的量，不是常量就是变量。常量在 模拟的全过程中保持不变。而变量都带有时间下标，如水平 变量WATER.K，WATER.J等。 DYNAMO中变量名的字符数不得超过6个，且第一个字 符必须为字母。任何具有字符数1—6的变量都能被识别。但 是，在DYNAMO已经明确定义了字符串的名称（也称为保 留字符）必须除外。也就是说对保留字符用户不能再用作变 量名，这些保留字符主要有L , R , A , S , N , C , T , PRINT , PLOT , RUN , NOISE , NOTE等。

3、运算符号 DYNAMO 采用通用代数符号,加法“＋”， 减法“－”，乘法“*”，除法“/”，可用括号。 运算顺序为先乘方、开方，再乘除，最后是加 减。括号最优先，此外同级运算总是先左后右。 4、DYNAMO语句 DYNAMO与其它计算机语言一样，基本单元 是语句。共有15种语句。

DYNAMO 语句一览表 方程语句 控制语句 其他语句 标识符 语句名称 L 状态方程式 SPEC 说明 * 标题 R 决策方程式 OPT 选择 NOTE 注释 A 辅助方程式 PRINT 制表 X 续行 S 附加方程式 PLOT 绘图 C 常数方程式 RUN 运行 N 初始值方程式 T 表方程

语句由标识符和语句体构成。标识符与语句 体之间至少留有一个空格。语句体由变量、数值、 算式、字符串等按规定书写而成。语句体中不能 有空格出现。每行有的版本规定不能超过72个字符，不够书 写时可另起一行，但左端必须加上续行符号“X”。 （1） 状态方程语句 例如：有顾客向商店购买商品，商店又向工厂定购产品。工 厂向商店交付商品，商店向顾客销售商品的过程中,在时刻 t 的商品库存量 y 显然与单位时间商品入库量 和出库量 的关系满足：

L Y . K = Y . J + DT * ( XIN . JK －XOUT . JK ) 用DYNAMO语言描述，就是现在时刻K的库存量等于过 去时刻J的库存量,加上由过去时刻J到现在时刻K的入库量与 出库量之差乘以单位时间DT 。即为状态方程语句： L Y . K = Y . J + DT * ( XIN . JK －XOUT . JK ) 状态方程语句的一般形式: L LEVEL . K = LEVEL . J + DT * ( INFLOW . JK－OUTFLOW .JK ) K时刻的水平等于J时刻的水平加上单位时间（即模拟步长） DT乘以JK 期间输入流量与输出流量之差。

L 语句中,水平变量必须由初值语句给赋予 初值，DT 可用SPEC说明语句或者常数方程式语 句给出。 （2）速率方程语句（亦称决策方程语句） 速率方程语句是用来计算速率变量的方程式，是描述状 态方程中水平变量在单位时间DT内增加或减少的量。 速率变量是一种决策变量，而决策的多样性，使速率变 量的速率方程语句没有固定的形式，要依据具体情况而定。 但一般说来它是水平变量、参数、常数等的函数，即 R RATE.KL = F(水平变量，参数，常量等) 例如: 设KL期间工厂的生产率XOUT.KL 与在时刻K的未供订 货量成正比，这是一个在生产能力一定的情况下比较合理的 假设。这种情况下工厂生产率的速率方程语句为：

R XOUT . KL = Z. K/C 其中：R—表识符，后面是语句体。 C—常数，它表示满足未供订货量的时间，是根据经 验设定的（当未供订货使生产能力改变后，可选C 取不同的 常数值。因为生产能力提高必然要缩短满足未订货的时间）。 Z.K—K时刻的未供订货量。 商店的定货速率方程语句可写为： R ONR.KL = AOR.K + (DIN.K-AIN.K)/AJT 在 K L期间的订货量等于 K 时刻的平均销售量与 K 时刻 的期望库存与实际库存之差除以实际库存调整到期望库存时 间之和。

一般速率变量都是由现在K时刻的某些变量，水平变量 和常数来求出将来（KL时间间隔）的速率变量。确定它们之 间的函数关系是极端重要，但很困难的工作。一般可以根据 历史经验，理论分析，逻辑分析等方法来获得。 （3） 辅助方程语句 辅助方程语句是计算辅助变量的DYMAMO方程式。一 般如果在速率方程式比较复杂或者DYMAMO语言书写所不 允许的情况下，则可以引入辅助变量和辅助方程语句加以解 决。辅助方程语句一般格式为： A 辅助变量. K = 算式，变量或数值 其中：A—标识符。辅助方程语句计算K时刻的辅助变量值。

例如：上述订货速率方程语句中DIN.K与AIN.K 等是辅助变量。辅助方程语句可表示为： A AOR.K = SMOOTH(ARR.JK，DRR) A DIN.K = AIR*AOR.K 式中：ARR—顾客向商店的订货速率。 DRR—平滑时间。 SMOOTH—一次指数平滑函数，其作用是对ARR.JK 数据进行一次指数平滑。 AIR—商店的存货常数，它表示期望库存与最近平滑 平均需求速率间的关系。

（4） 附加（Supplementary） 例如： S TOTAL.K = AIN.K+IAF.K 式中：S—标识符。 TOTAL.K—商品总量。 AIN.K—商店库存量。 IAF.K—工厂库存量。 （5） 常量方程语句 常量方程语句是给常量赋以常数值的语句。常量在一次 模拟运行中保持不变，但在不同次的运行中可以采用不同的 值。

例如： C T1 = 0, T2 = 0, T3 = 100 或 C T1 = 0/T2 = 0/T3 = 0 式中：C—标识符，语句体可以为不同变量赋值，但要用“，” 或“/”分开。 （6） 赋初值方程语句 赋初值语句是运行开始时为变量赋以初始值的语句。在 模拟开始时，所有水平变量及部分决策变量和辅助变量需要 具有初值。例如： L Y.K = Y.J + DT*(XIN.JK-XOUT.JK) N Y = 1000 式中：N—标识符，一般这种初值语句紧跟在需赋初值的状态 方程等语句后面。

（7） 表函数语句 表函数语句一般由两个语句来构成。一条辅助方程语句, 一条表方程语句。这种语句用来设置一种能描述两变量间的 非线性函数关系，即，某辅助变量是另一个称为输入变量的 非线性函数。这个辅助变量称为表函数，表函数的一般表示 形式为： A 表函数变量名 .K = TABLE（表变量名，输入变量名. K , X, Y, Z） T 表变量名 = W0/W1/…/Wn， 式中：TABLE—表函数符号，输入变量是表函数的自变量， 表函数变量是表函数的因变量。

应，并取相应的数值，X , Y , Z变量表示输入变 量的起始值， 终止值以及变量取值的等分间隔。 T—表变量名赋值。 辅助方程语句实现自变量与因变量之间的对 应，并取相应的数值，X , Y , Z变量表示输入变 量的起始值， 终止值以及变量取值的等分间隔。 例如：变量RDEN代表野兔密度，RNBF代表野兔纯出生率 （只/年），这两变量之间呈非线性关系。 RDEN 0.25 0.5 0.75 1.0 1.25 RNBF 1.50 2.40 2.20 1.10 0.00 -1.00

A RNBF.K = TABLE ( TRNBF , RDEN. K, 0, 1.25, 0.25) 在DYNAMO中可用表函数语句实现变量RNBF与RDEN 之间的非线性关系。 A RNBF.K = TABLE ( TRNBF , RDEN. K, 0, 1.25, 0.25) T TRNBF = 1.50/2.40/2.20/1.10/0.00/-1.00 时间标记表 左端变量类型 左端标号 右端变量 L A R S C N K J JK 不可 无 KL

(8) 控制语句 控制语句是用来控制模拟,打印数据,图形输出的语句。 ① SPEC说明语句 SPEC说明语句包含四个参数： DT—模拟步长，即两次模拟的时间间隔。 LENGTH—用来设置模拟的终止时刻。在DYNAMO中 除非特别说明，总是假定时间的起始值为0。 PLTPER—给定输出图形中相邻点的时间间隔。如不 声明DYNAMO默认为零。 PRTPER—给定以表格形式输出打印结果中相邻两数 据间的时间间隔。同样如不特别声明自动取零。

例如： SPEC DT=0.125/LENGTH=30/PRTPER=0/PLTPER=0.5 该语句表示，模拟步长为0.125个时间单位，模拟时间长 度为30个时间单位，不要求输出打印结果，输出图形相邻点 之间的间隔为0.5个时间单位。 ② PRINT制表语句 该语句用来将要输出的变量以表格的形式打印显示出来。 如： PRINT AOR , DIN。 计算机执行该语句，将把各时刻的平均销售量(AOR)与 各时刻的期望库存量按SPEC说明语句中参数规定，输出到打 印机或显示器。

PLOT ABC = A , DEF = D/LMN = L(0,*)/TUV = T , XYZ = 同比例作图的要用“，”号分开。用户要用特定符号来画图， 要用“=特定字符”的形式。如： PLOT ABC = A , DEF = D/LMN = L(0,*)/TUV = T , XYZ = X(0 , 1000) 计算机执行该语句，ABC，DEF两变量将按同一比作图, 并且ABC变量的曲线用字符“A”, DEF用字符“D”画出。LMN 将从下限0到DYNAMO选择的上限作图。而TUV , XYZ将以 相同的比例在0到1000的范围内作图,并分别用“T”和“X”字符 作出。

④ OPT选择语句 选择语句用来选择输出方式，选择屏幕上 坐标，图中曲线点迹的显示方式等。 ⑤ RUN运行语句 运行语句用于模拟的启动。或改变参数后再开始运行。 它一般放在全部程序的最后。 (9) 其它语句 标释语句 如 * BANK BALANCE MODEL, 是 DYNAMO程序的名 称及有关信息。它一般要放在全部模拟程序的最开始一行。 注释语句

对程序的说明，帮助人们理解程序，它是非执行语句。 (10) DYNAMO函数 这是为了模拟的方便，DYNAMO中内置了21种函数。 如，NOTE FACTORY SECTION, 其主要作用是用来 对程序的说明，帮助人们理解程序，它是非执行语句。 (10) DYNAMO函数 这是为了模拟的方便，DYNAMO中内置了21种函数。 固有函数表 函数类别 函数名称 形式或说明 普 通 函 数 平方根函数 SQRT（变量或数值） 指数函数 EXP 自然对数函数 LOGN 正弦函数 SIN 余弦函数 COS

PULSE（脉冲高，脉冲产生时刻，脉冲间隔） 阶跃函数 STEP（函数值，发生时刻） 斜坡函数 普 通 函 数 脉冲函数 PULSE（脉冲高，脉冲产生时刻，脉冲间隔） 阶跃函数 STEP（函数值，发生时刻） 斜坡函数 RAMP(p,q) 时间≤q , RAMP = 0（p斜率） 随机函数 NOISE（）产生（-1/2，1/2）区间内的均匀分布随机数 采样函数 SAMPLE（p,q,r） 最大值函数 MAX(p1 ,p2 ,…,pn) n个变量中取最大值 最小值函数 MIN(p1 ,p2 ,…,pn) n个变量中取最小值 选择函数 CLIP相当于FOTRAN中的IF 表格函数 TABLE（表格变量，变量，最小值，最大值，间隔） 合计函数 SUM1 SUM2 SUM3

宏 函 数 SMOOTH 一次指数平滑序列 DLINF3 信息流的延迟 DELAY1 DELAY2 物流，资金流，订货流的延迟 DELAYP 订货流，物资流资金流的延迟数量

5、有关说明 建立DYNAMO模拟程序只要依照流图,将各变量 用DYNAMO各种语句加以合理表达，并加上所需的 控制语句,说明语句，运行语句等即可获得到可以模拟运行的 模型。　 注意： DYNAMO模拟模型各种语句要在确定系统的水平变量、 速率变量、辅助变量、常数等的基础上，在确定这些变量之 间的关系的基础上才能进行。只有这样才能编写出正确语句, 正确的模型。另外,DYNAMO最大的特点是面向方程对象的 语言。除了个别语句外，编程者不需考虑执行顺序。因此,程 序书写简单，再加上输出图非常容易，因此，使用极其方便。

四、几个典型结构的DYNAMO模拟计算 目的: 一是通过几个典型结构进一步熟习模拟模型的建立。 三是了解三个基本结构的动态行为。 1、一阶正反馈回路 出生率R1 (人/年)增加, 总人口P增加; 总人口增加, 又使得 出生人口增加(出生率增加)。这就是说人口和出生率之间形成 了正的反馈回路。 给定人口的年增长率是2%，人口的初始值是100，则模拟 人口增长过程的DYNAMO程序是：

* POPULATION INCREASING SPEC DT = 1/LENGTH = 30/PRTRER = 1/PLTPER = 1 + 出生率R1 人口数 P R1 一阶正反馈回路流图 * POPULATION INCREASING SPEC DT = 1/LENGTH = 30/PRTRER = 1/PLTPER = 1 L P.K = P.J + DT * ( R1.JK－0 ) N P = 100 R R1 . KL=P . K * C1

PRINT P , R1 PLOT P = 1 RUN 模拟表 模拟结果示意图 模拟步长 (年) P R1 100 2 1 102 2.04 P（人口数） t / 年 模拟结果示意图 模拟表 模拟步长 (年) P R1 100 2 1 102 2.04 104.04 2.0808 …

设：X 表示初始库存量, Y 表示期望库存量 (假设期望库 存量大于初始库存量), Z 表示将当前库存量调整到期望库存 的时间。 2、一阶负反馈回路 设：X 表示初始库存量, Y 表示期望库存量 (假设期望库 存量大于初始库存量), Z 表示将当前库存量调整到期望库存 的时间。 + － 订货速率 R1 Z Y X R1 I D 一阶负反馈回路流图 库存差D 库存量I

该反馈回路将库存调整到期望库存量的机理如下： 当库存量增加，库存量与期望库存量的差额D就减少, 即 两者是负因果关系。于是D-R1-I-D就构成了负反馈回路。 现假定：X=1000,Y=6000, Z=5 (周) ,则描述该库存系统 动态行为的DYNAMO程序为： * One-order Negative Feedback System Spec DT = 1/LENTH = 20/PRTRER = 1/PLTPER = 1 L I.K = I.J + DT * (R1.JK－0) N I = X C X = 1000 R R1.KL = D.K/Z A D.K = Y – I.K C Y = 6000

C Z = 5 PRINT I , R1 , D PLOT I = 1 负反馈结果示意图 模拟步长(周) I R1 D 1000 5000 数量（件） D I T/周 6000 1000 负反馈结果示意图 模拟步长(周) I R1 D 1000 5000 1 2000 800 4000 2 2800 640 3200 3 3440 512 2560 4 3952 409 2048 ┅

3、两阶负反馈回路 由于从订货到入库有时间延迟，因而形成所谓“途中存 货”。这样，库存系统就会从原来的一阶负反馈回路变为两 阶负反馈回路。 二阶负反馈回路示意图

现假定初始库存量X=1000, 期望库存量Y=6000,调整库存 时间Z=5(周), 初始途中存货G=10,000, 定货商品的入库时间 W=10 (周),则描述该库存系统动态行为的DYNAMO程序是： * Second Order Negative Feedback SPEC DT = 1/LENGTH = 20/PLTPER = 1/PRTPER = 1 L G.K = G. J + (DT) * (R1.JK-R2.JK) N G = V C V = 10000 R R1.KL = D.K/Z A D.K = Y－I.K

C Z = 5 C Y = 6000 R R2.KL = G.K/W C W = 10 L I.K = I.J + (DJ) * (R2.JK-0) N I = X C X = 1000 PLOT I PRINT G , R1 , D , R2 , I RUN

步长 （周） G R1 D R2 I 10000 1000 5000 1 800 4000 2000 2 9800 600 3000 980 3 9420 404 2020 942 3980 4 8882 215.6 1078 888.2 4922 … T/周 1000 6000 数量（件） 模拟结果示意图图

五、系统动力学模拟的基本步骤 从上面的内容，特别是系统动力学方法解决实际系统问 题的全过程，我们可以看出以下两点： a.系统动力学提供了对社会系统进行模拟试验的有效手 段,它是一种计算方法。这种方法使用虽然不太复杂，但要取 得满意的结果，也需要花费相当的努力，需要多方面的人共 同合作。 b.系统动力学解决问题应遵循下列步骤： （一）明确系统模拟的目的 一般来说，系统动力学对社会系统进行模拟试验的主要 目的是认识和预测系统的结构和未来的行为，以便为进一步 确定系统结构和设计系统的最佳参数，以及制定合理的政策 等提供依据。

（二）确定系统边界 系统动力学建模的对象是封闭系统。必须根据系统模拟 的目的，确定系统的合理边界。这是因为系统动力学分析的 系统行为是基于系统内部种种因素而产生的，并假定系统外 部因素不从本质上影响系统行为，而且，外部因素也不受内 部因素的控制。这个假定实际中是容易满足的，因为只要把 影响系统大的本质因素都纳入系统后（通过各种方法包括常 数，函数表示的影响等) ,系统就是封闭系统。因此系统动 力学并不是只能用于封闭系统。 （三）因果关系分析 在边界分明的基础之上，就要明确系统内部各要素间的 因果关系，并用表示因果关系的反馈回路来描述形成因果关 系图。

要做到这一点，首先要求系统分析人员有丰富的实践经 验，对实际系统有敏锐的洞察力，还要要求集中各方面专家 的群体智慧。只有这样才能较为正确地制定各要素间的因果 关系即反馈回路。正是这些反馈回路的合理耦合使系统的复 杂行为得以恰当表现。 （四）建立系统动力学模型 系统动力学模型可包括两部分，即： 1、流图 流图是根据因果关系图，特别是其中的反馈回路，应用 专门为系统动力学规定的一套描述各种变量的符号绘制而成 。 在绘制的过程中，对因果关系及各变量间的关系是加深认识 过程，以及不断细化和补充的过程。

经过这样的过程，系统会被考虑的更加细致,更加完善， 更接近于模拟的现实，流图是一种图像模型，她能帮助我们 掌握复杂社会系统的结构及其行为的动态性，也便于人们对 系统进行讨论与沟通。 2、结构方程式 结构方程式是系统中各变量间定量关系的反映，是模拟 模型的基础。要对变量间的关系进一步过细研究才能完成这 些方程式的DYNAMO语言形式的正确表达。因此这一步也 是对以前不断完善与提高的过程。经过这一步，再加上一些 控制语句等就可得到模拟模型（即程序）。DYNAMO模拟 程序必须通过系统中的水平变量，速率变量，常数，辅助变 量等基础上来进行。所以分析各变量之间关系，确定它们之 间的函数关系是一个重要任务。

（五）计算机仿真试验 根据DYNAMO模拟模型在计算机上进 行模拟计算。 （六）结果分析 为了解模拟试验是否达到预期目的, 或者检验系 统结构是否有缺陷，必须要对模拟结果进行分析。 （七）系统模型的修正 根据模拟结果分析，对系统模型进行修正。修 正的内容包括：修正系统的结构，或修正系统的运 行参数、策略，或重新确定系统边界等等，以便最 终使模型能更真实反映实际系统。

例：野兔—山猫模型 根据学者威廉，沙福特的工作改编。某自然生态区，野 兔山猫共存，最初的山猫数为6只，野兔数为700只，而该系 统容纳野兔的最大量是2400只。某管理机构要想了解野兔山 猫的相互依存关系以及它们数量的动态变化规律。拟定采用 系统动力学模拟来达到其目的。通过分析，确定了该系统的 边界，即系统只包括野兔山猫两个主要因素，还有与其有关 的一些要素。这两要素的依存关系为：山猫以捕食野兔为生 存的食物来源。如果山猫数增多，则每年捕食的野兔数增 多，野兔的数量就会变得越来越少。这样因为食物减少，山 猫数就会又减少。随着山猫数的减少，野兔就又将逐渐增多。 野兔的增多不是无限的，而是受环境容纳量的限制。除此自 然相互依存，相互制约下，考虑到人们狩猎的需要，所以把 人们对野兔和山猫的猎捕因素也纳入封闭模型考虑。从而可 得该系统的因果系统图。

野兔—山猫模型的流图 FRABTR RTRAP RNBR RABBIT RKL CC RDEN TRNBF RRPLS FLNXTR LTRAP LNBR FLNXTR RRPLS RNBR RTRAP FRABTR TLNBF TRNBF RDEN TRKPL LSUBR LYNX RABBIT RKL

参数的含义如下： RABBIT— 野兔数（只），水平变量（只）； RNBR— 野兔纯出生率, （只/年）； INRAB— 模拟开始时的野兔数（只），本例为700只； TRNBF— 野兔纯出生率因子，由表函数给定 RDEN— 野兔密度 CC— 系统容纳野兔的能力（只），本例为2400只; RKL— 山猫捕食野兔的速率（只/年），速率变量； RTRAP— 野兔被人猎捕率（只/年）; FRABTR— 野兔被人猎捕系数（1/年）; TRKPL— 每头山猫一年内捕食的野兔数 (只/头·年,由表函 数给出); LYNX— 山猫数（头）, 水平变量; INLYNX— 模拟开始时的山猫数（头）, 本例为6头;

LNBR— 山猫纯出生率（头/年）; LTRAP— 山猫被人猎捕率（头/年）; LSUBR— 山猫的生存状况因子; RRPLS— 维持山猫生存年需捕食野兔数(只/头·年), 本例为 200只; FLNXTR—山猫被人猎捕系数（1/年）; 根据上述流图对各变量间的关系进行仔细研究、分析,可 得出各变量间的关系如下： 对于野兔： （1）野兔数.K = 野兔数.J + (DT) × (野兔出生率－猫捕食野 兔率－野兔被人猎捕率) 。 （2）初始野兔数700只。 （3）野兔出生率.K =野兔数.K × 出生率因子；

（4）猫捕野兔速率.KL = 山猫数.K × 每只猫一年猎捕野兔数 ； 野兔出生率因子与密度关系表 野兔密度 0.25 0.5 0.75 1 1.25 野兔出生率因子 1.50 2.40 2.20 1.10 0.00 -1.00 （4）猫捕野兔速率.KL = 山猫数.K × 每只猫一年猎捕野兔数 ； 每只猫年捕野兔数与野兔密度关系表 野兔密度 0．25 0．.5 0．75 1 1．25 每只猫年捕野兔数 0．00 150 250 325 375 400

（6）野兔被人猎捕数. KL = 野兔数.K×野兔被人猎捕系数。 该例计算为了简单，假定猎捕系数为 0，于是没有研究 （5）野兔密度.K = 野兔数.K/CC ；CC = 2400 （6）野兔被人猎捕数. KL = 野兔数.K×野兔被人猎捕系数。 该例计算为了简单，假定猎捕系数为 0，于是没有研究 得出该猎捕系数变化规律。 对于山猫： （7）山猫数.K = 山猫数.J + (DT)×（山猫纯出生率—山猫 被人猎捕率）；山猫数为6只作为初始值开始模拟。 （8）山猫纯出生率.KL = 山猫数.K×山猫出生率因子； 山猫出生率因子与山猫的生存状况的关系表 生存状况因子 0．5 1 1．5 2 山猫出生率因子 - 4．0 - 0．6 0．0 0．3

RABBIT-LYNX MODEL 少捕野兔数（本例为200只） (10) 山猫被人猎捕率.KL = 山猫数.K×山猫被人猎捕系数。 （9）而生存状况因子.K = 每只山猫一年捕野兔数.K/生存至 少捕野兔数（本例为200只） (10) 山猫被人猎捕率.KL = 山猫数.K×山猫被人猎捕系数。 在该模拟中，为了简单，假定山猫被人的猎捕系数为0。 有了各种方程、系数等，并考虑到要绘图输出山猫数 （用“L”字符画出），野兔数量（用“R”字符画出）的变化曲 线。还有要设定模拟的步长，总时间以及模拟总长度等，可 得下列模拟程序。 RABBIT-LYNX MODEL NOTE NOTE …………………………………………………………… NOTE 　NAME – RABBIT –LYNX MODEL NOTE 　AUTHOR –WILLIAM SHAFFER

NOTE ………………………………………………… NOTE RABBIT SECTOR L RABBIT．K = RABBIT．J + (DT) (RNBR． JK - RKL． X JK －RTRAP．JK) N RABBIT = INRAB NOTE RABBITS ( RABBITS ) C INRAB = 700 NOTE INTITAL RABBITS ( RABBITS ) R RNBR．KL = (RABBIT．K)(RNBF．K) NOTE RABBIT NET BIRTHS (RABBITS / YEAR) A RNBF. K = TABLE ( TRNBF, RDEN. K , 0 , 1.25 , 0.25 ) NOTE RABBIT NET BIRTHS FACTOR( 1 / YEAR) T TRNBF =1.50/2.40/2.20/1.10/0.00/-1.00

NOTE TABLE FOR RABBIT NET BIRTH FACTOR R 　 RKL.KL = (LYNX.K)(RKPL.K) NOTE RABBIRS KILLED BY LYNX (RABBITS / YEAR ) A 　　 RKPL.K = TABLE ( TRKPL, RDEN.K , 0 , 1.25 , 0.25 ) NOTE RABBIRS KILLED PER LYNX (RABBITS / LYNX – YEAR ) T 　　 TRKPL=000/150/250/325/375/400 NOTE TABLE FOR　 RABBITS KILLED PER LYNX A 　 　RDEN.K= RABBIT．K/CC NOTE RABBIT DENSITY (DIMENSIONLESS) C 　　 CC = 2400 NOTE CARRYING CAPACITY FOR RABBITS (RABBITS)

R 　 RTRAP．KL = (RABBIT．K)(FRABTR) NOTE RABBITS　 TRAPPED (RABBITS / YEAR) C 　 FRABTR = 0.0 NOTE FRACTION OF RABBITS TRAPPED ( 1/YEAR) NOTE NOTE LYNX SECTOR L LYNX．K =LYNX．J +(DT)(LNBR.．JK － LTRAP．JK) N LYNX = INLYNX NOTE LYNX (LYNX) C 　 INLYNX =6 NOTE INITIAL LYNX (LYNX) R 　 LNBR．KL = (LYNX.K)(LNBF.K)

NOTE LYNX NET BIRTHS (LYNX /YEAR) A 　 LNBF．K = TABLE (TLNBF．LSUBR．K, 0, 2, 0.5) NOTE LYNX NET BIRTH FACTOR (1 /YEAR) T TLNBF　＝　－4.0 /－0.6 /0.0 /0.3 /0.5 NOTE TABLE FOR LYNX NET BIRTH FACTOR A LSUBR．K =RKPL．K /RRPLS NOTE LYNX SUBSISTENCE RATIO (DIMENSIONLESS) C 　 RRPLS=200 NOTE REQUERED RABBITS PER LYNX FOR SUBSISTENCE NOTE (RABBITS /LYNX－YEAR)

R LTRAP．KL =(LYNX.．K)(FLNXTR) NOTE 　LYNX TRAPPED (LYNX /YEAR) C 　FLNXTR=0.0 NOTE FRACTION OF LYNX TRAPPED ( 1/ YEAR) NOTE NOTE CONTROL STATEMENTS SPEC 　 DT=0.125 /LENGTH=30 /PLTPER=0.5 PLOT 　 RABBIT =R /LYNX =L(5, 25) OPT 　 TXI = 20 RUN

六、结果分析 结果分析是模型能否反映实际系统的行为，模拟是否已达到 预期目的检验阶段或验收阶段，必须给予足够的重视。 系统动力学模型的修正是在结果分析的基础上进行的。 结果分析是模型能否反映实际系统的行为，模拟是否已达到 预期目的检验阶段或验收阶段，必须给予足够的重视。 除了上面已经讲过的数据的输入分析，输出分析中应采 取的必要措施外，还应更注意模型的有效性分析。 模型的有效性主要是指模型自身和它的行为与实际系统 和它的行为进行比较能达到一定的精度，能满足模拟的目的。 模型的有效性分析是指模型是否有效的整个分析修正过程。

就是要由建模者与熟悉所研究系统的专家对模型和它的输出作出模型是否正确，是否满足模拟的要求的判断。 有效性的逐渐逼近过程 修改 初始模型 模型第一次修改 模型第二次修改 模型与实际比较 修改后 第二次修改后 现实系统 在这个过程中，要注意两方面的检验： 一是主观检验 就是要由建模者与熟悉所研究系统的专家对模型和它的输出作出模型是否正确，是否满足模拟的要求的判断。

二是客观检验 就是要将所研究系统的有关数据和相应的模型输出数据 进行比较。比较中必然要用到统计检验，以及有关各种各样 的定性分析和定量分析的工具。借助它们做出模型是否有效 的判断。 特别要注意下列三方面的有效性： (1) 直观有效性 就是模型对用户是否具有吸引力，要使用户相信模拟模 型对他们决策的重要性，于是要对各种用户所关心的输入参 数对输出结果的灵敏度进行试验。看是否满足用户的要求。 (2) 模型假设的有效性 包括结构假设和数据假设。结构假设包括实际系统工作 的机制，建模时所用的各种简化与抽象。这些假设应通过

实际观察验证，或与实际了解系统的专家讨论确定。 数据假设是建立在可靠数据的收集和统计分析的基 础之上的，如：山猫-野兔模型中所用的表函数，还 有系统对应总容量等。要就数据的可靠性与实际专 家及有关人员讨论加以确认。 (3) 输入-输出变换的有效性 这是指模型是否能预测系统的未来。模型是输 入-输出的变换器。就是说模型接受输入参数值，把 输入变换成输出，所以输入输出变换的有效性要求 必须至少存在一组输入条件下的系统输出数据，以 便与模型输出比较，确认模型的有效性。