枝向量矩阵反馈环算法的MATLAB实现

1、枝向量矩阵反馈环算法的MATLAB实现系统动力学是通过建立流位､流率系来研究信息反馈系统的一门科学｡用系统动力学来处理复杂系统时一个关键的问题就是反馈环的计算｡一个社会经济系统往往是由许多反馈环构成的复杂系统,如美国的福瑞斯特等人建立的世界模型n就是由82条反馈环组成的｡因此,要对此系统的动态性､复杂性进行分析,就必须确定其反馈环的集合,这对于其结构分析､基模确立､模型调试､结果分析都非常重要｡针对

2、反馈环的计算也有不少方法,如常见的反馈环计算法有反馈环图示计算法､枝向量行列式反馈环计算法及枝向量矩阵反馈环计算法等1｡这些算法都还未编程,在计算机上难以实现｡其困难一方面在于算法本身,有的方法根本就没办法用计算机实现;另一方面是算法程序实现上的困难,像枝向量行列式反馈环计算法的行列式计算法则不同于数值为元素的行列式的计算｡对于这些问题,本文比较了几种常见的反馈环计算方法,并给出了运用枝向量矩阵求反馈环算法的MATLA改现｡算例表明,所设计的程序具有很好的应用效果&#65377

3、;1 反馈环计算的常见方法比较常见的反馈环计算方法有反馈环图示计算法､枝向量行列式反馈环计算法及枝向量矩阵反馈环计算法1｡图示计算法是一种通过计算入树模型的所有反馈环的图示法｡此方法思路清晰,但当流图非常复杂时,用此图示法就感到烦琐,而且无法用计算机实现,所以效果不佳｡下面从时间复杂度上来比较在计算出所有反馈环时的行列式算法和矩阵算法｡为了具有可比性,笔者将行列式算法中枝向量行列式元素的乘法和矩阵算法中枝向量矩阵元素的乘法分别作为两个算法的基本操作｡由行列式算法很

4、容易得到其时间复杂度O(nXn!)｡其中n是枝向量行列式的阶数｡矩阵算法能够求出当向入树模型中依次增加入树时新增的各阶反馈环,把入树模型中所有入树新增的反馈环加起来即可得到模型中所有的反馈环｡由此可以计算出矩阵算法求所有反馈环的时间复杂度为O(n(n+1)/2)2)｡其中n为枝向量矩阵的阶数｡若以矩阵算法的运算量为分子,以行列式算法的运算量为分母,求极限有可见,矩阵算法的时间复杂度小于行列式算法的的时间复杂度｡计算表明,随着n增大，特别是当n>4时，矩阵算法

5、的优势更加明显,这对于算法的计算机实现非常重要｡鉴于此给出了矩阵算法的计算机实现｡2 枝向量矩阵算法的实现凭借MATLAB虽大的数值计算和字符处理功能,实现了运用矩阵算法计算出入树模型中的所有反馈环24｡为了方便表达,用类MATLAB勺伪代码来书写相应的计算机实现算法｡枝向量矩阵有流率基本入树枝向量矩阵和强简化流率基本入树枝向量矩阵5,在矩阵算法的基础上分别给出由这两个枝向量矩阵计算所有反馈环的计算机实现算法,下面先给出一些相关的概念｡2.1 相关概念定义1以流率基本入树T1(t),T

