枝向量矩陣反饋環(huán)算法的MATLAB實現(xiàn)VIP

1、枝向量矩陣反饋環(huán)算法的MATLAB實現(xiàn)系統(tǒng)動力學(xué)是通過建立流位､流率系來研究信息反饋系統(tǒng)的一門科學(xué)｡用系統(tǒng)動力學(xué)來處理復(fù)雜系統(tǒng)時一個關(guān)鍵的問題就是反饋環(huán)的計算｡一個社會經(jīng)濟系統(tǒng)往往是由許多反饋環(huán)構(gòu)成的復(fù)雜系統(tǒng),如美國的福瑞斯特等人建立的世界模型n就是由82條反饋環(huán)組成的｡因此,要對此系統(tǒng)的動態(tài)性､復(fù)雜性進行分析,就必須確定其反饋環(huán)的集合,這對于其結(jié)構(gòu)分析､基模確立､模型調(diào)試､結(jié)果分析都非常重要｡針對

2、反饋環(huán)的計算也有不少方法,如常見的反饋環(huán)計算法有反饋環(huán)圖示計算法､枝向量行列式反饋環(huán)計算法及枝向量矩陣反饋環(huán)計算法等1｡這些算法都還未編程,在計算機上難以實現(xiàn)｡其困難一方面在于算法本身,有的方法根本就沒辦法用計算機實現(xiàn);另一方面是算法程序?qū)崿F(xiàn)上的困難,像枝向量行列式反饋環(huán)計算法的行列式計算法則不同于數(shù)值為元素的行列式的計算｡對于這些問題,本文比較了幾種常見的反饋環(huán)計算方法,并給出了運用枝向量矩陣求反饋環(huán)算法的MATLA改現(xiàn)｡算例表明,所設(shè)計的程序具有很好的應(yīng)用效果&#65377

3、;1 反饋環(huán)計算的常見方法比較常見的反饋環(huán)計算方法有反饋環(huán)圖示計算法､枝向量行列式反饋環(huán)計算法及枝向量矩陣反饋環(huán)計算法1｡圖示計算法是一種通過計算入樹模型的所有反饋環(huán)的圖示法｡此方法思路清晰,但當(dāng)流圖非常復(fù)雜時,用此圖示法就感到煩瑣,而且無法用計算機實現(xiàn),所以效果不佳｡下面從時間復(fù)雜度上來比較在計算出所有反饋環(huán)時的行列式算法和矩陣算法｡為了具有可比性,筆者將行列式算法中枝向量行列式元素的乘法和矩陣算法中枝向量矩陣元素的乘法分別作為兩個算法的基本操作｡由行列式算法很

4、容易得到其時間復(fù)雜度O(nXn!)｡其中n是枝向量行列式的階數(shù)｡矩陣算法能夠求出當(dāng)向入樹模型中依次增加入樹時新增的各階反饋環(huán),把入樹模型中所有入樹新增的反饋環(huán)加起來即可得到模型中所有的反饋環(huán)｡由此可以計算出矩陣算法求所有反饋環(huán)的時間復(fù)雜度為O(n(n+1)/2)2)｡其中n為枝向量矩陣的階數(shù)｡若以矩陣算法的運算量為分子,以行列式算法的運算量為分母,求極限有可見,矩陣算法的時間復(fù)雜度小于行列式算法的的時間復(fù)雜度｡計算表明,隨著n增大，特別是當(dāng)n>4時，矩陣算法

5、的優(yōu)勢更加明顯,這對于算法的計算機實現(xiàn)非常重要｡鑒于此給出了矩陣算法的計算機實現(xiàn)｡2 枝向量矩陣算法的實現(xiàn)憑借MATLAB雖大的數(shù)值計算和字符處理功能,實現(xiàn)了運用矩陣算法計算出入樹模型中的所有反饋環(huán)24｡為了方便表達,用類MATLAB勺偽代碼來書寫相應(yīng)的計算機實現(xiàn)算法｡枝向量矩陣有流率基本入樹枝向量矩陣和強簡化流率基本入樹枝向量矩陣5,在矩陣算法的基礎(chǔ)上分別給出由這兩個枝向量矩陣計算所有反饋環(huán)的計算機實現(xiàn)算法,下面先給出一些相關(guān)的概念｡2.1 相關(guān)概念定義1以流率基本入樹T1(t),T

