基于禁忌搜索的数字微流控生物芯片多目标综合优化算法

王 鹤

(河南工程学院机械工程学院 河南 郑州 451191)

0 引 言

数字微流控技术的快速发展，使得生化检验不再局限于只能通过传统实验室的方式来实现，它推动着生化检验向着微型化、集成化、自动化与便携化方向发展。样品消耗量极小、成本低、可重复使用等优点使其在诸多领域都有广泛的应用，如生化检验和医药分析等方面[1-2]。数字微流控技术是以离散的微液滴为单元，通过某种液滴驱动方式实现对液滴的产生、输运、合并、混合、分离、存储和检测等多种基本操控或处理。其中，驱动液滴完成操作的最常见方式之一就是介电湿润技术。按既定次序对电极施加电压，液滴在电湿润力的作用下在电极阵列上执行上述各种操作，以完成指定的生化检验分析[3-4]，如图1(a)所示。

(a)

(b)图1 数字微流控生物芯片

最小化生化检验完成时间通常是生化分析的目标之一，这是由于生物样本脆弱，容易失去活性，要想使其在芯片上长时间保持最佳的临床或实验室环境是相当困难的。因此，如果要保证生化检验分析结果的完整性，就要在一定的资源和时序约束下，尽可能地提高各样本或/和试剂操作的并行性，减少各液滴操作在芯片上的执行时间，以最小化生化检验的完成时间。资源约束是指数字微流控生物芯片一旦被制造出来，其尺寸就被固定，不可改变；时序约束是指在生化检测实验中，各个液滴操作之间具有一定的功能依赖关系，即各操作是具有一定的先后顺序的。由于数字微流控生物芯片与超大规模集成电路(VLSI)在设计方面有诸多相似之处，因此，美国杜克大学的krishnendu Chakrabarty教授等[5-8]将计算机辅助设计(Computer Aided Design，CAD)与数字微流控生物芯片系统相结合，采用自顶向下的分析算法来设计芯片，这种设计理念能够从整体上优化芯片结构，减少人为干预，并且提高芯片利用率。

所谓的数字微流控生物芯片的自顶向下的设计理念。简单而言，就是根据芯片使用者所提供的生化检测过程，利用序列图模型和某种算法对其进行架构级调度、几何级布局以及液滴的寻址。其中：架构级调度是指资源绑定和液滴操作的调度，几何级布局是指将与液滴操作相绑定的功能模块(如混合器、稀释器及存储单元等)在芯片上进行合理的几何位置布局，而液滴寻址是规划液滴在各个功能模块之间或储液池/废液池与功能模块之间的移动路线。这里的“功能模块”，其本质就是一种虚拟设备。通常每个功能模块都是由若干个电极组成的，具有液滴混合、稀释和存储等功能，而且具有重构性。图1(b)所示的混和/稀释器就属于一类功能模块。当同一时间内有多个液滴在片上执行操作时，任意两液滴之间都要保持至少一个电极的间距，这是为了避免在不同功能模块上执行操作的液滴之间发生意外混合。因此，在功能模块的外围往往需要设置一个能够完全包围该功能模块的隔离区，如图1(b)所示，网格区域就是一个隔离区，由若干隔离电极组成。由此可知，功能模块实际在芯片上占据的电极数量要比其自身所占数量大得多。

许多学者已从架构级调度、几何级物理布局和液滴寻址等多个方面对数字微流控生物芯片展开研究，但与液滴操作相绑定的功能模块在操作执行过程中，其位置通常是固定不变的，而本文根据实际空闲电极的数量和位置，适时改变某些功能模块在片上的位置，以提高片上电极的利用率，并结合改进的禁忌搜索算法来实现对生物芯片的架构级调度以及几何级布局。

1 片上生化检验系统建模

1.1 生化检验模型

数字微流控生化分析实验可以被看作是一组具有先后次序的液滴操作，可用如图2所示的有向无环图进行描述。

图2 多元生化分析实验有序图模型

