基于滿足性判定的布爾網絡環求解算法

郭文生，楊國武，李曉瑜，高 敏

·計算機工程與應用·

基于滿足性判定的布爾網絡環求解算法

郭文生1，楊國武1，李曉瑜1，高 敏2

(1. 電子科技大學計算機科學與工程學院 成都 610054；2. 加州大學洛杉磯分校電子工程學院 美國 加州 洛杉磯 CA 90034)

環是布爾網絡狀態轉換過程中的穩定態，在模式檢測、基因調控網絡和可達性分析等領域都有重要的意義。計算布爾網絡狀態轉換中的所有環是一個NP完全問題。該文基于全解布爾滿足性判定(SAT)算法，設計了一種求解所有小于等于指定步長環的算法。算法基于布爾網絡的狀態轉換函數和狀態環屬性生成合取范式形式(CNF)的問題集，通過融合沖突子句學習(CDCL)、非時序回退、阻塞子句和變量分類等技術，降低算法的計算復雜度。實驗結果表明，該算法能夠高效地計算指定步長的環。對于無法計算所有環的復雜網絡，指定步長計算環的方式將更有應用價值。

布爾網絡; 環; 滿足性判定; 阻塞子句

布爾網絡的穩定狀態包括兩種：不動點和環。不動點是環的一種特例，即環步長為1。布爾網絡環的計算是一個NP完全問題[1-2]。近年來，隨著布爾滿足性判定SAT算法的不斷優化，越來越多的布爾網絡問題都可以用SAT算法求解，如無邊界的符號模式檢測、可達性分析、量詞消除、基因調控網絡求環等[3-8]，但對于復雜布爾網絡的所有環的求解仍然具有非常高的計算復雜度。本文的目標是基于SAT全解算法計算所有小于指定步長的環。

1 相關技術

1.1 基本術語

使用SAT算法求解，絕大部分應用都需要轉換成通用的SAT問題描述形式，即CNF范式。CNF范式由所有子句的合取形式構成，其中每個子句是若干文字(變量或非變量)的析取。每個文字是變量本身或者是變量的非。CNF范式的格式為：

式中，F表示要求解的滿足性問題；Ci表示子句；N表示范式中包含的子句的數量；lk表示文字，lk∈{v,?v}(v是一個變量)；k表示子句中包含文字的數量(k≤n，n為變量的數量)。

1.2 基于沖突子句學習和非時序回退的DPLL算法

Davis-Putnam-Longeman-Loveland(DPLL)算法是完全SAT求解算法的基本架構。幾乎它是所有現代的完全求解算法的基礎，包括單文字子句消除、肯定-否定規則、原子公式消除等boolean constraint propagation(BCP)規則。而CDCL和非時序回退的技術能夠有效的降低搜索空間，進而降低SAT算法的時間和空間復雜度[9]。

BCP算法實現子句單元廣播，即根據已賦值的變量推導未賦值的變量值。如果BCP推導時出現沖突，則調用CDCL算法進行沖突學習并生成沖突子句。沖突子句可以避免重復的路徑搜索。CDCL算法基于蘊含圖(implication graph)技術實現沖突子句的學習和非時序的決策級回退。蘊含圖描述在搜索過程中的部分變量賦值及蘊含關系，選擇不同的決策變量選擇和不同的BCP順序將生成不同的蘊含圖，因此搜索過程中蘊含圖是動態變化的。

1.3 基于單解求解器的全解求解算法

計算滿足性判定問題F的所有滿足解，可以基于單解求解器通過如下步驟實現：

1) 判斷F是否滿足。若不滿足，則輸出不滿足；否則，計算一個滿足解x1, x2,…,xn，n是問題中的變量數。

2) 用vi表示滿足解中的一個變量賦值。若解xi=0，則對應的變量表示為?vi；否則表示為vi。因此，一個滿足解可以表示為由變量或變量的非形成的一個合取子句s。

3) 將滿足解形成的合取子句取反，即所有變量取反并形成一個析取子句c。

4) 將析取子句c加入到求解的問題集F，形成新的問題集FFc=∧。

5) 判斷新的問題集F′是否存在滿足解。若存在，則轉到步驟2)；否則輸出所有已計算的滿足解。

上述求全解的過程存在如下問題：

1) 每一次滿足解的計算都是以整個問題為基礎，對前一次求解過程中的搜索空間沒有記錄。

2) 沒有根據求全解的特點進行搜索算法優化。

