软件包中系统动力学函数.doc

第8章 Vensim PLE 软件包中系统动力学函数

系统动力学所以能处理复杂的系统问题，除提出流位流率系简化流率基本入树建模法去描述系统外，还有一个重要原因是其专用软件都设计了一系列通用的系统动力学函数。

第一节 数学、逻辑、测试函数

§ 8.1.1 数学函数
Vensim PLE备有五种普通数学函数供用户使用。
1．SIN（X）
定义1：SIN（X）为三角正弦函数，X须以弧度表示，其值小于8.35×105
当自变量是角度时,应通过乘以2π/360 转化为弧度。
2．EXP（X）
定义2：EXP(X) = eX ，e是自然对数的底，e=2.7182…，X的值必须小于36。
人们常用指数函数去描述系统，有了上面函数将会带来很大方便。
3. LN（X），变量X大于零。
即以e为底的对数函数，它与EXP（X）互为反函数，这样可以用EXP（X）和LN（X） 来计算非以e为底的幂函数和对数函数。
4. SQRT（X）=√X—，X必须是非负量。
5. ABS（X） = │X│，对X取绝对值。

§ 8.1.2 逻辑函数
逻辑函数的作用类似于其它计算机语言中的条件语句，Vensim PLE的逻辑函数有三种。
1. 最大函数MAX（P，Q）
MAX表示从两个量中选取较大者，P和Q是被比较的两个量，结果也是在这两个量中选取。
P 若P≥Q
定义1：若MAX（P，Q）=
Q 若P≤Q
其中P，Q是变量或常量，则MAX（P，Q）为最大函数。
可用MAX函数从多个量中选取较大者。如从P，Q，D三个量中选择较大者可用：MAX（D，MAX（P，Q））。
最小函数 Q 若P≥Q
定义2：若MIN（P，Q）=
P 若P≤Q
则MIN（P，Q）为最小函数。
１． MIN同MAX一样，可以从MIN（P，Q）
基本功能中派生出各种用法。
3. 选择函数IF THEN ELSE(C,T,F)
定义3：若IF THEN ELSE(C,T,F)
T C条件为真时
= （C为逻辑表达式）
F 否则
则IF THEN ELSE(C,T,F)为选择函数。
IF THEN ELSE函数常用于仿真过程中作政策切换或变量选择。有时也叫条件函数。
§ 8.1.3 测试函数
设计这一部分函数的目的主要是用于测试系统动力学模型性能用，所以称为测试函数。
在给出测试函数以前，我们必须重申一个概念，系统动力学的变量皆是时间TIME的函数，所以当仿真时间TIME发生变化时，各变量值都随之发生变化。不过，各变量与TIME的依赖关系存在差别，有的是以TIME为直接自变量，有的则是间接变量。测试函数以TIME为直接自变量，但在函数符号中常缺省。
1. 阶跃函数STEP（P，Q）
定义1：
Q 若TIME≤Q
STEP（P，Q）=
P 若TIME > Q
其中，P---阶跃幅度；Q---STEP从零值阶跃变化到P值的时间，则STEP（P，Q）为阶跃函数。

P

STEP

A*(C-B)

A

*(TIME-B)
跃 Q B TIME
图8.1.1

阶路STEP（P，Q）和斜坡RAMP（A，B，C）函数图
2. 斜坡函数RAMP（P，Q，R）
定义2：
0 若TIME≤Q

RAMP（P，Q，R）= P *（TIME-Q），若R≥TIME > Q

P *（R-Q），若TIME≥R
其中,P为斜坡斜率，Q为斜坡起始时间，R为斜坡结束时间，则RAMP（P，Q，R）为斜坡函数。
3. 脉冲函数PULSE(Q,R)
定义3: 若PULSE(Q,R)随TIME变化产生脉冲。
其中：
Q---第一个脉冲出现的时间
R---相邻两个脉冲的时间间隔
脉冲宽度为仿真步长，则PULSE(Q,R)为脉冲函数。
4．均匀分布随机函数RANDOM UNIFORM（A，B，S）
定义5：RANDOM UNIFORM（A，B，S）产生在区间（A，B）内的均匀分布随机数，S给定随机数序列就确定，S取不同的值产生随机数序列也不同。RANDOM UNIFORM（A，B，S）为均匀分布随机函数。
上面我们给出了四种测试函数，实际上还有前面数学函数SIN（X）等也可以作为测试函数。
一个系统动力学模型，可以通过改变常数再运行的办法，实现多种测试函数分别进行测试。