该有向无环图模型可表示为G(V,E)，节点集V={vi:i=0,1,2,…,18}，每个节点代表一个操作，而I和M则分别代表液滴生成操作和混合操作，其中:vj(j=1,2,…,8)对应8个样本或试剂液滴生成操作I，vk(k=9,10,…,14)分别对应1个稀释操作DL1和5个混合操作M，vl(l=15,16,17)则对应3个液滴检测操作D。另外，设置了两个没有任何液滴操作的空节点NOP，即v0和v18。两个液滴操作的前后依赖关系用边集E={(vi,vj):i,j=0,1,2,…,18}来表示，即操作vj必须在vi结束之后才能开始。每个操作vi的持续时间均用各个节点的权重ωi来表示，其中两个空操作节点的权重设置为0。系统参数决定了每个液滴生成操作的时间，而液体的流动特性对其影响可忽略不计[9]。在两个小液滴完成混合操作形成大液滴之后，通常还要对该大液滴进行分离操作，这样做的目的是为了确保在生化检验中各液滴的体积保持一致。因此液滴分离操作的时间均已包含在液滴混合操作的完成时间之内，不再另行单独计算。与液滴生成和混合等操作相比，液滴输运操作时间非常短，故忽略不计，因此在有向无环图中，任意两节点之间的边权重均设置为0。液滴操作不同，其所对应的资源需求和完成时间也不同，而资源需求包含可重新配置的资源需求P和不可重新配置的固定资源需求Q。在所有液滴操作中，液滴生成和液滴检测操作对应的资源需求属于Q，而液滴输运、合并、混合、稀释或分离等操作所对应的资源需求则属于P。

1.2 生物芯片系统综合问题描述

数字微流控生物芯片的系统综合问题非常复杂，可将其分解成四个部分来分析解决：(1) 资源绑定。在功能模块库中为某个液滴操作vi选取一个相应的功能模块Mp×q，那么该操作必须在此模块上执行。(2) 操作调度。在资源绑定和各约束(包含资源约束和时序约束)的前提下，明确各个液滴操作vi的起始时间。(3) 功能模块的几何布局。在芯片上为每个液滴操作vi绑定的功能模块找到适当的物理位置。(4) 液滴寻址。规划液滴在各个功能模块之间或储液池/废液池与功能模块之间的移动路线。根据实际情况，本文对生物芯片的系统综合的目标有两个，一是提高电极利用率；二是减小生化检验完成时间。这里，液滴寻址不作为本文研究的内容，只单从资源绑定、操作调度和功能模块的物理布局这三个角度展开研究。

在传统的数字微流控生物芯片的系统综合算法中，与液滴操作相绑定的功能模块在操作执行过程中，其位置通常是固定不变的。表1给出了液滴在不同类型功能模块上完成操作所对应的时间[10-11]。

表1 不同混合器/稀释器完成操作所需时间

假设液滴生成操作时间和光学检测时间分别为2 s和20 s。另外，设定每种样本或试剂的储液池数量Nr=1，每种酶测定所需检测器的数量Nd=1，两者均属于固定资源需求Q。在确定各模块位置之后单独执行液滴寻址，因此，这里并未将其列为本文研究内容。针对图2所示的生化分析实验，图3给出了一个系统综合方案，其中M2×2混合器2个，分别绑定给操作M1和M2；M2×4稀释/混合器各1个，分别绑定给操作DL1和M4；M2×3和M1×5混合器各1个，分别绑定给操作M3和M5。根据以上资源绑定方式以及图3所示的功能模块布置于片上的位置，图3(a)给出了最佳液滴操作调度方案。在该方案中，一旦将功能模块布置在芯片的某个位置，在与该功能模块相对应的液滴操作执行的过程中，该模块的位置保持不变。这里，调度方案选择用甘特图表示，每个长方形代表一个液滴操作，其长度表示该操作的执行时间。比如，混合操作M3于t=2 s开始，t=8.1 s结束，在M2×3混合器上执行，该操作持续时间为6.1 s。受空闲电极在芯片上分布的限制，操作M4必须要在操作M3结束之后才能开始执行，而且由于要等到操作M4结束之后，与操作M3和M4所对应的液滴才能开始执行操作M5，所以操作M3完成之时即t=8.1 s时，要暂时将液滴存储在存储器S上直至操作M4完成，具体见图3(c)和(d)所示。从图3(c)可知，空闲电极的总数量大于操作M4绑定的功能模块M2×4所需要的电极数量，但由于空闲电极的分布零散，造成无法在t=4.9 s时刻在片上布置混合器M2×4的位置，降低了电极的利用率，液滴操作的并行性被削弱，同时生化检验完成时间也大大增加。针对以上问题，本文利用功能模块可动态重构的特点，适时改变功能模块的位置，以提高电极的利用率，最小化生化检验完成时间。