針對上述問題，文獻[3]提出了阻塞子句(blocking clause)的方式。當單解SAT求解器獲得一個滿足解時，并不停止求解過程，而是將該滿足解對應的變量合取表示方式取反形成一個析取子句作為阻塞子句加到問題集中，繼續求解過程。通過這種方式可以保留前一次求解過程中所學習的沖突子句，并不再搜索已經搜索過的路徑空間。然而，由于滿足性問題的滿足解可能非常多，從而導致阻塞子句需要較大的內存空間，并且大量的阻塞子句也會增加問題的求解復雜度。文獻[4]為降低阻塞子句的數量，使用子句聚合方式對產生的阻塞子句集進行蘊含化簡，縮短阻塞子句的長度，進而降低阻塞子句的數量。文獻[10]提出了基于算法避免再次獲得重復滿足解的方式，降低內存開銷并保持求解復雜度不變。它沒有使用阻塞子句的方式避免計算重復解，而是基于決策變量與求得解的蘊含關系生成沖突子句，進而避免后續的計算出現重復的解。這種生成沖突子句的方式在本質上等同于生成阻塞子句的方式。同時，它提出了求解問題中變量劃分的方法，它將要求解問題中的變量劃分成兩類：重要變量集和不重要變量集。在搜索過程中先選擇重要變量進行決策賦值，當重要變量全部賦值完成后才對不重要變量進行決策賦值。然而該算法是基于CNF范式的求解，其重要變量和非重要變量沒有結合具體問題進行定義。文獻[11]為降低阻塞子句的數量并降低后續計算的搜索空間使用最短阻塞子句的方式計算全解，它基于覆蓋算法最小化滿足片段減少阻塞子句中的文字數，但該算法使用的覆蓋算法本身復雜度也比較高，對于問題子句數和變量數多的求解問題計算復雜度較大。

上述求全解算法只考慮了如何避免已求得的滿足解的問題，而計算不動點、基因調控網絡吸引區求解等問題在計算過程中不能僅避免已求得的滿足解，還需要避免已經計算的不動點或狀態環。

2 布爾網絡環求解算法

2.1 基于狀態轉換函數和環屬性創建求解問題

布爾網絡描述的是變量之間的相互依賴關系，每個變量的值只能是0或1。變量的下一個時刻的值由依賴節點的當前時刻的值確定，即：

式中，vi′表示節點i的下一個時刻的狀態值；v1，v2,…,vm表示節點i依賴的m個節點的當前時刻的值。

布爾網絡的環是布爾網絡狀態轉換過程中的穩定狀態。如果環的步長為N，則布爾網絡至少需要經過N個時刻的狀態遷移才能進入環狀態。因此，根據布爾網絡的狀態轉換函數，可以通過N步展開的方式生成包含N+1個布爾網絡狀態的轉換路徑，即：

式中，path表示布爾網絡狀態轉換圖中的一條路徑；s0表示網絡的初始狀態；sk=vk1,vk2,…,vkn表示第k時刻網絡的狀態。

同時，步長小于等于N的所有環的屬性p可以描述成為：

根據式(2)～式(4)可以得到計算步長小于等于N的環的求解問題集，并可以通過引入輔助變量的方式將其轉換成CNF形式。

2.2 基于SAT全解算法的環計算

基于阻塞子句、沖突蘊含圖生成阻塞子句和變量分類等SAT全解算法的實現技術，結合布爾網絡轉換關系和環屬性描述方式，本文提出了求解步長小于等于N的環的求解算法。其算法的基本流程為：Input: CNF范式F，變量集V，環步長N；

Output：全部滿足解all-solution或不滿足UNSAT；初始化:

L=NULL; //L是確定賦值的文字集合；

DecisionLevel = 0; //DecisionLevel是搜索的深度；Status = SAT;

生成CNF問題集F： cisionVars(F);

IF (BCP(F,L)==CONFLICT) THEN Return UNSAT；WHILE (Status == SAT)

Begin

preClauseCount =sizeof()F; //獲得原始問題的子句數量；

WHILE (sizeof(L) ＜ sizeof(V)) //sizeof()獲得變量的數量；

Begin

l=PickBranchLit(F,VVar(L)); //從未賦值的變量中選擇一個變量賦值，Var(L)是已賦值的變量集合；

Add(l, L); DecisionLevel++; //增加賦值變量到解集L；

IF (BCP(F, L) == CONFLICT)

Begin

DecisionLevel = CDCL(F, L);