6、2(t),Tn(t)中的枝向量為元素,依次排列構(gòu)成的向量(Ri(t),Aij(t),Lj(t)稱為枝向量｡其中Ri(t)為流率,Lj(t)為流位,Aij(t)為入樹中枝的輔助變量依次排列的組合;如果是多枝的變量構(gòu)成枝向量,Aij(t)隱含有流率與流位｡定義2對于枝向量(Ri(t),Aij(t),Lj(t),若Ri(t)就是Lj(t)所對應(yīng)的流率,則此枝向量對應(yīng)流圖中的一個反饋環(huán),把此枝向量稱為反饋環(huán)枝向量或反饋環(huán)向量｡定義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)對應(yīng)的兩枝｡定義4不包含加法的非零枝向量為單元枝向量,幾個單元枝向量做加法所組成的量稱為這幾個枝向量的復(fù)合枝向量｡以后若不特別聲明,枝向量包括零向量､單元枝向量和復(fù)合枝向量｡定義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)構(gòu)成的方陣：稱為此入樹T1(t),T2(t),Tn(t)模型的枝向量矩陣｡如果矩陣中的元素在入樹中不存在,把它記為零向量｡此枝向量矩陣的乘法法則與代數(shù)中給出的數(shù)矩陣乘法法則相同｡定義8將入樹模型枝向量矩陣A的全體aii(i=1,2,n)置為零,其他aij(i中j)不

9、變,所得矩陣稱為此入樹模型的對角置零枝向量矩陣｡將枝向量矩陣A的最后一行全體元素置為0,得到的矩陣記為A｡A稱為A的末行置零矩陣｡2.2 由流率基本入樹枝向量矩陣計算所有反饋環(huán)的算法通過流率基本入樹建模法6可以建立起模型的流率基本入樹,按照上述定義就可以得到相應(yīng)的枝向量矩陣｡要通過枝向量矩陣計算出模型的所有反饋環(huán),就要實現(xiàn)枝向量矩陣的乘法,而這就必須首先要實現(xiàn)枝向量矩陣中的元素(枝向量)的加法和乘法｡由定義3可以給出兩個枝向量做加法的具體算法｡算法1c=VPlus

10、(a,b)輸入:枝向量a和b;輸出:枝向量a和b做加法的和c｡(1) 如果a為零向量,則c=b,轉(zhuǎn)(4);否則轉(zhuǎn)(2)｡(2) 如果b為零向量,則c=a,轉(zhuǎn)(4);否則轉(zhuǎn)(3)｡(3) 把a和b構(gòu)成復(fù)合枝向量c,轉(zhuǎn)(4)｡(4) 結(jié)束｡由定義4設(shè)計了兩個枝向量(單元枝向量或零向量)做乘法的具體算法｡算法2c=DYMultiply(a,b)輸入:單元枝向量a,b;輸出:做枝向量乘法的結(jié)果c｡(1) 如果a為零向量或b為零向量則輸出結(jié)果c為零向量

11、,轉(zhuǎn)(4);否則轉(zhuǎn)(2)｡(2) 若b中的第一個流率恰是a中最后一個流位所對應(yīng)的流率且a､b中無相同變量,則將a中的變量依次放在b中的變量前面構(gòu)成枝向量c中的變量,轉(zhuǎn)(4);否則轉(zhuǎn)(3)｡(3) 若a中的第一個流率恰是b中最后一個流位所對應(yīng)的流率且a､b中無相同變量,則將b中的變量依次放在a中的變量前面構(gòu)成c中的變量,轉(zhuǎn)(4);否則c為零向量,轉(zhuǎn)(4)｡(4) 結(jié)束｡由于枝向量矩陣中的元素可能為復(fù)合枝向量,對于任意兩個枝向量之間的乘法可通過下面的具體算法實現(xiàn):算法

12、3c=Multiply(a,b)輸入:枝向量a,b;輸出:枝向量c｡(1) 如果a､b都是單元枝向量或零向量,則c=DYMultiply(a,b),轉(zhuǎn)(5);否則轉(zhuǎn)(2)｡(2)如果a為單元枝向量或零向量而b為復(fù)合枝向量,則將b分為兩部分,記為bl和b2｡其中,b1為b中的第一個單元枝向量,b2為b中除了第一個單元枝向量之外的其他枝向量構(gòu)成的枝向量｡c1=DYMultiply(a,b1),c2=Multiply(a,b2)｡c=VPlus(c1,c2),轉(zhuǎn)(5)&