(a) 由甘特图表示的调度方案

(b) t=2 s

(c) t=4.9 s

(d) t=11.95 s

(e) t=14.85 s图3 生化检验系统综合实例

2 功能模块动态移位

数字微流控生物芯片上的每个功能模块都是一个虚拟设备，不同电极组合得电可使得该功能模块在芯片上的位置发生变化，因此，功能模块具备可动态重构这一特性。利用这一特点，在液滴操作执行过程中，本文适时改变某功能模块位置，以最终达到提高片上电极的利用率和生化检验完成时间最小化这两个目标。

图4给出了基于以上思想的与图2生化分析所对应的新的系统综合方案。虽然图4(b)所示的空闲电极的总数大于操作M4绑定的功能模块M2×4所需要的电极数，但空闲电极的分布零散使得在t=4.9 s时刻无法为混合器M2×4找到合适的位置。因此，将混合操作M3绑定的功能模块M2×2在t=4.9 s时刻，整体向左移动三个纵列的位移，即可将空闲电极聚集，随之可将M4的功能模块M2×4布置在芯片上，如图4(c)所示。混合器M2×2在操作M3的执行过程中，产生这样的位置移动，大大提高了液滴操作的并行性，使得生化分析的完成时间由38.05 s(图3(a))减少至34.85 s(图4(a))，减少了8.4%。当然，一旦混合器M2×2的位置发生变化，就需要对其上的液滴在新旧两个位置之间进行了路径规划，但由于液滴输运操作的时间极短，因此该液滴在混合器M2×2的两个位置之间的移动时间可忽略不计。除此以外，图3(d)中包含的存储器S在图4所示的新方案中也可省略不用，减少了功能模块的使用数量。

(a) 由甘特图表示的调度方案

(b) t=2 s

(c) t=4.9 s

(d) t=8.1 s图4 改进的生化分析系统综合方案

下面采用改进的紧急搜索算法来实现以上的数字微流控生物芯片的系统综合以及功能模块的动态移位。

3 改进的禁忌搜索算法实现系统综合

禁忌搜索(Tabu search，TS)算法是一种亚启发式的随机搜索算法，是对局部领域搜索的一种扩展，是一种全局迭代寻优算法，是对人类记忆功能的一种模拟。TS算法的主要思想就是从一个初始可行解开始，对该解的邻域进行搜索并优选，同时设置禁忌表，用来存储已经搜索过的局部最优点，以便在后续搜索中有意避开(但并不是绝对禁止)禁忌表中所包含的局部最优点。也就是说，在搜索过程中，当最优候选解在禁忌表中时，则选择未被禁忌的次优候选解来替代当前解，因此该算法是通过引入一个灵活的存储结构和相应的禁忌准则来避免迂回搜索的。另外，当候选解优于最优解，但该最优解却被禁忌时，利用特赦准则来赦免这些最优解，从而保证探索的多样化以最终达到全局最优的目的[12-13]。但是通过算法的实际应用发现，初始解的好坏会对算法的搜索结果产生很大的影响，因此这里引入以长期记忆为基础的变异因子对传统的禁忌搜索算法进行优化，以此完成对数字微流控生物芯片的系统综合以及功能模块的动态移位。