第二节 表函数

§ 8.2.1 表函数表示形式
1. Vensim PLE中表函数表示形式
定义1 自变量与因变量的关系通过列表给出的函数叫表函数。例如下表就确定了一个表函数。

自变量X
0
1
1.5
2
2.5
因变量Y
0.5
1 1
2
5
10

表函数是系统动力学的一个重要特征,它用于建立两个变量之间的非线性关系,特别是软变量之间的关系。例如：员工士气对工作效率的影响程度。一般，两个变量先归一化或者先规整化，再根据经验给出大致的关系图来。这样设计的变量是无量纲量。在进入Vensim PLE 软件Equation Editor，即点去图标Y=X2 后，若方程还未定义，有AS Graph选项。选择此选项，会出现下面对话框。该对话框用于图形化定义，上例表函数可直接填入框中(图8.2.1)。包括自变量和函数值即因变量值列举，自变量和函数的最大值等。当自变量为非已知统计点时，可用线性插值法取其近似值。用鼠标左键在图形框中点按，会自动构成图形。
Vensim PLE软件中表函数表达形式还可通过选择方程类型TYPE中Lookup进行列举表出，即把表函数自变量，因变量最大值、最小值及一些自变量与因变量对应的点值列出。如上例描述的表函数可以在方程输入框写成：
[（0,0）-（10,10）]（0,0.5）（1,1）（1.5,2）（2,5）（2.5,10）其中[ ]中前面( )中0,0分另为自变量、因变量最小值，若自变量小于最小值，因变量取最小值，后面（ ）中10，10分别为自变量、因变量最大值，若自变量超出最大值，因变量取最大值，[ ]后面五个（ ）是已知自变量和因变量对值点，若自变量值不在给出点中，则自动用线性插值法求因变量对应值。自变

