枝向量矩陣反饋環算法的ＭＡＴＬＡＢ實現

摘要:找到系統中的全部反饋環是系統動力學模型分析的基礎#65377;針對這一問題，比較了幾種常見的反饋環計算方法，得出在時間復雜度上矩陣算法優于行列式算法#65377;在計算機上用MATLAB實現了基于流率基本入樹和強簡化流率基本入樹枝向量矩陣計算所有反饋環的矩陣算法，分析了算法的復雜性，并給出了相應的算例#65377;

關鍵詞:系統動力學;反饋環;枝向量矩陣;MATLAB實現

中圖分類號:TP301.6文獻標志碼:A

文章編號:1001-3695(2007)04-0251-04

系統動力學是通過建立流位#65380;流率系來研究信息反饋系統的一門科學#65377;用系統動力學來處理復雜系統時一個關鍵的問題就是反饋環的計算#65377;一個社會經濟系統往往是由許多反饋環構成的復雜系統，如美國的福瑞斯特等人建立的世界模型Ⅱ就是由82條反饋環組成的#65377;因此，要對此系統的動態性#65380;復雜性進行分析，就必須確定其反饋環的集合，這對于其結構分析#65380;基模確立#65380;模型調試#65380;結果分析都非常重要#65377;針對反饋環的計算也有不少方法，如常見的反饋環計算法有反饋環圖示計算法#65380;枝向量行列式反饋環計算法及枝向量矩陣反饋環計算法等^[1]#65377;這些算法都還未編程，在計算機上難以實現#65377;其困難一方面在于算法本身，有的方法根本就沒辦法用計算機實現;另一方面是算法程序實現上的困難，像枝向量行列式反饋環計算法的行列式計算法則不同于數值為元素的行列式的計算#65377;對于這些問題，本文比較了幾種常見的反饋環計算方法，并給出了運用枝向量矩陣求反饋環算法的MATLAB實現#65377;算例表明，所設計的程序具有很好的應用效果#65377;

1反饋環計算的常見方法比較

常見的反饋環計算方法有反饋環圖示計算法#65380;枝向量行列式反饋環計算法及枝向量矩陣反饋環計算法^[1]#65377;圖示計算法是一種通過計算入樹模型的所有反饋環的圖示法#65377;此方法思路清晰，但當流圖非常復雜時，用此圖示法就感到煩瑣，而且無法用計算機實現，所以效果不佳#65377;下面從時間復雜度上來比較在計算出所有反饋環時的行列式算法和矩陣算法#65377;

為了具有可比性，筆者將行列式算法中枝向量行列式元素的乘法和矩陣算法中枝向量矩陣元素的乘法分別作為兩個算法的基本操作#65377;由行列式算法很容易得到其時間復雜度O(n×n!)#65377;其中n是枝向量行列式的階數#65377;矩陣算法能夠求出當向入樹模型中依次增加入樹時新增的各階反饋環，把入樹模型中所有入樹新增的反饋環加起來即可得到模型中所有的反饋環#65377;由此可以計算出矩陣算法求所有反饋環的時間復雜度為O((n(n+1)/2)2)#65377;其中n為枝向量矩陣的階數#65377;若以矩陣算法的運算量為分子，以行列式算法的運算量為分母，求極限有

可見，矩陣算法的時間復雜度小于行列式算法的的時間復雜度#65377;計算表明，隨著n增大，特別是當n>4時，矩陣算法的優勢更加明顯，這對于算法的計算機實現非常重要#65377;鑒于此給出了矩陣算法的計算機實現#65377;

2枝向量矩陣算法的實現

憑借MATLAB強大的數值計算和字符處理功能，實現了運用矩陣算法計算出入樹模型中的所有反饋環^[2~4]#65377;為了方便表達，用類MATLAB的偽代碼來書寫相應的計算機實現算法#65377;枝向量矩陣有流率基本入樹枝向量矩陣和強簡化流率基本入樹枝向量矩陣^[5]，在矩陣算法的基礎上分別給出由這兩個枝向量矩陣計算所有反饋環的計算機實現算法，下面先給出一些相關的概念#65377;