6、2(t),Tn(t)中的枝向量为元素,依次排列构成的向量(Ri(t),Aij(t),Lj(t)称为枝向量｡其中Ri(t)为流率,Lj(t)为流位,Aij(t)为入树中枝的辅助变量依次排列的组合;如果是多枝的变量构成枝向量,Aij(t)隐含有流率与流位｡定义2对于枝向量(Ri(t),Aij(t),Lj(t),若Ri(t)就是Lj(t)所对应的流率,则此枝向量对应流图中的一个反馈环,把此枝向量称为反馈环枝向量或反馈环向量｡定义3枝向量加法(Ri(t),Aij(t),Lj(t)+(Rt(t),Atp(t),Lp(t)仅表示在流率基本入

7、树中存在(Ri(t),Aij(t),Lj(t)和(Rt(t),Atp(t),Lp(t)对应的两枝｡定义4不包含加法的非零枝向量为单元枝向量,几个单元枝向量做加法所组成的量称为这几个枝向量的复合枝向量｡以后若不特别声明,枝向量包括零向量､单元枝向量和复合枝向量｡定义5枝向量乘法(Ri(t),Aij(t),Lj(t),Rj(t),Atp(t),Lp(t)Rt(t)是 Lj(t)的流？弓是伊礁船0站恐形尴嗤？变量(Rt(t),Atp(t),Lp(t),Rp(t),Aij(t),Lj(t)Ri(t)是Lp(t)的流

8、？弓是伊礁船0站恐形尴嗤？变量0其他情况定义6由入树Ti(t)(包括流率基本入树及简化､强简化流率基本入树)的枝向量(Ri(t),Ai1(t),L1(t),(Ri(t),Ai2(t),L2(t),(Ri(t),Ain,Ln(t)(i=1,2,n)构成的方阵：称为此入树T1(t),T2(t),Tn(t)模型的枝向量矩阵｡如果矩阵中的元素在入树中不存在,把它记为零向量｡此枝向量矩阵的乘法法则与代数中给出的数矩阵乘法法则相同｡定义8将入树模型枝向量矩阵A的全体aii(i=1,2,n)置为零,其他aij(i中j)不

9、变,所得矩阵称为此入树模型的对角置零枝向量矩阵｡将枝向量矩阵A的最后一行全体元素置为0,得到的矩阵记为A｡A称为A的末行置零矩阵｡2.2 由流率基本入树枝向量矩阵计算所有反馈环的算法通过流率基本入树建模法6可以建立起模型的流率基本入树,按照上述定义就可以得到相应的枝向量矩阵｡要通过枝向量矩阵计算出模型的所有反馈环,就要实现枝向量矩阵的乘法,而这就必须首先要实现枝向量矩阵中的元素(枝向量)的加法和乘法｡由定义3可以给出两个枝向量做加法的具体算法｡算法1c=VPlus

10、(a,b)输入:枝向量a和b;输出:枝向量a和b做加法的和c｡(1) 如果a为零向量,则c=b,转(4);否则转(2)｡(2) 如果b为零向量,则c=a,转(4);否则转(3)｡(3) 把a和b构成复合枝向量c,转(4)｡(4) 结束｡由定义4设计了两个枝向量(单元枝向量或零向量)做乘法的具体算法｡算法2c=DYMultiply(a,b)输入:单元枝向量a,b;输出:做枝向量乘法的结果c｡(1) 如果a为零向量或b为零向量则输出结果c为零向量

11、,转(4);否则转(2)｡(2) 若b中的第一个流率恰是a中最后一个流位所对应的流率且a､b中无相同变量,则将a中的变量依次放在b中的变量前面构成枝向量c中的变量,转(4);否则转(3)｡(3) 若a中的第一个流率恰是b中最后一个流位所对应的流率且a､b中无相同变量,则将b中的变量依次放在a中的变量前面构成c中的变量,转(4);否则c为零向量,转(4)｡(4) 结束｡由于枝向量矩阵中的元素可能为复合枝向量,对于任意两个枝向量之间的乘法可通过下面的具体算法实现:算法

12、3c=Multiply(a,b)输入:枝向量a,b;输出:枝向量c｡(1) 如果a､b都是单元枝向量或零向量,则c=DYMultiply(a,b),转(5);否则转(2)｡(2)如果a为单元枝向量或零向量而b为复合枝向量,则将b分为两部分,记为bl和b2｡其中,b1为b中的第一个单元枝向量,b2为b中除了第一个单元枝向量之外的其他枝向量构成的枝向量｡c1=DYMultiply(a,b1),c2=Multiply(a,b2)｡c=VPlus(c1,c2),转(5)&

13、amp;#65377;否则,转(3)｡(3)若a为复合枝向量而b为单元枝向量或零向量,则参照(2)中的D作同本效t理；否则，转(4)｡(4) 此时a和b都为复合枝向量,对它们进行如下处理:将b分为两部分,记为bl和b2｡其中,b1为b中的第一个单元枝向量,b2为b中除了第一个单元枝向量之外的其他枝向量构成的枝向量｡将a用同样的方法分为两部分,分别记为a1､a2｡c1=DYMultiply(a1,b1),c2=Multiply(a1,b2);c3=Multiply(

14、a2,b1),c4=Multiply(a2,b2)｡将c1､c2､c3和c4做枝向量加法得到c,转(5)｡(5) 结束｡实现了枝向量的基本运算后,就可将它们用到枝向量矩阵的乘法运算中｡下面分两步来实现由枝向量矩阵计算出入树模型的所有反馈环,首先给出向入树模型中增加入树时新增反馈环的算法,然后把所有入树新增的反馈环汇总得到并输出所有的反馈环｡(1) 当向入树模型中增加入树时,同时也会相应地增加反馈环的数量｡下面给出当增加一

15、棵入树时新增反馈环的具体算法｡算法4result=CreatLoop(A,k)输入:枝向量矩阵A和A的左上角子枝向量矩阵的阶数k;输出:result为1xk的枝向量矩阵，用来保存计算得到的1k阶反馈环｡result=cell(1,k);%初始化result若k=1,则result1,1=A1,1,转;否则转｡result1,1=Ak,k;%对角上的元素构成一阶反馈环%乍A的左上角k阶枝向量矩阵Ak%把枝向量矩阵Ak的对角置零 %取Ak的第k行元素构成行枝向量矩阵Xk %做Ak的末行置零矩阵A0A0=Ak;%初始化A0fori=1

16、:k;A0k,i=0;end峨第1n-1阶枝向量矩阵乘法得到2k阶新增的反馈环｡fori=2:k%初始化Yk用来保存计算出的相应反馈环Yk=cell(1,k);forp=1:k;Yk1,p=0;end?r%ffiXk和A0做枝向量矩阵乘法%Xk只有最后一个元素含有反馈环结束｡(2) 将所有的入树依次添加到模型中,通过调用算法4就可以计算出所有的反馈环｡计算并输入树模型中的所有反馈环的具体算法如下:算法5AllLoop(A)输入:枝向量矩阵A;输出:枝向量矩阵A所能构成的所有的反馈环｡n=length(A);

17、loop=cell(n,n);%loop用来存放枝向量矩阵A所能构成的所有反馈环fo门=1:n?rNewLoop=CreatLoop(A,i);?r%ffiNewLoop赋给矢!阵loop中的第i行forj=1:i;loopi,j=NewLoop1,j;endend%俞出各阶反馈环?rsum=0;%sum用来统计生成的反馈环的总数?rforj=1:nnum=0;%nurmg计生成的同一阶的反馈环的总数fori=1:n若i刁，则如果loopi,j不是零向量,则输出loopi,j中包含的所有j阶反馈环,每输出一个反馈环,使num=num+1｡end若num=0,则生成的反馈环

18、中不含有j阶反馈环；否则输出所有的num条j阶反馈环，并使sum=sum+nun｡end输出显示共有sum条反馈环｡结束｡1.3 由强简化流率基本入树枝向量矩阵计算所有反馈环的算在文献1,5中提到了强简化流率基本入树模型,并给出了强简化流率基本入树的枝向量矩阵｡为了方便,下面将通过强简化流率基本入树所得的枝向量矩阵简称为强简化矩阵,只要对算法15作些修改就可实现利用强简化矩阵求出强简化流率基本入树模型的所有反馈环｡下面首先给出一些相关的概念｡定义9删除流率基本入树

19、模型T1(t),T2(t),Tn(t)中全部非重复辅助变量顶点,并仍按原方向连成关联弧的过程称为流率基本入树模型的强简化变换,经强简化变换所得的模型称为原模型的强简化流率基本入树模型｡定义10若一入树T(t)存在P枝入树枝结构相同,则保留其中一枝(若有一枝含有以重复变量顶点vk(t)为根的极大子入树,则保留该枝),该枝称为保留枝,并将其余P-1枝全部删除;同时将与根关联的重复变量顶点vk(t)改为Pvk(t),或将与该根关联的流位变量L(t)改为PL(t),称这一过程为并枝｡定义11在强简化流率基本入树模型中,经过并枝得到枝向量P(Ri(t),Aij

20、(t),Lj(t)｡其中,(Ri(t),Aij(t),Lj(t)为保留枝向量,Ri(t)为流率,Lj(t)为流位,Aij(t)为入树中重复的辅助变量;P为并枝的条数,将它称为枝向量P(Ri(t),Aij(t),Lj(t)的系数,把该枝向量称为强简化流率基本入树的枝向量,简称强简化向量｡其他相关概念请参照文献1,5｡需要说明的是,求强简化流率基本入树模型的反馈环时,反馈环向量前面的系数表示反馈环的条数｡求强简化流率基本入树模型的所有反馈环与求流率基本入树模型的所有反馈环的不同主要是强简化向量含有系数&

21、#65377;实际上可将流率基本入树模型的枝向量的系数看做1,这样就可通过修改原算法得到一个比较通用的算法｡首先将算法2的原步骤(2)的算法改为:单元枝向量a､b的和c的系数(记为coe)等于这两个单元枝向量系数的乘积｡如果coe为1,则c前面不添加数字系数｡然后将算法5的步骤中的第二重循环改为若i>j且loopi,j不是零向量,则分解出loopi,j中的nloop个单元反馈环向量存储在DYLoops中求得DYLoopsk的系数为kcoe;%DYLoopsk为一单元反馈环向量输出有kcoe条j阶的DY

22、Loopsk;通过上述修改,算法15的通用性就增强了,它既可以用来求流率基本入树模型的所有反馈环也可以用来解决强简化流率基本入树模型的所有反馈环的求解｡1.4 复杂性分析算法14是算法5的子函数｡可以看出,算法1和算法2的时间复杂度为O(1);算法3运用了递归的思想,其时间复杂度取决于输入参数a､b,在最坏情况下,若a､b分别含有n1和n2个单元枝向量,则其时间复杂度为O(n1,n2);算法4的时间复杂度为O(k3);算法5通过算法14得到所有反馈环,其运算量的最大项为所以算法5的时间复杂度为O(n(n+1

23、)/2)2)｡3 算例为了便于验证比较,算例采用文献1中的模型｡其中一个典型的模型是“王禾丘能源生态系统工程主导结构流率基本入树模型”｡下面以此模型为例(说明:由于枝向量中的变量均以时间t为自变量,为了实现方便,可只写变量的名称,将t省略)｡3.1 由流率基本入树枝向量矩阵计算所有反馈环根据此模型的流率基本入树可建立起相应的基本入树的枝向量矩阵，把它记为A｡A为五阶矩阵,为了方便，把A记为(aij)5｡那么A为运行命令AllLoop(A)可从输出中得到:1阶反馈环

24、5条;2阶反馈环7条;3阶反馈环6条;4阶反馈环3条;5阶反馈环2条;共有反馈环23条｡利用这些输出就可以很方便地进行流图的反馈分析｡例如,输出的2阶反馈环向量(R2,MF1,MF,RATIO,L1,R1,POR,L2)对应模型中的反馈环为山地薪柴林面积变化量R2(t)(hm2/年)-山地薪柴提供量MF1(t)(kg)+薪柴需求量MF(t)(kg)-实产沼气能占生活用能比RATIO(t)-人口数L1(t)(人)+人口变化量R1(t)(人/年)+人口的环境影响因子POR(t)-山地薪柴林面积L2(t)(hm2)+R2(t)｡3.2 由强简化流率基本入树枝向量矩阵计算所有反馈环此模型的强简化流率基本入树