If (DecisionLevel ＜ curDecisionLevel)

Begin

removeClausesFrom(preClauseCount);

c=GenBlockClause(D);

FFc=∧;

End

IF (DecisionLevel==0) Status = UNSAT, Break;

ELSE BackJump(DecisionLevel);

End

End

curDecisionLevel = DecisionLevel;

s =L;//獲得一個滿足解；

c=GenBlockClause(D);

FFc=∧;

BackJump(DecisionLevel ? 1);

End

IF (Sizeof(All-Solution) ＞ 0) Return All-Solution;

ELSE Return UNSAT;

PickBranchLit( )函數實現決策變量的選擇并為其賦值。根據轉換CNF過程中存在大量輔助變量的特性，選擇決策變量時采用了變量分類的方式減少決策變量的選擇空間。setDecisionVars( )函數用于設置變量的屬性，即該變量是否可以作為決策變量。GenBlockClause( )函數根據當前的決策變量賦值獲得阻塞子句。參數D表示決策級為k時的所有決策變量賦值對應的文字集合D={d1, d2,…,dk}，則GenBlockClause( )函數生成的阻塞子句為c= ?d1∨?d2∨…∨?dk。

2.3 基于變量分類的決策文字選擇優化

為提高搜索速度并降低搜索空間，本文將待選擇變量分為兩類：決策變量和蘊含變量。決策變量是布爾網絡狀態轉換過程中對應的狀態變量，蘊含變量是布爾網絡轉換成CNF范式過程中引入的輔助變量。本文中布爾網絡的描述形式采用的是類似于berkeley logic interchange format (BLIF)的格式，轉換關系示例如表1所示。

表1 布爾網絡轉換關系描述示例

表1能夠準確的描述布爾網絡節點之間的依賴關系，并根據該依賴關系得到對應的布爾網絡狀態轉換函數。在同步布爾網絡中，給定一個布爾網絡的轉換步數N后，可以根據轉換關系展開得到布爾網絡的狀態轉換圖，并將其轉換成CNF形式。然而在轉換過程中將會引入大量的輔助變量。表1中的每一個轉換定義項都需要引入一個輔助變量。下面以表1中節點2的下一個狀態轉換函數為例描述轉換過程：

1)若節點的當前狀態為1，則表示為vi( i∈{1,2,3})，否則表示為?vi。

3) 根據表1中(.n 2 2 1 3)的描述，其轉換函數v2′=f( v1, v3)，同時其后的兩行分別表示v2′=?v1∧v3=1和v2′=v1∧?v3=1，即節點2的轉換關系函數可表示為：

4) 針對轉換關系函數中的每一個合取項引入一個輔助變量y1, y2，并令y1=?v1∧v3,y2=v1∧?v3，則節點2的轉換關系函數可表示為：

5)步驟4)的轉換關系轉換成CNF范式的子句集為：

根據上述轉換過程可以看出，輔助變量的賦值是由節點狀態值所蘊含的，因此在SAT求解過程中只要確定了節點狀態值，則輔助變量的賦值能夠根據蘊含關系確定，因此在PickBranchLit( )函數選擇決策變量的過程中不需要選擇輔助變量作為決策變量，進而降低SAT求解的搜索空間。

2.4 基于蘊含圖的決策變量生成阻塞子句優化

根據計算全解和增加子句的過程，可將要求解的問題分成兩部分：原始的問題F和阻塞子句集B。SAT問題F的一個滿足解si( x1, x2,…,xn)可以用合取子句ci=l1∧l2∧…∧ln表示，如果xi=1，則子句中對應的文字li=vi；否則li=?vi。因此，所有已獲得的m個滿足解的集合將形成一個析取disjunctive normal form(DNF)范式，其格式為：

基于滿足解集，可以得到阻塞子句集BS=?。在計算是否存在1m+個滿足解時，需要將B加入到問題集F中，即FFB=∧，以避免獲得重復的解。顯然，S中包含的子句越多且每個子句中的文字數越多，則計算1m+個滿足解的時間和空間復雜度越大。

優化的目標之一是降低阻塞子句中文字的數量。沖突子句學習的SAT求解算法會在求解過程中基于選擇的決策變量進行蘊含推理，進而生成局部蘊含圖。蘊含圖示例如圖1所示。

