基于模拟退火算法的微流控芯片结构优化设计*

石运琪，王英志，唐雁峰

（长春理工大学电子信息工程学院，吉林 长春 130022）

0 引 言

连续流微流控生物芯片被称为片上实验室（lab on chip，LOC），有着小型化、自动化、高通量、检测试剂消耗少等优点，有利于生物检测、化学实验的快速性和准确性。在应用方面，微流控芯片可以进行免疫测定、临床病理、快速测定原型、DNA提取和测序［1］等重要工作。微流控芯片的主体结构由上、下2 层片基组成聚甲基丙烯酸甲酯（polymethyl methacrylate，PMMA）、聚二甲基硅氧烷（polydimethylsiloxane，PDMS）、玻璃等材料，包括微通道、混合器件、微阀、反应器、进出口等结构，其中，微阀不仅可以在微流道交叉点处控制实验试剂的通过与停止，更决定着整体芯片的制造成本，精密的微流控芯片可能包含成百上千个微阀。

微流控芯片的设计具有极强的复杂性和易错性，传统设计几乎都是技术人员手动完成，过程繁琐且复杂，这一过程耗时耗力，极大影响了微流控芯片技术的发展，芯片布局自动化设计尤为重要。文献［2］中所提出的算法收敛速度过慢，且器件放置方法单一，过程冗余，无法适应布局条件过于庞大的设计；文献［3］中在布局布线时没有考虑流道交叉点的重要性，建模时基于基点设计没有以芯片面积为前提，降低了算法发现更有效器件放置的可能性；文献［4］中将器件放置和流道布线视为相互独立的设计，且过程复杂度过大，损害了芯片设计的质量；文献［5］采用的是普遍的启发式技术，通常以任意方式便利搜索空间，不能保证根据不同质量标准确定好的设计，并且算法只最小化了所构造的steiner树直径，可能导致芯片设计失败；微流控芯片设计中，芯片面积是主要权重之一，文献［6］中没有对其考虑，提高了芯片布局方案的缺陷；文献［7］中忽略了器件布局的关键性，在芯片自动化设计中只进行了流道布线，造成了算法整体运行时间过长的不足。

本文在考虑器件布局和流道布线的基础上，以芯片总面积、总流道长度、流道交叉点为维度，对微流控芯片进行自动化设计，当生物检测实验过于复杂时，本文通过快速序列对（fast sequence pair，Fast SP）算法进行器件初始布局，再结合改进的模拟退火（simulated annealing，SA）算法进行芯片布图设计。

1 微流控芯片工作原理

微流控芯片通过半导体的加工技术来构建整体流道系统，将生物实验和分析过程从实验室转移到由众多反应器件和流道构成的微小芯片上，加入反应试剂，通过外部硬件设施（微机械泵、电水力泵）驱动，试剂为连续流体，于芯片上形成反应，如图1所示。

图1 工作原理示意

2 微流控芯片的结构设计算法

2.1 结构设计

结构设计建模如图2所示，其流程主要包括：1）生成生化实验的流程；2）进行器件整理汇总；3）器件布局；4）流道布线；5）确定微阀数量、流道总长度、芯片总面积。

图2 微流控芯片结构设计建模

通过本文算法，将器件的集合、器件的连接关系、器件管道之间的约束作为前提条件，在一定工作时间内达到微流控芯片总面积、微流道总长度的解相对优秀的同时，减少微流道的交叉数量（减少芯片微流道、混合器件、反应器件等结构内的微阀数量），实验最终目的是降低微流控芯片设计的实验复杂度，提高芯片质量，同时降低芯片制造成本。

2.2 算法设计

本文使用Fast SP算法［8］表示器件布局的初始方案，根据前提约束条件，通过SA 算法调节器件之间间距，采用Dijkstra算法布线。在布局的初始阶段，如果直接进行微通道布线，会导致实验冗余，继而可能引起实验失败，芯片制造质量下降等问题，因此，本文算法在得到器件初始布局后即计算芯片总面积，根据前提条件中的器件连接关系进行连接，算出器件之间的曼哈顿距离，将其作为流道总长度的权值，再计算出器件连接后线段的交叉点数量作为最终布线后交叉点的权值。

2.3 Fast SP