在本文提出的系统综合算法中，以有向无环图模型G(V,E)、m×n生物芯片电极阵列、功能模块库L等作为算法输入，完成资源绑定、操作调度和功能模块的物理布局，实现提高片上电极利用率和生化检验完成时间最小化这两个目标。对于由TS确定的资源绑定，我们使用列表调度来确定操作vi的调度。列表调度是以优先级列表PL为基础，其中包含了所有已准备好被调度的操作，这样操作的所有前驱操作均已完成。每次迭代，优先级最高的操作vi被调度，并且该操作被调度之前，其绑定的功能模块Mp×q已被放置在芯片上，而各操作的优先级也是由TS确定的。一旦确定该操作的功能模块Mp×q在片上的位置，列表调度就会根据Mp×q与源模块之间的曼哈顿距离计算液滴寻址时间。TS将每个操作vi绑定到一个随机选择的功能模块Mp×q，由此产生算法的初始解，并且利用关键路径优先级函数给出各操作的优先级。

3.1 功能模块动态移位算法实现

以图3(c)为例，我们给出功能模块动态移位算法的具体实现过程。在t=4.9 s时刻，操作DL1已结束，这时优先级列表中包含了三个已准备好被调度的操作，即M1、M2和M4，三者具有由高到低的优先级排序，而且前两者均绑定到模块M2×2，后者绑定到模块M2×4。列表调度在t=4.9 s时刻首先为操作M1的功能模块M2×2在片上找到合适的物理位置。选择物理位置的基本原则就是尽量将功能模块集中放置在芯片上，这样的话，剩余空闲电极分布也会比较集中，便于并行布置其他功能模块，提高操作的并行性。在t=4.9 s时刻，操作M3的功能模块将空闲电极分成了四个交叠的矩形区域a1、b1、c1和d1，如图5(a)所示。这四个区域分别用矩形的左下角坐标和右上角坐标表示，即(0,0,3,8)、(3,8,8,11)、(7,3,11,8)和(0,0,8,4)。除区域a1以外，其他三个区域均可容纳下操作M1的功能模块M2×2。假设该功能模块被放置到区域d的左下角位置，如图5(b)所示。之后更新空闲电极分布区域为a2、b2、c2和d2，即(0,4,3,8)、(3,8,8,11)、(7,3,11,8)和(4,0,8,4)。同样除区域a2以外，其他三个区域均可容纳下操作M2的功能模块M2×2，但是区域d是正好能容纳下M2，因此，基于为功能模块选择物理位置的基本原则，将M2放置在区域d，这样不会造成剩余空闲电极分布过于分散，以增大剩余空闲电极容纳下操作M4的几率，如图5(c)所示。但是，我们发现，即使这样选择M2功能模块的物理位置，虽然剩余空闲电极总数大于操作M4的功能模块M2×4所需要的电极数量，但是由于空闲电极分布分散，实际上是无法为操作M4的功能模块M2×4找到合适的物理位置的。因此，本文在操作执行过程中移动某功能模块位置，尽可能减小空闲电极的分散度。

(a) t=4.9 s时刻功能模块初始状态

(b) 放置M1绑定的功能模块

(c) 放置M2绑定的功能模块

(d) M3绑定的功能模块动态移位

(e) 放置M4绑定的功能模块图5 功能模块动态移位实现过程