圖1中待求解的CNF范式由4個子句構成F=C1∧C2∧C3∧C4。xi=a@l表示決策級為l時xi的值為a。子句數據庫為C1=x1∨x2，C2=?x2∨?x3，C3=x3∨x4∨x5，C4=?x5∨x6。決策級為1時，決策變量及分配的值為{x1=0@1}，蘊含分配的值為{x2=1@1, x3=0@1}。決策級為2，決策變量及分配的值為{x4=0@2}，蘊含分配的值為{x5=1, x6=1}。最后獲得可滿足解{x1=0, x2=1,x3=0, x4=0, x5=1, x6=1}。

圖1 SAT求解過程中的蘊含圖示例

根據圖1可以看到，當前的可滿足解s={?v1, v2,?v3,?v4, v5, v6}是由兩個決策變量的賦值d={x1=0@1,x4=0@2}確定的。因此，為避免下一次計算出現重復的解，則生成的阻塞子句是?v1∧?v4的非，即v1∨v4。顯然，使用決策變量生成阻塞子句的方式極大的降低了阻塞子句中文字的數量。

優化目標之二是減少阻塞子句的數量。在基于沖突子句的SAT求解過程中，當獲得一個滿足解時，采用時序回退的算法計算下一個滿足解。即全解搜索方式采用深度優先搜索算法。在深度優先搜索時，如果沖突回退的決策級比當前的決策級小，則表明當前決策級之下已經不存在滿足解，則可以刪除當前產生的所有阻塞子句，并基于回退的決策級backlevel生成一個阻塞子句c=?d1∨?d2∨…∨?dbacklevel，表示backlevel之后的所有滿足解都已經被計算，并使用一個阻塞子句c避免所有搜索過的滿足解空間。

3 實驗結果分析

基于Minisat求解器[12]實現了不動點和步長小于等于K的周期性環的SAT全解算法period cycles algorithm based on decision node and important variable(PC-DI)。并選擇了一個基于SAT計算所有環的求解器BNS[5]作為比較工具。實驗使用與文獻[5]的測試案例和N-K隨機布爾網絡作為測試基準。N-K網絡使用R環境中的BoolNet包隨機生成8個最大環長不超過100的布爾網絡。測試平臺的配置為Intel Xeon CPU@3.3 GHz 8-Core(測試時只用了一個核)，內存為128 GB，操作系統是Ubuntu 12.04。

實驗分析了兩方面的性能：所有周期性環求解性能和指定步長周期性環求解性能。

為了使用本文的算法計算所有的周期環，首先使用BNS計算出所有環并獲得最大的環長，然后以最大環長為指定步長運行本文的求解器，測試結果如表2所示。

表2 所有環求解計算結果

圖2 指定步長環的運行時結果

實驗結果顯示本文算法在計算環的最大步長小于100的真實案例和隨機案例中都具有較大的性能優勢。

基于表2的測試結果，選擇其中4個BNS計算時間在100 s以上的案例(Random2, Random4, Random6, Random8)進行指定步長環計算的分析。對每一個測試案例指定不同的步長進行測試。指定的步長初始值為1，每次增加4個步長，最大指定步長設置為69(大于所有測試案例的環的最大步長)。測試結果如圖2和圖3所示。

圖2比較的是指定最大環步長時計算所有周期環的時間，圖3比較的是指定最大環步長時計算出的周期環的數量。Random2、Random4、Random6和Random8的所有周期環分別是5 632個、24個、2 280個和478個。圖2顯示Random4和Random8的運行時間增加比較平緩，Random2和Random6分別在指定步長為61和69時運行時間增加較大。結合圖3可以看出，Random2和Random6分別在指定步長為61和69時周期性環數顯著增加。

圖3 小于等于指定步長環的總環數

實驗結果顯示，本文算法對于環長分布比較均勻的布爾網絡具有較低的計算復雜度，并且指定的步長與運行時間成線性增加的關系。同時，對于運行時間較長且環長分布不均勻的案例(步長大的環占所有環的比例較高)，其主要的時間開銷集中在步長大的周期環計算過程中。

基于實驗結果可以看到，當指定的步長小于最大狀態環步長時，算法計算復雜度較低。因此本文算法可以通過指定步長的方式求解更復雜的布爾網絡的狀態環。對于無法在有限的時間內計算所有周期環的案例，本文將具有更大的優勢。

4 結 束 語