2.1相關概念

定義1以流率基本入樹T1(t)，T2(t)…，Tn(t)中的枝向量為元素，依次排列構成的向量(Ri(t)，Aij(t)，Lj(t))稱為枝向量#65377;其中Ri(t)為流率，Lj(t)為流位，Aij(t)為入樹中枝的輔助變量依次排列的組合;如果是多枝的變量構成枝向量，Aij(t)隱含有流率與流位#65377;

定義2對于枝向量(Ri(t)，Aij(t)，Lj(t))，若Ri(t)就是Lj(t)所對應的流率，則此枝向量對應流圖中的一個反饋環，把此枝向量稱為反饋環枝向量或反饋環向量#65377;

定義3枝向量加法(Ri(t)，Aij(t)，Lj(t))+(Rt(t)，Atp(t)，Lp(t))僅表示在流率基本入樹中存在(Ri(t)，Aij(t)，Lj(t))和(Rt(t)，Atp(t)，Lp(t))對應的兩枝#65377;

定義4不包含加法的非零枝向量為單元枝向量，幾個單元枝向量做加法所組成的量稱為這幾個枝向量的復合枝向量#65377;以后若不特別聲明，枝向量包括零向量#65380;單元枝向量和復合枝向量#65377;定義5枝向量乘法

(Ri(t)，Aij(t)，Lj(t)，Rj(t)，Atp(t)，Lp(t))Rt(t)是Lj(t)的流率且兩個枝向量中無相同變量

(Rt(t)，Atp(t)，Lp(t)，Rp(t)，Aij(t)，Lj(t))Ri(t)是Lp(t)的流率且兩個枝向量中無相同變量0其他情況