这里采用贪心搜索策略来决定哪个功能模块需要移动，而且对正在执行某一操作的功能模块进行位置移动的话，还需要规划其上的液滴从模块初始位置到另一位置的移动路径。模块位置移动的路径相当于其上液滴的移动路径，路径距离根据该模块前后两个位置左下角之间的曼哈顿距离计算。我们对这种情况下的液滴寻址时间设置一个约束，即一个时间步长。因此，一旦能够容纳下即将被调度的操作绑定的功能模块或达到上述时间约束，就停止对功能模块的动态移位。当然，如果动态移位之后发现没有足够的空间可以容纳即将被调度的操作功能模块，那么恢复被移动功能模块的初始位置。在每次迭代中，贪心搜索策略选择最佳移动，所谓的最佳移动就是在功能模块移位之后，原来分散的两个矩形区域之间的曼哈顿距离可以达到最小，而且通过矩形区域的合并增大相邻空闲电极的数量。图5(c)显示了片上功能模块所有可能的移动，即操作M3的模块整体向左移动三个电极，或者分别向右或向上移动四个电极。贪心法选择操作M3的模块整体向左移动三个电极，也就是矩形区域a3同时向右移动三个电极，两者交换位置，这是最佳移动，因为经过这样的移动之后，图5(c)中的矩形区域a3和c3之间的曼哈顿距离只有3，并且两者合并为c4，(4,4,11,8)，包含28个相邻的空闲电极，足以容纳下操作M4的功能模块M2×4，如图5(d)、图5(e)所示，算法终止。

3.2 改进的禁忌搜索算法

算法的主要思想为：依照传统禁忌搜索算法进行搜索，获得局部最优解，将其作为父代进行变异，得到迄今为止的最优解，并将其作为下一次迭代的初始解；如此循环迭代可实现对数字微流控生物芯片的系统综合优化。

(1) 邻域 紧急搜索算法是一种基于邻域搜索技术的元启发式算法，通过对当前解xnow进行设计转换在其邻域N(xnow)生成一组邻域解，由此选出候选解x。目前广泛使用的邻域结构有交换和插入两种，根据数字微流控生物芯片系统的自身特点，本文选择交换邻域结构，对当前解主要进行两种移动，一是资源重新绑定，即随机选择液滴操作vi，解除与其当前绑定的功能模块Mp1×q1，同时将该操作绑定到另一个功能模块Mp2×q2；二是操作优先权的交换，即随机选择两个液滴操作，交换两者的优先权。

(2) 禁忌表 禁忌表是禁忌搜索算法的关键，存储内容包括禁忌对象和禁忌长度两项，决定了禁忌表中各对象的任期以及搜索范围。本文使用两个禁忌表，每种移动都对应有一个禁忌表。如果随机选择操作vi被资源重新绑定到功能模块Mp2×q2，那么它们就会以(vi,Mp2×q2)的形式列入到与资源重新绑定对应的那个禁忌表中；如果随机选择操作vi和vj交换优先权，那么这两个操作以(vi,vj)的形式被列入与操作优先权交换对应的禁忌表中。因此，被列入禁忌表中的解在禁忌长度内不允许被再次访问，但为了防止算法中的所有对象都被禁忌这种情况出现，也为了避免遗失优良状态，使用“特设准则”来加速实现全局优化。

(3) 特赦准则 所谓的特赦准则，是指如果邻域候选解集中存在一个新解，优于当前最优解，即使这个新解是被禁忌表禁止的，算法也会选择这个更优的新解，即这个新解被解禁。

以图2生化分析中各操作资源重新绑定移动为例，图6给出了禁忌搜索的三种邻域解。液滴生成和液滴检测操作属于不可重新配置的固定资源需求Q，所以不在资源重新绑定的范围内。对图6(b)-(d)的三个邻域解分析发现，虽然图6(b)对应的解是被禁忌表禁止的，但这个解优于当前最优解，所以算法会解禁该解，并将其替代当前最优解。

(a) 当前最优解及禁忌表

(b) 操作M3的功能模块由M2×2重新绑定到M2×4，这属于禁忌，但优于当前最优解

(c) M4的功能模块由M2×4重新绑定到M2×2，这属于非禁忌移动，而且优于当前最优解

(d) M5的功能模块由M1×5重新绑定到M1×4，这属于非禁忌移动，但该解比当前最优解差图6 禁忌搜索邻域(资源重新绑定)

