十 系统动力学理论基础
（一）系统动力学概述
系统动力学（System Dynamics，简称SD）由美国麻省理工学院福瑞斯特（Jay W.Forrester）教授创始于1956年，用来分析研究复杂信息反馈系统。系统动力学初期主要应用于工业企业管理的生产与雇员情况波动、市场股票与市场增长的不稳定性等问题。1961年，福瑞斯特出版了工业动力学学科的第一本专著《工业动力学》，为系统动力学的发展奠定了基础。系统动力学经过不断应用和发展已逐渐成熟，其理论与应用研究涉及管理、社会等各种学科和领域，认为系统的内部动态结构和反馈机制决定了系统的行为模式和特性。[71]系统动力学采用定性与定量分析相结合的方式，从系统内部机制与微观结构着手，通过剖析系统进行建模分析，利用计算机模拟技术探索系统内部结构及其动态行为关系，并找寻解决系统内在问题的相关对策建议。因此，系统动力学模型特别适合分析结合社会、经济等非线性的复杂大系统问题。
系统动力学与工业工程学、运筹学、计量经济学等方法相比，有以下几个突出的特点：擅长处理长周期和长期性问题，适合进行数据缺少条件下的研究，擅长处理高阶、非线性、时变问题，常被用来进行情景分析。
（二）系统建模原理
1.一阶系统的结构行为
系统动力学在时域中，采用状态空间法描述系统结构，进而分析和研究系统。系统向量形式的状态方程为：
X=f（X，U，t），X∈Rm，U∈Rr
其中，R为欧式空间；向量X为m维，U为r维。
对于一阶定常自由系统，m=0，r=0时不变系统，则有：
dx/dt=f（x）=a0+a1x+a2x3+…
如果等式右端保留一次项和常数项，上式变为：
dx/dt=f（x）=a0+a1x
其原函数为：
其中，c为常数。
可见，当a1＞0时，系统呈指数增长特征，自然界中的动植物在无约束的条件下表现为此种指数特性。当a1＜0时，系统呈指数减少，自然界中的放射性射线强度的衰减、被照杀或药杀的细菌消亡过程、人口死亡速率大于出生速率的过程都表现为指数减少。
如果等式右端保留二次项、一次项和常数项，dx/dt=f（x）=a0+a1x+a2x3+…则变为：
dx/dt=f（x）=a0+a1x+a2x2
其原函数为：
其中，c为常数。
可见，a1＞0，a2＞0，t→∞，x趋于定值，此时系统呈S形增长特性。
综上所述，对于一阶系统，无论控制作用多么复杂，系统或者呈指数增加，或者呈指数衰减，或者呈渐进增长。一旦趋于某个既定的目标，则出现平衡且会永远保持下去。因此一阶系统不会发生超调，更不会发生振荡。
2.二阶系统的结构行为
二阶系统比一阶系统更为复杂，一般在一个系统中包含两个独立的状态变量，并且这两个状态变量在同一个回路中。系统动力学用状态空间法在时域中表述系统的结构和研究系统的功能与行为。系统向量形式的状态方程如下：
X=f（X，U，t），X∈Rm，U∈Rr，其中，R为欧氏空间。
为了简化叙述，这里以二阶定常自由系统为例说明二阶系统的描述问题，这时系统为m=2，r=0时不变系统，其向量方程可表示为：
其中，
表示dx/dt组成的列向量；A为状态转移矩阵。这个向量方程可以从数学角度进行分析，研究其行为特性。下面以库存系统为例，对二阶系统的行为特性进行分析（见图1-13）。
图1-13 二阶库存系统存量流量
在上述库存系统中，如果有在途库存，系统就从一阶系统变为二阶系统。之所以称其为二阶系统，是因为在系统的最大回路中包含了两个水平变量TL和IL。其中，IL为实际库存量，TL为在途库存量，LT为平均订货提前期，Constant为常量，其他变量与一阶库存系统相同。
这里为了便于分析，假设Constant=0，该系统的数学描述为：
将上式对t求导，整理可得：
可见上式为二阶线性非齐次微分方程，其特征方程为：
该特征方程有两个根r1和r2的取值，可得上述齐次方程的通解如下：
（1）r1≠r2时，
（2）r1=r2=r时，IL=（c1+c2t）en
（3）r=α+βi时，IL=［c1cos（βt）+c2sin（βt）］eat
从齐次方程的通解来看，二阶系统的行为要比一阶系统复杂得多，不但可能出现超调，也可能出现振荡现象。
系统的行为特性同特征方程的变化有关，即同库存偏差调节时间与订货提前期的相对取值有关。即
（1）AT＞4IT，特征方程有两个不等实根r1和r2，r1≠r2，库存的变化规律如式（1）所示；
（2）AT=4IT，特征方程有两个相等的实根r1=r2=r，库存的变化规律如式（2）所示；
（3）AT＜4IT，特征方程有两个虚根r=α+βi，库存的变化规律如式（3）所示。
在此只考虑了一个简单的线性系统，而且也只分析了其对应的奇次方程解的情况。可见，在二阶系统中，由于库存调节时间AT和订货提前期LT的比例关系不同，系统的行为特征也不同，如何控制订货使库存达到所期望的目标，与一阶系统相比要复杂得多。
（三）系统动力学流率基本入树建模法
南昌大学贾仁安教授及其研究小组于1998年创立了流率基本入树建模法[72]，该方法以还原论思想为指导，将图论中生成树理论应用于动态复杂系统的反馈结构分析。此方法将所研究的整个系统按研究目的划分为不同的子系统，然后设定每个子系统内部的流位、流率和辅助变量，抓住系统反馈结构变量中最基本的流率变量，用一组以流率变量为根的树模型来刻画系统内各变量之间的因果关系，最后通过引入嵌运算构建系统网络存量流量图。流率基本入树建模法有两大好处：
（1）有利于分部分、分子系统进行规范化建模，提高线段性思考的集中度与精确度；有利于用整体论与还原论相结合的思想方法对问题进行有效研究；有利于仿真方程的建立。
（2）为利用代数的方法研究动态反馈复杂性系统问题提供了可能性。有了流率基本入树模型，通过将入树的枝转化为枝向量，构造枝向量行列式和枝向量矩阵，就可以利用代数方法进行系统的动态反馈复杂性分析，从而实现图论与线性代数在研究系统反馈动态复杂性问题中的结合。
1.流率基本入树建模法的基本概念
定义2.1 若t∈T，一个动态有向图T（t）=（V（t），X（t））中，存在一个点v（t）∈V（t），使T（t）中的任何一点u（t）∈V（t），有且仅有一条由u（t）至v（t）的有向道路，则此有向图T（t）称为一棵入树，且v（t）称为树根，满足入度d-u（t）=0的u（t）称为树尾，从树根至树尾的一条有向道路称为一根树枝。
定义2.2 在系统动力学存量流量图中，以流率为树根、以流位为树尾的入树T（t）称为流率入树。流率入树T（t）中含流位的个数称为入树的阶数，从树尾沿一枝至树根所含流位的个数称为这枝的枝阶长度。流率入树最大枝阶长度称为该入树的阶长度。
定义2.3 各枝阶长度为1的流率入树称为流率基本入树。
定义2.4 不真包含在任何其他流率基本入树的流率基本入树称为极大流率基本入树。
定义2.5 存量流量图中任何一个子图称为半子存量图，满足含流位L（t）有其流率R（t）（或流出率R1（t），或流入率R2（t））的半子存时流量图称为子存量流量图。
定义2.6 已知t∈T，半子存量流量图G1（t）=（Q1（t），E1（t），F1（t）），G2（t）=（Q2（t），E2（t），F2（t））则
（1）作G1（t）∪G2（t）且保持F1（t）和F2（t）确定的映射关系。
（2）若流率RP（t）及其对应的流位Lp（t）在Gi（t）（i=1，2）中，则在（1）的基础上再增加一条弧，构成因果链：RP（t）→Lp（t），同时给出实际意义下的因果链极性。
由（1）和（2）所得到一个新的半子存量流量图G（t），定义这种运算为嵌运算，嵌运算记为
，则
。
嵌运算满足以下性质：
交换律：
结合律：
2.流率基本入树建模法的建模步骤
在引入上述基本概念的基础上，给出流率基本入树建模法的基本步骤。
步骤1：通过系统分析，建立流位流率系：
{（L1（t），R1（t））（L2（t），R2（t）），…，（Ln（t），Rn（t））}。
步骤2：分别建立以流率变量Ri（t）为根、以流位变量Lj（t）为尾的，且流位变量直接或通过辅助变量控制流率变量的流率基本入树，可得图1-14所示的流率基本入树模型。
图1-14 流率基本入树模型
图1-14中Aij（t），Bij（t）（其中i，j=1，2，…，n）可能是多个辅助变量构成的有向链。
步骤3 对这些基本入树模型T1（t），T2（t），…，Tn（t）作嵌运算，即顶点与顶点并、弧与弧并、流率与对应流位相连，则可得流位流率系下的存量流量图模型。
建立流率基本入树和直接建立存量流量图模型是建立系统结构模型的两个等价方法。
同一流位流率系下的网络存量流量图G（t）与流率基本入树模型T1（t），T2（t），…，Tn（t）具有当且仅当关系。即在同一流位流率系下，由网络存量流量图G（t）分解可得入树模型T1（t），T2（t），…，Tn（t），由入树模型T1（t），T2（t），…，Tn（t）可得网络存量流量图G（t）。