量、因变量的最小值、最大值可依据实际背景来确定，列出的对应值点作为已知点可从历史数据中计算或分析给出。Vensim PLE专用软件对表函数的增减性、取值间隔均匀性没有严格要求，但使用者可根据实际问题给出取值间隔、分段满足增减性的表函数。在Vensim PLE中建立的入树或流图内一个表函数必须有三部分完成，即一个自变量X,一个因变量Y及一个Y关于X的因子表，其因果关系为：
Y变量 方程可写为：
↗ ↖ Y变量=X因子表（X变量）
X因子表 X变量
上例表函数中
X因子表=[（0,0）-（10,10）]（0,0.5）（1,1）（1.5,2）（2,5）（2.5,10）
表函数的建立方法将在§8.2.2介绍。
2. Micro DYNAMO及PD PLUS中表函数表示形式与Vensim PLE软件不同的是在Micro DYNAMO及PD PLUS中有特定不同类型，其表示含义可由定义给出并固定下来。现使用 Vensim 软件的读者，可以不阅读下面内容。
⑴ Micro DYNAMO两类表函数
定义1：若TABLE（TY，X.K，XLOW，XHIGH，XINCR）
中： TY---表量名（因变量已给值）
X---自变量
XLOW---自变量X的最小值
XHIGH---自变量X的最大值
XINCR---自变量X的取值间隔
自变量取值为XLOW至XHIGH间以等间隔XINCR取X1，X2，…… Xm m个值，且
m =（XHIGH-XLOW）/XINCR + 1
对应于X1，X2，…… Xm 的TY的值在DYNAMO方程中以T方程：
T TY= E1/E2/……/E m 给出。
当X0∈（XLOW，XHIGH），但X0≠Xi（i=1,2,…,m）时，其变量值按线性插法给出，当X的值超出[XLOW，XHIGH]范围时，因变量取对应的端点值，并给出警告信息。
则 TABLE(TY，X.K, XLOW, XHIGH, XINCR）称为第一类表函数。
T TY= E1/E2/……/E m称为其表量语句，又称为T语言。

例1：已知两变量X和Y，因变量Y随自变量X变化的关系的曲线所示(图8.2.2)
自变量X从X=-3开始，按等距离取7个点得表8.2.1。设Y为辅助变量，用第一类表函数语句表示的DYNAMO语句为：
A Y.K=TABLE（TY,X.K.-3.3.1）
T TY =-20/0/10/16/20/24/30。

X
-3 -2 -1 0 1 2 3
Y
-20 0 10 16 20 24 30
表8.2.1

注1：该例在Vensim PLE中变量关系图为：

Y变量 方程可写为：
↗ ↖ Y变量=X因子表（X变量）
X因子表 X变量

X因子表 = [（-3,-20）-（3,20）](-3,-20),(-2,0),(-1,10),(0,16),(1,20),(2,24),(3,30) （曲线如图8.2.3）

当X值超出[-3,3]时，Y取对应的端点值，不给出警告信息。
定义2：若TABHL（TY，X.K，XLOW，XHIGH，XINCR）中随X的取值范围超出[XLOW,XHIGH]时,因变量取对应的端点值,但不给出错误信息外,其它内容与第一类表函数相同,则TABHL(TY，X.K，XLOW，XHIGH，XINCR）称为第二类表函数。
注2：在PD PLUS语言中，表函数中T语句中斜线改为了逗号，且定义了下面两类表函数，不妨称为第三、第四类表函数。
第三类表函数为：TABXL（TY，X.K，XLOW，XHIGH，XINCR），此表函数，除X取值范围超出[XLOW，XHIGH]范围时，因变量取端点的趋势外推值外，其它内含与第一类表函数相同。
第四类表函数为：TABPL（TY，X.K，XLOW，XHIGH，XINCR），此表函数除利用多项式使曲线在各点起滑连接代替线段连接外，其它含义与第 一类表函数相同，但表量语的数值后要加上m个零，如，对前面第一类表函数的例子，改为TABPL，则写为：
A Y.K=TABPL（TY,X.K.-3.3.1）
T TY =-20，0，10，16，20，24，30，0，0，0，0，0，0，0。

§8.2.2 表函数建立方法介绍

前面介绍的Vensim PLE专用软件，DYNAMO语言表达表函数的方式，是比较简单、容易接受的，但一个表函数的建立却不是如此容易的，往往是一个定性与定量相结合反复分析的结果，要建立一个具体表函数，肯定必须考虑所涉及的自变量、因变量的实际背景，再仔细研究其包含的一般数学问题及一般统计问题，进行深层次的量化分析，最后得出能反应变量间一般关系规律的量表作为表函数才能用于SD模型，日常生活中，可以经常碰到时间间隔相同的统计年报表、季报表、月报表，这些都是表函数。但是，在建立一个实际系统的SD模型时，这些统计报表很难作为一个完整的表函数直接放入模型中。其一，表函数不一定以时间为自变量，实践表明，表函数的自变量很多是同模型中其它一个或若干个变量的因变量；其二，统计报表没有未来若干年的预测数据，而SD模型的目标重点在对未来规律的仿真。根据以上分析和以往的经验，建立一个具体表函数，必须涉及下面基本步骤：
一、 确定变量变化范围及取值间隔：
首先根据实际背景，初步统计或估计因变量，自变量变化范围，即最小值、最大值，再依据获取数据的难易程度及灵敏度、精确度的要求来确定变量点间的取值间隔，由于Vensim PLE中自变量间隔不一定要求均匀，在定取值间隔时可依实际背景把变量范围看作一个阶段或若干阶段来定，在不同阶段，取值间隔可不一样，目的是必须准确反应变量间变化规律。
二、 确定函数的变化趋势
根据变量间因果关系极性来确定函数的增减规律，对整个变量范围要进行分析，可能某阶段呈递增态，另一个阶段呈现递减态，有的阶段不明显，根据变化幅度大小可以调整各阶段取值间隔点间隔及点密度。
三、 找出特殊点与特殊线
针对一个实际系统，建立SD模型，调用的表函数往往涉及到一些特殊取值点，如极值点、参照点、临界点，若有这些起点，在用Vensim PLE时，最好能直接给出在表函数的表达形式中，比如一个表函数的因变量Y是模型中其它变量乘积因子，最好找出Y=1对应的自变量X，把（X，1）直接放入表函数表达式中；如Y是模型中其它变量的和式，最好找出Y=0对应的自变量X，把（X，0）直接放入表函数表达式中；如Y=Y0 是表函数单调性改变的点，最好找出相应自变量X0，把（X0，Y0）直接放入表函数表达式中，等等。虽然表函数一般情况下是非线性的，但不排除在某阶段呈线性，若在某阶段呈线性，可在表函数表达式直接给出此阶段的两端点，（X1，Y1），（X2，Y2）。（此两点代表确定的那条直线）。特殊线是指Y=X，Y=Y0 或Yˊ= X0，的一些线段，反映了表函数的某种特征，若有，最好能直接给出。对于如何确定特殊点、特殊线，必须据实际背景，采用各种系统工程方法来定，可参考后面表函数实例。
四、 确定斜率
这里讲的斜率主要是指非特殊点的变化情况，在一个表函数中除找出特殊点外，更重要的是非特殊点的变化规律，在Vensim PLE中，表函数表达式内必须列出反应变化规律的若干非特殊点，这些点通过分析历史数据、预测数据得到，直接关系到表函数是否有效，其分析方法可根据实际背景借用各种相关理论及技术手段。
为了让读者更直观的了解表函数的建立过程，下面给出二个现有SD模型的表函数进行简单分析。
例1：
如图8.2.4为一个城市建设简化模型中剩余可建面积对其事业单位新建面积的影响因子ELBC关于已占有土地比LFO的表函数，即ELBC=f（LFO），两个变量均为无量纲，此函数无法用基本初等函数表示，但可采用定性定量相结合的方法建立表函数。

第一步：确立自变量LFO的变化范围及间距。
假设建模研究时，已建土地比是0.1,又不考虑土地余留，故LFO变化范围为[0.1,1]。对于[0.1,1]区间如何等份？根据整个模型的精度要求及日常实际，分成9个等份。间隔取0.1为宜。因若在仿真中LFO需取两位小数，也可由表函数的线性取值办法解决(因一般统计土地比只保留到两位小数)。
第二步：确定函数增减性。
当LFO∈[0.1,0.4]时，ELBC随LFO递增
当LFO∈[0.4，1]时，ELBC随LFO递减
前式表明年新增长面积随城市建筑的增多而增加。这是因为，在城市发展的早期阶段，大量土地有待开发，且已有企事业单位的建成会为更多企事业单位的建立创造了有利的条件。如砖厂促进了建筑加快；路面的铺砌使材料运输更方便、更迅速；水、气、电企事业单位建立能提供生活基本保障。另外，经实地分析，当已占土地比不超过40%时，建地挑选也有更大余地。出于上述定性分析得前式成立。
根据定性分析，与其它城市建立的历史事实横向比较，已占面积LFO超40%以后，由于各种主要类型企事业单位已基本建立，市场供给网已基本建立，好的建筑环境为数不多，建筑土地资金费相应增加，这些都制约着开发地区的年建筑面积的新增，则有后式成立。
第三步：确定特殊点。
1、 由历史统计数据，当土地面积比LFO=0.1时,年新增建筑面积比为60%,则
0. 7×ELBC=0.6,
ELBC=0.86
得点(0.1,0.86)
2、 有企业建设年新建面积方程分析，当LFO=0.2时，ELBC=1
3、 确定极大值点LFO=0.4时，ELBC的值，这是一个预测值，确定这种值，系统动力学本身未提供有效的方法，一般建模者常借用其它预测方法来帮助解决，如借用特尔菲方法，趋势外推，时间序列法，回归分析法、GM模型等方法。最简单的是专家咨询法。通过定性分析的ELBC最大值为1.1。得特殊点（0.4，1.1）。由实际情况，显然还有特殊点（1，0）。
第四步：确定斜率
也就是确定非特殊点对应的ELBC值。这些数据来自两部分，一部分是历史数据，另一部分是预测数据。据系统分析综合两部分结果，参考有关资料得到了图8.2.5的表函数。

例2：
如图8.2.5为我们建立的珠海市宏观经济SD模型中建城区绿地面积第三产业影响因子关于第三产业指数的表函数。这里第三产业指数是第三产业增加值的函数，根据珠海市的实际情况，珠海市1997年的建城区绿地面积占有率作为标准年，即当第三产业增加值为珠海市1997年实际值时，第三产业指数为1，这样得到了表函数的特殊点（1，1），根据珠海市第三产业的发展规模和速度，第三产业指数的取值范围定为0.5到10，根据珠海市过去的绿地面积占有率和城市的发展理念，对比1997年的情况，建城区绿地面积第三产业影响因子的取值范围定为0.9到2，珠海市是一个高度重视环境绿化的城市，近几年其绿地占有率逐年上升，可以预知随着第三产业的发展，绿地面积占有率还会提高，最后达到一个较稳定的数据，通过以上分析和专家咨询，综合珠海市的发展规划得到了上面图8.2.5表函数。该表函数可以随时根据珠海的发展进行调整，适当的时候可纳入珠海市的宏观调控计划模拟中。此SD模型在珠海市已使用了两年，得到了各方面的肯定，表函数的功能得到了充分发挥。

第三节 延迟函数
§8.3.1 Vensim PLE中延迟函数表示及使用方法
一、 定延迟函数的概念
义1 量变化需要经过一段时间的滞后才能得到响应，这种现象称为延迟。刻划延迟现象的函数称为延迟函数。
延迟是系统动力学中一个重要概念，因为在系统中存在大量延迟现象，例如培训的学员要经过一段时间才能发挥作用；投资要经过一段时间才能成为新的增生产能力；人得病，有潜伏期；污染物排放到江河之中，要经过扩散才能使江河发生污染等。另外，延迟函数的构造丰富了系统动力学理论。
二、 延迟函数的分类
发生的物流流线上的延迟称为物流延迟；发生在信息流线上的延迟称为 信息延迟。
根据以上概念，原则上所有的物流和信息流，在其流线上都会出现延迟，但我们在建模时，应抓住主要延迟进行设计，才能使复杂与精确性得到统一。由物流和信息流的不同，延迟函数分为物流指数延迟函数和信息延迟函数，在Vensim PLE中其函数名有固定的表示形式，分别为DELAY1、DELAY2、DELAY3和SMOOTH、SMOOTH3等。

三、 使用方法
延迟函数在Vensim PLE软件中直接给出，其函数在用到时可直接调用。图8.3.1是仓存Invent.dml模型的简化流率基本入树，其中变量接到的订单ORDRCV是关于延迟时间DEL的延迟函数变量，其方程为：
接到的订单ORDRCV=DELAY3（订单率ORDRS，延迟时间DEL）
流位方程为：
仓存INV=INTEG（接到的订单ORDRCV-货运率SHIP，期望的仓存DSINV）
在此该延迟变量可以不出现在入树或流图中，而可以直接放在流位方程中，这样在入树或流图中保证了流率对流位的直接作用对应关系，此时流位方程可直接写为：
仓存INV=INTEG（DELAY3（订单率ORDRS，延迟时间DEL）-货运率SHIP，期望的仓存DSINV），其运行结果和前面一样。下面是仓存invent.dml模型的所有方程，其中TEST是测试变量由测试函数组成，上机对TEST1、TEST2、TEST3、TEST4中的某个赋值1，其余的值仍为零，观察运行结果可以通过测试函数了解该模型变量的实增、稳态的下降、振动和随机扰动，这样能帮助我们弄清楚模型的反馈结构及其动态行为之间的联系（有兴趣的读者可以上机试试）。
(01) 标准货运NSHIP=
100
Units: 货运单位/周

(02) 测试输入量TEST=
STEP(10, 2 )*TEST1+TEST2*RAMP(20, 2 , 20 )+TEST3*PULSE(2, 200)+TEST4*RANDOM UNIFORM
(-5, 5 , 0.25)
Units: **undefined**

(03) 仓存调整时间IAT=
2
Units: 周

(04) 仓存调整INVADJ=
(期望的仓存DSINV-仓存INV)/(仓存调整时间IAT)
Units: 货运单位/周

(05) 仓存INV= INTEG (
接到的订单ORDRCV-货运率SHIP,
期望的仓存DSINV)
Units: 货运单位

(06) 订单率ORDRS=
平均货运率AVSHIP+仓存调整INVADJ
Units: 货运单位/周

(07) 货运率SHIP=
标准货运NSHIP+测试输入量TEST
Units: 货运单位/周

(08) 接到的订单ORDRCV=
DELAY3(订单率ORDRS, 延迟时间DEL )
Units: 货运单位/周

(09) 平均货运率AVSHIP=
SMOOTH(货运率SHIP, 信息延迟时间TAS )
Units: 货运单位/周

(10) 期望的仓存DSINV=
3*标准货运NSHIP
Units: 货运单位

(11) 信息延迟时间TAS=
2
Units: 周

(12) FINAL TIME = 25
Units: Week
The final time for the simulation.

(13) INITIAL TIME = 0
Units: Week
The initial time for the simulation.

(14) 延迟时间DEL=
3
Units: 周

(15) SAVEPER = 0.5
Units: Week
The frequency with which output is stored.

(16) TEST1=
0
Units: **undefined**

(17) TEST2=
0
Units: **undefined**

(18) TEST3=
0
Units: **undefined**

(19) TEST4=
0
Units: **undefined**

(20) TIME STEP = 0.25
Units: Week
The time step for the simulation.

为了更好地理解延迟函数的内在涵义，本书采用Micro PYNAMO语言程序来刻画。从§8.3.2节开始阐述各种延迟函数的内部结构原理及有关理论分析，这部分内容对于需深入掌握延迟函数的读者必须阅读。
13