江苏大学系统动力学课 第五讲.ppt

Vensim建模基础,江苏大学工商管理学院李文元,Contents,概述,变量与方程,延迟、平滑和平滑函数,函数,1,2,3,4,准确度与运行时间单位的选择,5,5.1概述,DYNAMO，取名来自 Dynamic Models(动态模型)的混合缩写。顾名思义，DYNAMO命名的涵意在于建立真实系统的模型，藉助计算机进行系统结构、功能与动态行为的模拟。 DYNAMO和系统动力学的关系，可追溯到50年代系统动力学发展的初期。DYNAMO的前身称SIMPLE (Simulation of Industrial Management Problems with Lots of Equation),以库存系统为例，为简单起见，考虑输入、输出速率为常数的情况。 假定每月发货与入库各为100与80件，则库存INV每月减少20件，其动态行为是线性的，以图形表示就是随时间变化的直线。 可用数学式表达： INV现在INV过去十(时间间隔)x(纯速率) 若库存量在5个月前为l200件，则： INV现在1200件+(5月)*(80件/月-100件/月) 1200+5*(-20) 1200-100 1100(件),5.1概述,5.1概述,2019/11/1,5.1概述,2019/11/1,5.1概述,当速率随时间变化时，可以把连续的时间分割成小的时间间隔，并假定在各小间隔内速率是固定的，然后，藉助计算机逐段地一一加以计算。若计算的时间间隔足够小，速率变动不大，则此结果将与从微分方程获得的精确解(如果可能求得的话)十分接近。,2019/11/1,5.1概述,在DYNAMO中，给变量带上时间下标以区别在时间上的先后。 英文字母K表示现在，J表示刚刚过去的那一时刻，L表示紧随当前的未来的那一时刻。DT表示J与K或K与L之间的时间长度。,5.1概述,库存方程可用DYNAMO表示如下： INV.K = INV.J + DT*(ORRE.JK - SH.JK) 式中： INV.K 库存现有量； INV.J DT前的库存量； DT 计算的时间间隔； ORRE 在JK间隔内收到的订货量 SH 在JK间隔内的发货量。,4.2变量与方程,5.2.1变量总类 在Vensim中有如下几种变量 Level variable:状态变量，是最终决定系统行为的变量，随着时间变化，当前时刻的值等于过去时刻的值加上这一段时间的变化量。 Rate vaviable:速率变量，是直接改变状态变量值的变量，反映状态变量输入或输出的速度，本质上和辅助变量没有区别。 Auxiliary variable:辅助变量，辅助变量值有系统中其他变量计算获得，当前时刻的值和历史的值都是相互独立的,2019/11/1,5.2变量与方程,Constant variable:常量，常量值不随时间变化 Exogenous variable：外生变量，随时间变化，但是这种变化不是系统中其他变量引起的。,2019/11/1,5.2变量与方程,5.2.2变量方程 （1）状态变量,2019/11/1,LvS(t)是t 时刻状态变量值 rateS(t)是该状态变量变化的速率,5.2变量与方程,如果写成离散的方程形式 LEVEL.K=LEVEL.J+DT*(INFLOW.JK-OUTFLOW.JK) 式中： LEVEL 水平(状态)变量； INFLOW 输入速率(变化率)； OUTFLOW 输出速率(变化率)； DT 计算间隔(从J时刻到K时刻)。,2019/11/1,5.2变量与方程,2019/11/1,5.2变量与方程,2019/11/1,式中：Pop为人口数量（人）；births为出生率（人/年）；deaths为死亡率（人/年）；immigration为迁入率（人/年）；emigration为迁出率（人/年）,5.2变量与方程,（2）速率变量,2019/11/1,rateS（t)是状态变量变化的速率；LvS（t)是t时刻状态变量值；aux（t),exo(t)分别是t时刻辅助变量和外生变量；const是常数,4.2变量与方程,Births和deaths表示每年出生和死亡的人口数量，有births=Pop（t）BR ；deaths=Pop（t）/AL 式中： Pop为人口数量（人）；births为出生率（人/年）；deaths为死亡率（人/年）；BR为出生率系数（1/年）AL为平均寿命（年）,2019/11/1,5.2变量与方程,（3）辅助方程 辅助方程是在反馈系统中描述信息的运算方程,2019/11/1,5.2变量与方程,2019/11/1,5.2变量与方程,5.2.3表函数 模型中往往需要用辅助变量描述某些变量间的非线性关系，显然简单地由其它变量进行代数组合的辅助变量已不能胜任。比较方便是能够以图形方式给出这种非线性关系。Vensim建模中可以通过表函数（lookup function）来实现这一功能。,2019/11/1,5.2变量与方程,建立表函数的步骤 确定表函数中自变量和因变量 明确变量的取值变化范围，通常通过历史数据和历史数据来确定 考虑曲线的形状和斜率，在什么范围内事平缓的，什么范围内曲线坡度高，要符合自变量和因变量之间影响关系，其中正斜率表示正反馈，负斜率表示反馈 选择合适的曲线端点、驻点、拐点等。 做出图形，通常X轴表示自变量，Y轴表示因变量,2019/11/1,5.2变量与方程,Vensim中表函数的数学描述形式是 lookup name,2019/11/1,5.2变量与方程,考虑一个企业建设的例子。企业建设BC是变化率，方程如下： BC.KL = NCF*BS.K*ELBC.K 式中： BC 企业建设(个/年)； NCF 额定建设系数(1/年)； BS 企业(个)； ELBC 土地对企业建设的影响(无量纲)。,2019/11/1,5.2变量与方程,讨论辅助变量ELBC方程的建立问题。 假定，ELBC为土地占用系数LF0的函数，LFO变化范围为从零至1。当LPO=0时，表示土地未被占用，LFO=1则表示土地被全部占用。 ELBC与LFO的非线性关系,可用图表表示。,5.2变量与方程,LPO=0，土地未被占用 LFO=1，土地被全部占用,5.3延迟、平滑和平滑函数,系统动力学模型中包含的物质流与信息流可能存在延迟。 在现实生活中，延迟的现象比比皆是： 厂家向顾客交付订货； 感染潜伏期病症； 播种与收获； 投资与回报； 。 延迟是信息反馈系统结构中颇为重要的一个角色。,2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,5.3延迟、平滑和平滑函数,加入的新状态变量为处于潜伏期人口INC，其输入速率是感染率INF，输出速率是疾病症候显现率SYMP,有方程组如下： INC(t0 )= TSS*INF(t0 ) SYMP(t) = INC(t) /TSS 式中TSS为潜伏期，比如流感的潜伏期为3天。,2019/11/1,5.3延迟、平滑和平滑函数,上述方程式可用Vensim中的DELAYl函数代替，功能相同，但简明方便得多(可以取消INC变量，即INC为隐含水平变量) SYMP（t） = DELAYl(INF（t）,TSS),2019/11/1,5.3延迟、平滑和平滑函数,(1)DELAY1：一阶延迟 DELAY1是指一阶物质延迟。它隐含了一个状态变量（INC）。以DELAY1替代一组方程，使用方便，缺点是该状态变量被隐含，不能输出的它的结果。,2019/11/1,DELAY1(INF,TSS),5.3延迟、平滑和平滑函数,(2) DELAY1I：一阶延迟 DELAY1I形式是DELAY1I( in , dtime , init )，含义是DELAY1I( input, delay time, initial value) DELAY1I=LV(t)/DT,2019/11/1,式中：IV和DT分别表示initial value和delay time,5.3延迟、平滑和平滑函数,如果式中IV等于input初始值，则DELAY1I就简化为DELAY1，形式为DELAY1(input,delay time) 此时，系统自动根据上面的方程组计算出LV的初始值，使输入input和延迟输出DELAY1相等。,2019/11/1,5.3延迟、平滑和平滑函数,(3）三阶延迟,2019/11/1,DELAY3(INF,TSS),5.3延迟、平滑和平滑函数,DELAY3I形式是DELAY1I( in , dtime , init )，含义是DELAY3I( input, delay time, initial value) 反映的方程组为: DL=delay time/3,2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,5.3延迟、平滑和平滑函数,如果式中IV等于input初始值，则DELAY3I就简化为DELAY3，形式为DELAY3(input,delay time),2019/11/1,5.3延迟、平滑和平滑函数,(4)DELAY FIXED DELAY FIXED函数是DELAY FIXED( in , dtime , init )。这个函数的作用是将input值延迟一段时间后等值输出，在延迟这一段时间内输出initial value.这里的delay time 可以是一个表达式，但必须满足在0时刻能计算出初值。,2019/11/1,5.3延迟、平滑和平滑函数,对于DELAY FIXED函数只可以在其内部建立表达式，但是其本身不可以作为表达式的一部分，不可以与其他表达式共同构成更大的表达式。 D=DELAY FIXED(B,T,B)+1：错误 DM=DELAY FIXED(MAX(A,B),C,A）正确,2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,平均延迟时间average delay time=30,5.3延迟、平滑和平滑函数,2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,物质延迟的阶次,阶次的含义在延迟结构中，指的是延迟环节内部包含的水平变量数。 下图给出了不同阶次延迟环节的输入量发生突变时，(即为阶跃输入时)，其相应的输出量的曲线。这是一簇曲线，包括1，2，3，6和12阶延迟的响应。 曲线簇表明，l阶与3阶的延迟特性彼此差别很大。1阶延迟表现出简单的指数形增长的特性，2阶延迟开始表现出S形增长特性，3阶时的S形增长特性已较明显，6阶与12阶的S形增长特性也就更加突出了。 随着阶数的增加，延迟环节的响应的增长模式本质相同，其错开程度取决于延迟时间。显而易见，l阶与3阶的曲线差别很大，增长模式全然不同，但3阶曲线与6阶甚至12阶曲线相比则无本质差别，只是程度上差异而已，同样是S形模式。,因此，vensim中仅备有DELAYl与DELAY3函数，而无更高阶次的单个函数。,5.3延迟、平滑和平滑函数,5.3.2平滑函数与信息平滑 与物质在系统中流动存在延迟类似，信息在系统中传递也存在延迟。 商品信息的传递一般都带有延迟的特性。这种信息传递的延迟，往往是系统结构中不可避免的组成部分。 平均或平滑信息导致延迟。系统动力学中描述信息的延迟可用平滑函数。,2019/11/1,5.3延迟、平滑和平滑函数,以生产经营管理中的现象为例，企业领导人决不会认定某日销售额突增的信息作为长远的趋势，把它作为库存、生产安排与招工等问题决策的依据。 决策者总是力图从销售信息中排除随机的因素，找出真实的趋势。换言之，对销售信息应求其在一段期间内的平均值。这种“平均”与“平滑”的处理方式在系统动力学的模型中屡见不鲜。 Vensim提供这一功能的函数称为平滑。 信息的平滑或平均实质上是一种积累过程。它们可以包含一个或多个水平变量。,2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,5.3延迟、平滑和平滑函数,在Vensim中，用SMOOTH函数可以表示上述方程组,2019/11/1,5.3延迟、平滑和平滑函数,被平滑的变量可以是水平、速率和辅助变量。平滑时间STIME通常为常数，但也可以是变量。 直观地看，平滑时间系指变量VAR经积累达到指数加权滑动平均值所需的时间。,2019/11/1,平滑函数对输入量的响应持性,若变量VAR为阶跃函数(突增后保持恒定)，其平滑值SVAR将渐渐趋于此恒定值。 若VAR是一个脉冲，SVAR不能达到VAR的幅值，并按另一指数式的寻的特性下降。 若VAR是一振荡的输入量，其平滑值SVAR亦将随着振荡，但幅值要小得多。 因此平滑函数具有平滑原变量的激烈起伏功能。经平滑得到的平均值正是所期望的真实趋势。,5.3延迟、平滑和平滑函数,5.3.3信息延迟 （1）一阶信息延迟函数SMOOTHI SMOOTHI的函数形式是SMOOTHI( in , stime , inival )含义是SMOOTHI( input, delay time, initial value),2019/11/1,5.3延迟、平滑和平滑函数,方程组是,2019/11/1,5.3延迟、平滑和平滑函数,(2)三阶信息延迟函数SMOOTH3I SMOOTHI的函数形式是SMOOTH3I( in , stime , inival )含义是SMOOTH3I( input, delay time, initial value),2019/11/1,5.3延迟、平滑和平滑函数,2019/11/1,5.4函数,5.4.1数学函数,2019/11/1,5.4函数,2019/11/1,5.4函数,2019/11/1,5.4函数,5.4.2逻辑函数 （1）MAX(A,B)取A，B中较大者 (2) MIN(A,B)取A，B中较小者,2019/11/1,5.4函数,(3)条件函数 IF THEN ELSE( cond , ontrue , onfalse )=,2019/11/1,5.4函数,5.4.3测试函数,2019/11/1,5.4函数,方程 发货率SHIP=正常发货率NSHIP+测试函数输入值TEST 发货率平滑时间TAS=2 平均发货率AVSHIP=SMOOTH(发货率SHIP, 发货率平滑时间TAS ) 库存INV= INTEG (订货率ORDERS-发货率SHIP, 300) 库存调节时间LAT=2,2019/11/1,5.4函数,库存调节率INVADJ=(期望库存DSINV-库存INV)/库存调节时间LAT 延迟时间DEL=3 期望库存DSINV=期望库存覆盖时间DIC*正常发货率NSHIP 期望库存覆盖时间DIC=3 正常发货率NSHIP=100 测试函数输入值TEST=STEP(10, 2) 订货率ORDERS=DELAY3(平均发货率AVSHIP+库存调节率INVADJ, 延迟时间DEL ),2019/11/1,5.4函数,（1）阶跃函数 阶跃函数是幅值在特定时刻发生突变的函数，常被用来突然改变变量数值 STEP( height , stime ) Height表示阶跃的幅度，stime表示阶跃发生的时刻。在时刻steptime前STEP函数值是0；当时间等于或大于steptime时，STEP的值等于height,2019/11/1,5.4函数,2019/11/1,发货率SHIP的正常值NSHIP为100件/周，阶跃函数于第2周(STRT=2)，由0突增加10件/周，使SHIP增加了10，达110件/周并保持不变。 虚线为发货率SHIP随时间变化曲线；实线为库存INV相应的变化特性。 当发货率突增10后，库存量从期望值DSINV300件，迅速下降，与此同时，订货率ORDRS增加，以补充库存的减少。然而从订货到厂家交货存在延迟(模型中设定为平均3周)，所以库存量继续下降，大约到第6周达谷底；之后才回升。 由于整个系统存在延迟环节与负反馈的调节作用，库存量呈现为围绕其期望值衰减振荡特性。,5.4函数,(2)斜坡函数 斜坡函数是一种连续增长或下降的时间的线性函数，其形式是RAMP( slope , start , finish ) 式中：slope表示线性函数的斜率，start表示斜坡函数的起始时刻，finish表示函数结束的时刻，在star时刻前，ramp函数值取0，在start时刻到finish时刻之间，RAMP=slop*(time-start);在finish时刻后，RAMP=slop*(finish-start),2019/11/1,5.4函数,Test= RAMP( 20, 2, 25) RAMP函数代表对发货率线性增长的需求。INV几乎跌落至0，但对发货率无任何影响。 可以预料，若斜坡函数的斜率值再大些，INV将出现负值。显然这是不合理的。在现实生活中，库存INV的下降情况需对发货率SHIP施加负反馈的信息，这一点在此简化库存模型中未予考虑,2019/11/1,5.4函数,（3）脉冲函数 脉冲函数提供瞬时冲击的方法，在现实生活中描述变量的独立变化现象，也就是变量在每一次变动后回复原值。两种脉冲函数： PULSE( start , duration ) PULSE TRAIN( start , duration , repeattime , end ),2019/11/1,5.4函数,PULSE( start , duration )是一个单脉冲函数，函数值在（start，start+duration）时间内等于1，其他时刻为0. PULSE(2 , 0.25)*10,2019/11/1,5.4函数,PULSE(2 , 0.5)*10,2019/11/1,5.4函数,2019/11/1,5.4函数,PULSE TRAIN( start , duration , repeattime , end )是一个周期性脉冲函数，脉冲从时刻start开始，函数值为1，延续时间为duration，然后再start+ repeattime时刻重复一次前面的脉冲，这样周期性脉冲直到end时刻结束，相邻两个脉冲的时间间隔为repeattime，在脉冲时间之外函数值均为0，需要注意的是，如果repeattime时间小于duration，则脉冲会从start一直持续到end。如果duration小于时间步长，则按1个步长计算 Test=PULSE TRAIN(2,0.25 ,5 , 25)*10,2019/11/1,5.4函数,2019/11/1,5.4函数,（4）正弦函数 正弦函数可用来测试模型对于振荡(比如正弦)输入变量的响应。 正弦函数的DYNAMO表达式为： A*SIN(6.283*TIME/B) 式中： A 振荡幅度； B 振荡周期； 6.283 2pi(pi为圆周率)。,2019/11/1,5.4函数,2019/11/1,TEST = AMP*SIN(6.283*TIME/PER) AMP = 10 PER = 5(左图),5.4函数,2019/11/1,TEST = AMP*SIN(6.283*TIME/PER) AMP = 1