Fast SP算法相比于传统的序列对算法降低了解码时间的复杂度，在融合了最长公共子序列（longest common subsequence，LCS）算法［9］后，加快了评估序列对布局的时间，LCS算法有巨大的实用性［10～12］。从序列对到1 个布图，给出1个序列对，需要得到：1）布图的起始点，2）每个器件的宽和高，3）布图的布局方向；才可以生成1 个布图。起始点是指开始的位置，即第一个器件放置的位置，计算器件之间相对位置时，需要器件的尺寸，布局方向能将芯片面积达到最小化。本文中，采用从左下到右上的布局方案。

从序列对生成布图和寻找2个序列中的最长加权公共子序列是相关的，每一个器件的x坐标就相当于计算LCS（S+，S-），确定每个器件的y坐标就相当于计算LCS（，S-），为S+的逆序列。

器件约束关系：序列对（S+，S-）中的每个序列单元表示每一个器件，给定2个器件a和器件b。在序列S+和序列S-中，如果a 均在b 之前，则a 在b 左边；若在序列S+中，a在b之前，在序列S-中，a在b之后，则a在b的上面。

序列对计算算法：

输入：序列对＜S+，S-＞，n个器件的宽度（长度）widths［n］（heights［n］）

输出：x（y）坐标xcoords（ycoords），布图结果的W×H尺寸

1）for（i＝1 ton）

2）weights［i］＝widths［i］∥器件宽度的权重

3）（xcoords，W）＝LCS（S+，S-，weights） ∥x坐标，总宽度W

4）for（i＝1 ton）

5）weights［i］＝heights［i］∥器件高度的权重

7）（ycoords，H）＝LCS（，S-，weights）∥y坐标，总高度H步骤1，2 对n个器件的宽度进行初始化weights。步骤3寻找S+和S-的最长加权公共子序列，计算每个器件的x坐标。步骤4，5对n个器件的高度进行初始化heights。步骤6得到S+的逆序列。步骤7 根据得到与S-的加权最长加权公共子序列算法计算每个器件的y坐标。

LCS算法：

输入：序列S1和S2，n个器件的权重weights［n］

输出：每个模块的位置positions，总长度L

1）计算2个序列S1和S2的加权LCS 的宽度，利用矢量block_order记录每个器件在S2中的索引。

2）将矢量lengths初始化为0，储存每个器件的最大量（长或宽）。

3）定义变量block为S1中的当前器件，index是当前器件在S2中的索引。所有在block左边的快都被安排到长度为lengths［block］区间内，而block随后放置。

4）更新布局的总长度lengths，lengths［n］表示n个器件都被确定后的总长度。

5）如果lengths［j］超过当前长度，则更新。

提取S1＝S+，S2＝S-，weights ＝widths 时的LCS 算法，定义布图的x坐标。提取S1＝SR+，S2＝S-，weights ＝heights时的LCS算法，定义布图的y坐标。

例子：从序列对到一个布图

输入：1）随机给出序列对S+＝＜acdbe ＞，S-＝＜cdaeb ＞；2）布局方向：自左下开始；3）布局起点（0，0）；4）器件abcde。

求解过程：widths［a b c d e］＝［8 4 4 4 4］

heights［a b c d e］＝［4 3 5 5 6］

求出x的坐标：

S1＝S+＝＜acdbe ＞，S2＝S-＝＜cdaeb ＞

weights［a b c d e］＝widths［a b c d e］＝［8 4 4 4 4］

block_order［a b c d e］＝［3 5 1 2 4］

lenghts ＝［0 0 0 0 0］

步骤（1）迭代i＝1：block ＝a

index ＝block_order［a］＝3

positions［a］＝length［index］＝lengths［3］＝0 ∥（计算当前器件位置）

t_span ＝positions［a］+weight［a］＝0 +8 ＝8 ∥（确定当前器件长度）

从index ＝3到n＝5 更新lengths矢量；lengths ＝［0 0 8 8 8］。步骤（2）～步骤（5）迭代开始，重复上述操作，x坐标：positions［a b c d e］＝［0 8 0 4 8］，求得布图的宽度W＝lengths［5］＝12

求出y坐标：

S1＝SR+＝＜ebdca ＞，S2＝S-＝＜cdaeb ＞

weights［a b c d e］＝heights［a b c d e］＝［4 3 5 5 6］

block_order［a b c d e］＝［3 5 1 2 4］

lenghts ＝［0 0 0 0 0］。步骤（1）～步骤（5）依照寻找x坐标的方式迭代5次，布图的高度H＝lengths［5］＝9