(4) 变异因子 为进一步扩大算法的搜索范围，参考遗传算法，将当前局部最优解作为父代进行变异。这里的变异就是将液滴操作的资源绑定和优先权分别进行交换。由于当前解在邻域内的优先权交换就是在两个液滴操作之间进行的，如果在变异阶段仍然采取两个操作之间交换进行变异，那么扩大搜索范围的能力会受到影响。因此，这里我们选择对两个以上的液滴操作进行资源绑定和优先权进行交换来完成变异。下面以图2所示的生化分析为例，给出变异的主要思想。就资源绑定变异来讲，由于液滴生成和液滴检测操作属于Q，所以只有DL1、M1至M5这六个操作的资源绑定可以进行变异：

① 从上述6种液滴操作中随机选择3个(可依具体问题确定随机选择的数目)，如表2所示，随机选择M1、M3和M4；

表2 变异算子

② 3个液滴操作的资源绑定完全交换，生成两个子代；

③ 在两个子代中找出最优解，如果变异后的子代最优解优于父代，那么将该子代最优解替代当前迭代最优解，并且将其作为下一次禁忌搜索迭代中液滴操作资源绑定的初始解。

因此，改进的禁忌搜索算法实施步骤如下：

(1) 初始化各个控制参数；

(2) 随机生成初始解，禁忌表设置为空集,并计算初始解所占用的电极数Ne以及完成生化检验的总时间T；

(3) 对当前解xnow进行设计转换在其邻域N(xnow)生成一组邻域解，并选出候选解xi，计算该解所占用的电极数Nei以及完成生化检验的总时间Ti；

(4) 判断xi是否满足特赦准则且Nei≥Neinow且Ti≤Tinow，如果均满足，则更新全局最优状态，最优解x0=xi，并且更新禁忌表；

(6) 判断是否达到设定的迭代步数，若是，则结束算法，并且输出全局最优解,反之，将算法转至(3)，重复上述步骤。

4 仿真结果及分析

为了验证算法的有效性和可靠性，本文通过对人体体液体外诊断分析实验的仿真，将其他算法与该算法进行了相应的比较。人体体液体外诊断实验可参照文献[14]，共有28个液滴操作，其中包含13个液滴生成操作、6个混合操作、3个稀释操作和6个液滴检测操作。

首先分别采用0-1整数线性规划算法(ILP)[15-16]与本文提出的禁忌搜索算法，实验仿真在三个不同面积的芯片上实施，两种方法对应的计算机仿真时间以及其最优解所对应的生化分析完成时间的对比如表3所示。可以看出，利用本文提出的禁忌搜索算法在三种不同面积的芯片上均可在5 min内完成计算仿真，时间远远低于采用ILP算法的仿真时间。

表3 ILP与TS的比较

其次，定义去除功能模块动态移位的禁忌搜索算法为TS*，用本文提出的改进的禁忌搜索算法TS与TS*分别在3 min和5 min仿真时限内，在三种不同面积的芯片上执行人体体液体外诊断实验，两种算法在每种情况下均运行50次，对比两者的仿真结果，如表4所示。

表4 TS与TS*的比较

由表4可知，芯片尺寸越大，操作并行性越高，生化实验完成时间越短；在一定的时间范围内，仿真时间越长，生化实验完成时间也越短。但是从算法50次仿真得到的最优解的平均值来看，在相同的仿真时限内，芯片尺寸越小，功能模块的动态移位这一因素对增大电极利用率、提高液滴操作的并行处理和减小生化实验完成时间的影响能力越高，比如，在3 min时限内，用TS在10×10、11×11和12×12芯片完成生化实验的时间分别比TS*减少了10.58%、7.12%和4.19%。

5 结 语

本文利用功能模块的动态重构特性，在某个操作执行的过程中，适时改变其绑定的功能模块在片上的位置，增大液滴操作的并行处理，同时结合改进的禁忌搜索算法来实现功能模块的动态移位以及数字微流控生化检验的架构级调度和几何级布局，以达到提高电极利用率，最小化生化检验完成时间的目的。通过人体体液体外诊断实验的仿真，对多个算法进行比较，仿真结果验证了本文算法的有效性和可行性。而且该算法同样也可用于其他生化实验的实施，对数字微流控生化检验的系统综合具有一定的参考价值。