13、amp;#65377;否則,轉(zhuǎn)(3)｡(3)若a為復(fù)合枝向量而b為單元枝向量或零向量,則參照(2)中的D作同本效t理；否則，轉(zhuǎn)(4)｡(4) 此時a和b都為復(fù)合枝向量,對它們進行如下處理:將b分為兩部分,記為bl和b2｡其中,b1為b中的第一個單元枝向量,b2為b中除了第一個單元枝向量之外的其他枝向量構(gòu)成的枝向量｡將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,轉(zhuǎn)(5)｡(5) 結(jié)束｡實現(xiàn)了枝向量的基本運算后,就可將它們用到枝向量矩陣的乘法運算中｡下面分兩步來實現(xiàn)由枝向量矩陣計算出入樹模型的所有反饋環(huán),首先給出向入樹模型中增加入樹時新增反饋環(huán)的算法,然后把所有入樹新增的反饋環(huán)匯總得到并輸出所有的反饋環(huán)｡(1) 當(dāng)向入樹模型中增加入樹時,同時也會相應(yīng)地增加反饋環(huán)的數(shù)量｡下面給出當(dāng)增加一

15、棵入樹時新增反饋環(huán)的具體算法｡算法4result=CreatLoop(A,k)輸入:枝向量矩陣A和A的左上角子枝向量矩陣的階數(shù)k;輸出:result為1xk的枝向量矩陣，用來保存計算得到的1k階反饋環(huán)｡result=cell(1,k);%初始化result若k=1,則result1,1=A1,1,轉(zhuǎn);否則轉(zhuǎn)｡result1,1=Ak,k;%對角上的元素構(gòu)成一階反饋環(huán)%乍A的左上角k階枝向量矩陣Ak%把枝向量矩陣Ak的對角置零 %取Ak的第k行元素構(gòu)成行枝向量矩陣Xk %做Ak的末行置零矩陣A0A0=Ak;%初始化A0fori=1

16、:k;A0k,i=0;end峨第1n-1階枝向量矩陣乘法得到2k階新增的反饋環(huán)｡fori=2:k%初始化Yk用來保存計算出的相應(yīng)反饋環(huán)Yk=cell(1,k);forp=1:k;Yk1,p=0;end?r%ffiXk和A0做枝向量矩陣乘法%Xk只有最后一個元素含有反饋環(huán)結(jié)束｡(2) 將所有的入樹依次添加到模型中,通過調(diào)用算法4就可以計算出所有的反饋環(huán)｡計算并輸入樹模型中的所有反饋環(huán)的具體算法如下:算法5AllLoop(A)輸入:枝向量矩陣A;輸出:枝向量矩陣A所能構(gòu)成的所有的反饋環(huán)｡n=length(A);

17、loop=cell(n,n);%loop用來存放枝向量矩陣A所能構(gòu)成的所有反饋環(huán)fo門=1:n?rNewLoop=CreatLoop(A,i);?r%ffiNewLoop賦給矢!陣loop中的第i行forj=1:i;loopi,j=NewLoop1,j;endend%俞出各階反饋環(huán)?rsum=0;%sum用來統(tǒng)計生成的反饋環(huán)的總數(shù)?rforj=1:nnum=0;%nurmg計生成的同一階的反饋環(huán)的總數(shù)fori=1:n若i刁，則如果loopi,j不是零向量,則輸出loopi,j中包含的所有j階反饋環(huán),每輸出一個反饋環(huán),使num=num+1｡end若num=0,則生成的反饋環(huán)

18、中不含有j階反饋環(huán)；否則輸出所有的num條j階反饋環(huán)，并使sum=sum+nun｡end輸出顯示共有sum條反饋環(huán)｡結(jié)束｡1.3 由強簡化流率基本入樹枝向量矩陣計算所有反饋環(huán)的算在文獻1,5中提到了強簡化流率基本入樹模型,并給出了強簡化流率基本入樹的枝向量矩陣｡為了方便,下面將通過強簡化流率基本入樹所得的枝向量矩陣簡稱為強簡化矩陣,只要對算法15作些修改就可實現(xiàn)利用強簡化矩陣求出強簡化流率基本入樹模型的所有反饋環(huán)｡下面首先給出一些相關(guān)的概念｡定義9刪除流率基本入樹