输出：布图面积大小为：W＝12 ×9。各器件左下角坐标为：a（0，5）b（8，6）c（0，0）d（4，0）e（8，0），结果如图3。

图3 器件示例布局结果

2.4 改进SA算法

本文对SA算法进行改进也是从这2个方面入手［13～16］。

1）SA在温度冷却控制方法主要分2种：快速降温方式（RSA）：T＝T0／log（1 +N）；一般降温方式（CSA）：T＝qT0+k。针对传统SA起始降温速度过快，不能得到全局最优解的缺点，本文所提出的改进SA 算法从降温速率上进行改进，降温速率函数如下

式中T0为初始温度，N为外循环迭代次数，T为当前温度。

由图4可知，在迭代过程中，本文所提出的改进SA 算法在高温阶段下降缓慢，实现了算法的全局搜索，有利于最优解的产生。

图4 3 种降温函数曲线

2）降温速度慢易造成全局收敛速度缓慢，利用Fast SP算法优先提出方案与SA 算法相结合的方法进行融合改进，加快收敛速度，对Fast SP算法方法改进如下：利用逆转变换，选择2个器件单元后，逆转这2个单元之间所有的单元；选择3个器件单元，将2个器件单元之间的单元调换到第3个单元后面。

改进SA 算法流程：利用器件的集合和器件连接关系为前提，以Fast SP算法求得的器件布局方案为输入，为达到良好布局布线结果，将器件mi 与它右侧或上方的器件mj之间的距离定义为rx和ax。将所有器件集合间距的宽和高定义为WX和HY。集合WX和HY之间元素的约束条件为［emin，emax］。WX和HY构成初始解S。继而进行SA算法环节，需要设置的主要控制参数有初始温度为T、外循环次数N、结束温度Tend、当前温度T0、以及链长L（每个T时的迭代次数）。

以上述已知条件为前提：1）初始化WX和HY，emin＜rx，rx＜emax。2）设置状态变量S＝（S+，S-，WX，HY），设置初始温度T，当初始温度大于最低温度时，迭代开始。3）调节状态变量S→S´，随机生成变量，生成和，与emin，emax进行比较，当＜emin时，令emin＝；当＞emax时，令emax＝。同理求得。4）利用Metropolis 准则df＝E（S´）-E（S），如果df＜0，则概率1接受新的布局，否则概率exp（-df／T）接受新的布局。5）利用新的降温速率函数进行降温，若T0小于结束温度，则停止迭代输出当前结果。

对SA处理器件布局进行仿真实验，如图5，实验样本来自现有算法［17］。

图5 2 种算法下器件布局设计的平均收敛速度

本文提出改进SA算法不但收敛性快，相较于传统SA算法，且在面临较多器件时，仍能得到较优的解。

本文设置参数：emin＝2，emax＝4，T＝1 200，Tend＝10-4，L＝100。评估芯片布局质量函数

式中A为布局后芯片面积，B为流道交叉点数量，C为流道线段长度，C2为总线段的平方。将α的权重值定为2，将β的权重值定位200，将γ的权重值定位40。

利用Dijkstra 算法［18，19］找出布局后的最短路径，算法的输入主要包括：1）非负边权重图G（WX，HY）。2）一个开始源节点s。3）一个目标结束节点t。引入2个集合：M（包含已求得最短路径的点和最短路径长度），N（未求出的最短路径的点和该点到源节点距离）。初始化2个集合，从N集合中找出路径最短的点，加入M集合，更新N集合中器件，循环迭代，直至遍历结束，找到微流控芯片最佳布线结果。

3 实验结果与分析

本文选取6组测试实例，实验结果表明，本文的综合算法明显优于手工布局的方法，与现有算法［17］进行实验结果时，芯片总面积平均减小19.7 %，流道总长度平均减小14.1%，流道交叉点数量平均减少25.8%。由于芯片布局问题属于非确定性多项式（non-deterministic polynomial，NP）完全问题，是多项式复杂程度的非确定性问题，在器件达到一定数量时，无法得到最优解，由表1 可看出，本文的综合算法可以在一定时间内得到较高质量的解，从而满足设计需要。

4 结 论

根据实验结果表明，本文所提出的算法可以满足微流控生物芯片的自动化设计需求，通过芯片面积的减少、流道的缩短、交叉点数量的下降，验证了本文方法在微流控芯片在结构设计上的有效性和优越性。