本文提出了一種高效的周期性狀態環求解算法。為降低決策變量的選擇空間提升搜索速度，該算法將布爾網絡轉換成CNF過程中引入的輔助變量標記為非決策變量。通過引入基于蘊含圖技術建立阻塞子句的方式，該算法能夠降低阻塞子句的數量并實現非時序題集的方式，可以充分利用SAT全解求解算法計算狀態環。實驗結果表明，該算法具有較低的計算復雜度，并且能夠對復雜的網絡進行指定步長的狀態環求解。

[1] ZHAO Q. A remark on "scalar equations for synchronous boolean networks with biological applications" by C. Farrow, J. Heidel, J. Maloney, and J. Rogers[J]. IEEE Transactions on Neural Networks, 2005, 16(6): 1715-1716.

[2] T AKUTSU, S KOSUB, A A MELKMAN, et al. Finding a periodic attractor of a Boolean network[J]. Transactions on Computational Biology and Bioinformatics, 2012, 9(5): 1410-1421.

[3] MCMILLAN K L. Applying SAT methods in unbounded symbolic model checking[C]//14th International Conference on Computer Aided Verification. Denmark: Springer, 2002.

[4] CHAUHAN P, CLARKE E M, KROENING D. Using SAT based image computation for reachability analysis[M]. The Netherlands: Springer, 2003.

[5] DUBROVA E, TESLENKO M. A SAT-based algorithm for finding attractors in synchronous Boolean networks[J]. Transactions on Computational Biology and Bioinformatics (TCBB), 2011, 8(5): 1393-1399.

[6] ZHENG De-sheng, YANG Guo-wu, LI Xiao-yu, et al. An efficient algorithm for computing attractors of synchronous and asynchronous Boolean networks[EB/OL]. [2014-01-01] http://journals.plos.org/plosone/article?id=10.1371/journal.p one.0060593.

[7] GUO W, YANG G, WU W, et al. A parallel attractor finding algorithm based on Boolean satisfiability for genetic regulatory networks[EB/OL]. [2014-01-03] http://journals. plos.org/plosone/article?id=10.1371/journal.pone.0094258.

[8] BRAUER J, KING A, KRIENER J. Existential quantification as incremental SAT[C]//23rd International Conference on Computer Aided Verification. Snowbird UT , USA: Springer, 2011.

[9] LIMA J F. Boolean satisfiability solvers: Techniques, implementations and analysis[J]. Electrónica e Telecomunica??es S. A., 2013, 5(2): 218-225.

[10] GRUMBERG O, SCHUSTER A, YADGAR A. Memory efficient all-solutions SAT solver and its application for reachability analysis[C]//5th international confrence on Formal Methods in Computer-Aided Design. Texas, USA: Springer, 2004.

[11] YU Y, SUBRAMANYAN P, TSISKARIDZE N, MALIK S. All-SAT using minimal blocking clauses[C]//27th International Conference on VLSI Design and 13th International Conference on Embedded Systems. Mumbai, India: IEEE, 2014.

[12] EEN N, SORENSSON N. Minisat 2.2.0[CP/DK]. [2014-02-02]. http://minisat.se/downloads/minisat-2.2.0.tar.gz.

編 輯 葉 芳

SAT-Based Algorithm for Finding Cycles in a Boolean Network

GUO Wen-sheng1, YANG Guo-wu1, LI Xiao-yu1, and GAO Min2

(1. School of Computer Science and Technology, University of Electronic Science and Technology of China Chengdu 610054; 2. School of Electronic Engineering, University of California-Los Angeles Los Angeles Califonia USA 90034)

A cycle which is a steady state in Boolean network state transition is important in modeling checking, genetic regulatory networks and reachability analysis. Obtaining all the cycles in Boolean network state transition is a NP complete problem. In this paper, we proposes an algorithm for solving all the cycles which have less than or equal to the assigned step length based on all-solution Boolean satisfiability (SAT) algorithm. By extending the state transition function of a Boolean network and the property of cycles, this algorithm creates a problem set in conjunctive normal form (CNF) and incorporates techniques such as conflict-driven clause learning (CDCL), non-chronological backtracking, blocking clause and variable classification to decrease computational complexity. Experiment result shows that this algorithm is efficient for solving the cycles. For large scale complex network that is hard to get all the cycles, this algorithm delivers more practical value.

blocking clause; boolean networks; boolean satisfiability; cycles

TP301.6

A

10.3969/j.issn.1001-0548.2015.06.015

2014 ? 05 ? 09；

2014 ? 09 ? 12

國家自然科學基金(61272175)

郭文生(1976 ? )，男，博士生，副教授，主要從事形式化方法和軟件測試技術方面的研究.