19、模型T1(t),T2(t),Tn(t)中全部非重復(fù)輔助變量頂點,并仍按原方向連成關(guān)聯(lián)弧的過程稱為流率基本入樹模型的強簡化變換,經(jīng)強簡化變換所得的模型稱為原模型的強簡化流率基本入樹模型｡定義10若一入樹T(t)存在P枝入樹枝結(jié)構(gòu)相同,則保留其中一枝(若有一枝含有以重復(fù)變量頂點vk(t)為根的極大子入樹,則保留該枝),該枝稱為保留枝,并將其余P-1枝全部刪除;同時將與根關(guān)聯(lián)的重復(fù)變量頂點vk(t)改為Pvk(t),或?qū)⑴c該根關(guān)聯(lián)的流位變量L(t)改為PL(t),稱這一過程為并枝｡定義11在強簡化流率基本入樹模型中,經(jīng)過并枝得到枝向量P(Ri(t),Aij

20、(t),Lj(t)｡其中,(Ri(t),Aij(t),Lj(t)為保留枝向量,Ri(t)為流率,Lj(t)為流位,Aij(t)為入樹中重復(fù)的輔助變量;P為并枝的條數(shù),將它稱為枝向量P(Ri(t),Aij(t),Lj(t)的系數(shù),把該枝向量稱為強簡化流率基本入樹的枝向量,簡稱強簡化向量｡其他相關(guān)概念請參照文獻1,5｡需要說明的是,求強簡化流率基本入樹模型的反饋環(huán)時,反饋環(huán)向量前面的系數(shù)表示反饋環(huán)的條數(shù)｡求強簡化流率基本入樹模型的所有反饋環(huán)與求流率基本入樹模型的所有反饋環(huán)的不同主要是強簡化向量含有系數(shù)&

21、#65377;實際上可將流率基本入樹模型的枝向量的系數(shù)看做1,這樣就可通過修改原算法得到一個比較通用的算法｡首先將算法2的原步驟(2)的算法改為:單元枝向量a､b的和c的系數(shù)(記為coe)等于這兩個單元枝向量系數(shù)的乘積｡如果coe為1,則c前面不添加數(shù)字系數(shù)｡然后將算法5的步驟中的第二重循環(huán)改為若i>j且loopi,j不是零向量,則分解出loopi,j中的nloop個單元反饋環(huán)向量存儲在DYLoops中求得DYLoopsk的系數(shù)為kcoe;%DYLoopsk為一單元反饋環(huán)向量輸出有kcoe條j階的DY

22、Loopsk;通過上述修改,算法15的通用性就增強了,它既可以用來求流率基本入樹模型的所有反饋環(huán)也可以用來解決強簡化流率基本入樹模型的所有反饋環(huán)的求解｡1.4 復(fù)雜性分析算法14是算法5的子函數(shù)｡可以看出,算法1和算法2的時間復(fù)雜度為O(1);算法3運用了遞歸的思想,其時間復(fù)雜度取決于輸入?yún)?shù)a､b,在最壞情況下,若a､b分別含有n1和n2個單元枝向量,則其時間復(fù)雜度為O(n1,n2);算法4的時間復(fù)雜度為O(k3);算法5通過算法14得到所有反饋環(huán),其運算量的最大項為所以算法5的時間復(fù)雜度為O(n(n+1

23、)/2)2)｡3 算例為了便于驗證比較,算例采用文獻1中的模型｡其中一個典型的模型是“王禾丘能源生態(tài)系統(tǒng)工程主導(dǎo)結(jié)構(gòu)流率基本入樹模型”｡下面以此模型為例(說明:由于枝向量中的變量均以時間t為自變量,為了實現(xiàn)方便,可只寫變量的名稱,將t省略)｡3.1 由流率基本入樹枝向量矩陣計算所有反饋環(huán)根據(jù)此模型的流率基本入樹可建立起相應(yīng)的基本入樹的枝向量矩陣，把它記為A｡A為五階矩陣,為了方便，把A記為(aij)5｡那么A為運行命令A(yù)llLoop(A)可從輸出中得到:1階反饋環(huán)

24、5條;2階反饋環(huán)7條;3階反饋環(huán)6條;4階反饋環(huán)3條;5階反饋環(huán)2條;共有反饋環(huán)23條｡利用這些輸出就可以很方便地進行流圖的反饋分析｡例如,輸出的2階反饋環(huán)向量(R2,MF1,MF,RATIO,L1,R1,POR,L2)對應(yīng)模型中的反饋環(huán)為山地薪柴林面積變化量R2(t)(hm2/年)-山地薪柴提供量MF1(t)(kg)+薪柴需求量MF(t)(kg)-實產(chǎn)沼氣能占生活用能比RATIO(t)-人口數(shù)L1(t)(人)+人口變化量R1(t)(人/年)+人口的環(huán)境影響因子POR(t)-山地薪柴林面積L2(t)(hm2)+R2(t)｡3.2 由強簡化流率基本入樹枝向量矩陣計算所有反饋環(huán)此模型的強簡化流率基本入樹