定義6由入樹Ti(t)(包括流率基本入樹及簡化#65380;強簡化流率基本入樹)的枝向量(Ri(t)，Ai1(t)，L1(t))，(Ri(t)，Ai2(t)，L2(t))，…，(Ri(t)，Ain(t)，Ln(t)) (i=1，2，…，n)構成的方陣:

稱為此入樹T1(t)，T2(t)，…，Tn(t)模型的枝向量矩陣#65377;如果矩陣中的元素在入樹中不存在，把它記為零向量#65377;

此枝向量矩陣的乘法法則與代數中給出的數矩陣乘法法則相同#65377;

定義8將入樹模型枝向量矩陣A的全體aii(i=1，2，…，n)置為零，其他aij(i≠j)不變，所得矩陣稱為此入樹模型的對角置零枝向量矩陣#65377;將枝向量矩陣A的最后一行全體元素置為0，得到的矩陣記為A#65377;A稱為A的末行置零矩陣#65377;

2.2由流率基本入樹枝向量矩陣計算所有反饋環的算法

通過流率基本入樹建模法^[6]可以建立起模型的流率基本入樹，按照上述定義就可以得到相應的枝向量矩陣#65377;要通過枝向量矩陣計算出模型的所有反饋環，就要實現枝向量矩陣的乘法，而這就必須首先要實現枝向量矩陣中的元素(枝向量)的加法和乘法#65377;

由定義3可以給出兩個枝向量做加法的具體算法#65377;

算法1c=VPlus(a，b)

輸入:枝向量a和b;

輸出:枝向量a和b做加法的和c#65377;

(1)如果a為零向量，則c=b， 轉(4);否則轉(2)#65377;

(2)如果b為零向量，則c=a， 轉(4);否則轉(3)#65377;

(3)把a和b構成復合枝向量c，轉(4)#65377;

(4)結束#65377;

由定義4設計了兩個枝向量(單元枝向量或零向量)做乘法的具體算法#65377;

算法2c=DYMultiply(a，b)

輸入:單元枝向量a，b;

輸出:做枝向量乘法的結果c#65377;

(1)如果a為零向量或b為零向量則輸出結果c為零向量，轉(4);否則轉(2)#65377;

(2)若b中的第一個流率恰是a中最后一個流位所對應的流率且a#65380;b中無相同變量，則將a中的變量依次放在b中的變量前面構成枝向量c中的變量，轉(4);否則轉(3)#65377;

(3)若a中的第一個流率恰是b中最后一個流位所對應的流率且a#65380;b中無相同變量，則將b中的變量依次放在a中的變量前面構成c中的變量，轉(4);否則c為零向量，轉(4)#65377;

(4)結束#65377;

由于枝向量矩陣中的元素可能為復合枝向量，對于任意兩個枝向量之間的乘法可通過下面的具體算法實現:

算法3c=Multiply(a，b)

輸入:枝向量a，b;

輸出:枝向量c#65377;

(1)如果a#65380;b都是單元枝向量或零向量，則c=DYMultiply(a，b)，轉(5);否則轉(2)#65377;

(2)如果a為單元枝向量或零向量而b為復合枝向量，則①將b分為兩部分，記為b1和b2#65377;其中，b1為b中的第一個單元枝向量，b2為b中除了第一個單元枝向量之外的其他枝向量構成的枝向量#65377;

②c1=DYMultiply(a，b1)，c2=Multiply(a，b2)#65377;

③c=VPlus(c1，c2)，轉(5)#65377;否則，轉(3)#65377;

(3)若a為復合枝向量而b為單元枝向量或零向量，則參照(2)中的①~③作同樣處理;否則，轉(4)#65377;

(4)此時a和b都為復合枝向量，對它們進行如下處理:

①將b分為兩部分，記為b1和b2#65377;其中，b1為b中的第一個單元枝向量，b2為b中除了第一個單元枝向量之外的其他枝向量構成的枝向量#65377;將a用同樣的方法分為兩部分，分別記為a1#65380;a2#65377;②c1=DYMultiply(a1，b1)，c2=Multiply(a1，b2);

c3=Multiply(a2，b1)，c4=Multiply(a2，b2)#65377;

③將c1#65380;c2#65380;c3和c4做枝向量加法得到c，轉(5)#65377;

(5)結束#65377;

實現了枝向量的基本運算后，就可將它們用到枝向量矩陣的乘法運算中#65377;下面分兩步來實現由枝向量矩陣計算出入樹模型的所有反饋環，首先給出向入樹模型中增加入樹時新增反饋環的算法，然后把所有入樹新增的反饋環匯總得到并輸出所有的反饋環#65377;

(1)當向入樹模型中增加入樹時，同時也會相應地增加反饋環的數量#65377;下面給出當增加一棵入樹時新增反饋環的具體算法#65377;

算法4result=CreatLoop(A，k)

輸入:枝向量矩陣A和A的左上角子枝向量矩陣的階數k;

輸出:result為1×k的枝向量矩陣，用來保存計算得到的1~k階反饋環#65377;

①result=cell(1，k); %初始化result

②若k==1，則result{1，1}=A{1，1}，轉⑨;否則轉③#65377;

③result{1，1}=A{k，k};%對角上的元素構成一階反饋環

④%作A的左上角k階枝向量矩陣Ak

%把枝向量矩陣Ak的對角置零

⑥%取Ak的第k行元素構成行枝向量矩陣Xk

⑦%做Ak的末行置零矩陣A0

A0=Ak; %初始化A0

for i=1:k;A0{k，i}=0;end

⑧%做第1~n-1階枝向量矩陣乘法得到2~k階新增的反饋環#65377;

for i=2:k

%初始化Yk用來保存計算出的相應反饋環

Yk=cell(1，k); for p=1:k; Yk{1，p}=0; end

%把Xk和A0做枝向量矩陣乘法

%Xk中只有最后一個元素含有反饋環

⑨結束#65377;

(2)將所有的入樹依次添加到模型中，通過調用算法4就可以計算出所有的反饋環#65377;計算并輸入樹模型中的所有反饋環的具體算法如下:

算法5AllLoop(A)

輸入:枝向量矩陣A;

輸出:枝向量矩陣A所能構成的所有的反饋環#65377;

①n=length(A); 

loop=cell(n，n); 

%loop用來存放枝向量矩陣A所能構成的所有反饋環

②for i=1:n

NewLoop=CreatLoop(A，i);

%把NewLoop賦給矩陣loop中的第i行

for j=1:i; loop{i，j}=NewLoop{1，j};end

end

③%輸出各階反饋環

sum=0;%sum用來統計生成的反饋環的總數

for j=1:n

num=0;%num統計生成的同一階的反饋環的總數

for i=1:n

若i≥j，則

如果loop{i， j}不是零向量，則輸出loop{i， j}中包含的所有j階反饋環，每輸出一個反饋環，使 num=num+1#65377;

end

若num==0，則生成的反饋環中不含有j階反饋環;否則

輸出所有的num條j階反饋環，并使 sum=sum+num #65377;

end

④輸出顯示共有sum條反饋環#65377;

⑤結束#65377;

2.3由強簡化流率基本入樹枝向量矩陣計算所有反饋環的算法

在文獻[1，5]中提到了強簡化流率基本入樹模型，并給出了強簡化流率基本入樹的枝向量矩陣#65377;為了方便，下面將通過強簡化流率基本入樹所得的枝向量矩陣簡稱為強簡化矩陣，只要對算法1~5作些修改就可實現利用強簡化矩陣求出強簡化流率基本入樹模型的所有反饋環#65377;下面首先給出一些相關的概念#65377;

定義9刪除流率基本入樹模型T1(t)，T2(t)，…，Tn(t)中全部非重復輔助變量頂點，并仍按原方向連成關聯弧的過程稱為流率基本入樹模型的強簡化變換，經強簡化變換所得的模型稱為原模型的強簡化流率基本入樹模型#65377;

定義10若一入樹T(t)存在P枝入樹枝結構相同， 則保留其中一枝(若有一枝含有以重復變量頂點vk(t)為根的極大子入樹，則保留該枝)，該枝稱為保留枝， 并將其余 P-1 枝全部刪除;同時將與根關聯的重復變量頂點vk(t)改為Pvk(t)，或將與該根關聯的流位變量L(t)改為PL(t)，稱這一過程為并枝#65377;

定義11在強簡化流率基本入樹模型中，經過并枝得到枝向量P(Ri(t)，Aij(t)，Lj(t))#65377;其中，(Ri(t)，Aij(t)，Lj(t))為保留枝向量，Ri(t)為流率，Lj(t)為流位，Aij(t)為入樹中重復的輔助變量;P為并枝的條數，將它稱為枝向量P(Ri(t)，Aij(t)，Lj(t))的系數，把該枝向量稱為強簡化流率基本入樹的枝向量，簡稱強簡化向量#65377;

其他相關概念請參照文獻[1，5]#65377;需要說明的是，求強簡化流率基本入樹模型的反饋環時，反饋環向量前面的系數表示反饋環的條數#65377;

求強簡化流率基本入樹模型的所有反饋環與求流率基本入樹模型的所有反饋環的不同主要是強簡化向量含有系數#65377;實際上可將流率基本入樹模型的枝向量的系數看做1，這樣就可通過修改原算法得到一個比較通用的算法#65377;首先將算法2的原步驟(2)的算法改為:

①單元枝向量a#65380;b的和c的系數(記為coe)等于這兩個單元枝向量系數的乘積#65377;

②如果coe為1，則c前面不添加數字系數#65377;

然后將算法5的步驟③中的第二重循環改為

若i>j且loop{i，j}不是零向量，則分解出loop{i，j}中的nloop個單元反饋環向量存儲在DYLoops中

求得DYLoops{k}的系數為kcoe; 

%DYLoops{k}為一單元反饋環向量

輸出有kcoe條j階的DYLoops{k};

通過上述修改，算法1~5的通用性就增強了，它既可以用來求流率基本入樹模型的所有反饋環也可以用來解決強簡化流率基本入樹模型的所有反饋環的求解#65377;

2.4復雜性分析

算法1~4是算法5的子函數#65377;可以看出，算法1和算法2的時間復雜度為O(1);算法3運用了遞歸的思想，其時間復雜度取決于輸入參數a#65380;b，在最壞情況下，若a#65380;b分別含有n1和n2個單元枝向量，則其時間復雜度為O(n1，n2);算法4的時間復雜度為O(k3);算法5通過算法1~4得到所有反饋環，其運算量的最大項為

所以算法5的時間復雜度為O((n(n+1)/2)2)#65377;

3算例

為了便于驗證比較，算例采用文獻[1]中的模型#65377;其中一個典型的模型是“王禾丘能源生態系統工程主導結構流率基本入樹模型”#65377;下面以此模型為例(說明:由于枝向量中的變量均以時間t為自變量，為了實現方便，可只寫變量的名稱，將t省略)#65377;

3.1由流率基本入樹枝向量矩陣計算所有反饋環

根據此模型的流率基本入樹可建立起相應的基本入樹的枝向量矩陣，把它記為A#65377;A為五階矩陣，為了方便，把A記為(aij)5#65377;那么A為

運行命令AllLoop(A)可從輸出中得到:1 階反饋環5條;2 階反饋環7條;3 階反饋環6條;4 階反饋環3條;5 階反饋環2條;共有反饋環23條#65377;利用這些輸出就可以很方便地進行流圖的反饋分析#65377;例如，輸出的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)#65377;

3.2由強簡化流率基本入樹枝向量矩陣計算所有反饋環

此模型的強簡化流率基本入樹枝向量矩陣B為

在MATLAB中輸入矩陣B:

運行命令AllLoop(B)可以從輸出中得到經強簡化的反饋環:1 階反饋環5條;2 階反饋環7條;3 階反饋環6條;4 階反饋環3條;5 階反饋環2條#65377;共有反饋環23條#65377;

從輸出中可以看出，強簡化反饋環向量更加簡明而且有的還含有系數#65377;其中強簡化反饋環向量的系數對應著相應數目的反饋環，如2(R3，L1，R1，1，L3)對應著流率基本入樹模型中的相應的兩條2階反饋環為

實際上，強簡化反饋環向量通過還原均可找到相對應的流率基本入樹枝向量矩陣的反饋環向量#65377;

比較可知，第3.1#65380;3.2節的輸出結果與文獻[1]中計算出的反饋環的結果是相同的#65377;文獻[1]是通過圖示法或代數求解的方法得到的，當模型很大時，通過圖示法或代數求解的計算量是可想而知的，有時也是不現實的#65377;通過計算機可快速準確地解決原模型中的反饋環求解問題#65377;

4結束語

在系統動力學中，對于復雜流圖中含有多少個反饋環，有時建模者也算不清，各反饋環中變量之間的具體如何變化就更無從談起，而這點對于建模系統分析#65380;模型調試#65380;結果分析很重要#65377;上述問題的解決對系統動力學的發展很有意義#65377;雖然也有不少計算反饋環的方法，但大都難以或根本沒辦法在計算機上實現，對于一個龐大復雜的模型，其計算量可想而知#65377;筆者從模型的枝向量矩陣出發，給出了可以同時求流率基本入樹模型和強簡化流率基本入樹模型的所有反饋環的比較通用的算法，并在計算機用MATLAB得到了實現，這大大提高了建立和分析模型的效率#65377;

本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文。