0 引言

约束可满足问题(constraint satisfiability problem, CSP)是人工智能科学领域中一个重要的研究方向，许多实际的优化问题都可以建模转化为CSP，这使得CSP通过运用不同编码方式的可满足性判定问题(satisfiability problem，SAT)更加拟合现实问题中的约束关系。CSP在电路设计^[1]、调度^[2]与规划^[3]等领域有着广泛的应用，因此对于CSP的研究探索就显得极其有意义。Xu等^[4]提出了一种新的随机CSP模型，名为RB模型，该模型是一种具有精确相变的增长域实例模型。此后，有学者在RB模型的基础上提出了多个改进模型，如K-CSP^[5]、d-K-CSP^[6]与p-RB^[7]等模型。这几种模型都具有一定的局限性。研究人员还分别在局部搜索^[8]、元启发式^[9-11]和完备算法^[12-14]中对RB模型实例的求解进行了研究。由RB模型生成随机实例的求解难度，在变量规模仅为10²量级时就变得极有挑战性^[15]。

元启发式算法通常被用于处理连续型问题，而RB模型实例却是典型的离散型问题。本文以海洋捕食者算法(marine predators algorithm, MPA)^[16]为基础，提出了一种离散化的海洋捕食者算法(dispersed marine predators algorithm，DMPA)，并将其与退火策略中的最小冲突启发式算法(simulated annealing with min-conflicts heuristic，SMCH)相融合求解RB模型实例。实验结果表明，相较于文献[9]中改进的模拟退火算法(revised simulated annealing algorithm，RSA) 和退火遗传算法(genetic simulated annealing algorithm，GSA)，求解概率与求解精度均有显著提升，算法的时间效率最大提升29%。在各组变量规模逼近可满足相变点时，DMPA-SMCH算法与其他启发式算法对比依然有高概率解。

1 相关知识 1.1 RB模型

一个RB模型实例是由n个彼此互不相同的变量集X={x₁, x₂, …, x_n}组成，每个变量x_i都由对应的值域D={d₁, d₂, …, d_n}和一个有限制的约束C={C₁, C₂, …, C_m}组成。RB模型的一类随机CSP实例表示为I∈RB(k, n, α, r, p)，其中对于每个实例的符号定义^[4]如表 1所示。

生成一个实例I∈RB(k, n, α, r, p)的具体步骤如下。

1) 设置模型RB的主要参数n，p，α和r。其中n是变量数量，p是约束密度，r和α是两个自定义的常数。

2) 重复随机选择m个约束条件，每个约束都是通过从n个变量中独立随机地选择k个变量形成。

3) 对于每个约束C，从完整的d_n^k个赋值中随机选取pd_n^k种不同k元赋值，形成一个不相容赋值的集合Q_a，当约束C的变量取值不属于这些不相容赋值集合Q_a时，该约束满足。

给出控制参数p，临界值p_cr，可确定RB模型的可满足性相变点，并给出定理1。

定理1   设p_cr= $ 1-\mathrm{e}^{-\frac{\alpha}{r}}$，若α> $ \frac{1}{k}$，r>0为两个常数，且k、α和r满足不等式$ k \mathrm{e}^{-\frac{\alpha}{r}}$≥1，则RB模型有解的概率p(Sat)可表示为

$ \lim _{n \rightarrow \infty} p(\text { Sat })=\left\{\begin{array}{lc} 1, & p<p_{c r} \\ 0, & p>p_{c r} \end{array}\right. $

(1)

利用上述定理公式，根据RB模型的参数取值可确定相变点p_cr。更准确地说，变量数接近无穷大时，若p < p_cr，RB模型生成的随机CSP实例以高概率接近于1可解；若p>p_cr，则CSP实例高概率接近0无解。对于一个随机实例I(k, n, α, r, p)，取k=2，n=6，α=0.8，r=1，p=0.1，可根据生成的RB模型实例画出约束网络图，如图 1所示。

图 1中，方框表示约束，圆圈表示变量，约束集C={C₁, C₂, …, C₁₁}，约束个数m=11, 变量X={x₁, x₂, …, x₆}，根据以上RB模型参数可计算出值域个数为d=4，值域D={0, 1, 2, 3}，从图中可以看到每个约束C_a与两个变量x_a相关。Q_a是由RB模型实例的生成步骤C产生的不相容赋值集合，分别为：

$ \begin{aligned} & Q_1=\{(2, 2), (3, 2)\} ; Q_2=\{(0, 2), (2, 3)\} ; \\ & Q_3=\{(0, 2), (2, 3)\} ; Q_4=\{(0, 2), (2, 3)\} ; \\ & Q_5=\{(0, 2), (2, 3)\} ; Q_6=\{(1, 1), (0, 1)\} ; \\ & Q_7=\{(2, 1), (3, 2)\} ; Q_8=\{(3, 2), (0, 0)\} ; \\ & Q_9=\{(2, 1), (3, 0)\} ; Q_{10}=\{(2, 0), (2, 2)\} ; \\ & \quad Q_{11}=\{(1, 3), (1, 2)\} 。\end{aligned} $

当X={1, 3, 1, 0, 1, 1}时，为上述RB模型随机实例的一组可行解。

1.2 最小冲突启发式

最小冲突启发式算法(minimum conflict heuristic algorithm，MCH)^[17]是求解CSP的一类局部搜索算法，通常迭代地修改单个变量的赋值，使违反约束的数量最小化，MCH算法求解RB模型实例的伪代码如算法1所示。

算法1   MCH算法

输入：RB模型实例I(C, Q, d, n)，最大迭代次数m。

输出：若找到解s，输出s，否则输出违反约束数V与对应的最优赋值。

1) 初始化一组随机赋值解s；

2) 若s满足所有的约束C, 返回解s；

3) 若不满足，找出不满足的约束冲突集Q_c，且从冲突集Q_c中随机选取一个冲突变量；

4) 为选定的冲突变量在域中重新赋值，使得在当前分配下，冲突数量最小；

5) 判断赋值后的s，若不满足，重复步骤3)和4)，直到循环结束，返回违反约束数V与对应的最优赋值。

2 MPA算法的离散化 2.1 定义目标函数

对于实例I，设X=(x₁, x₂, …, x_n)是一组完全赋值，对此，将目标函数F定义为给定的潜在解决方案违反的约束个数，具体目标函数为

$ F=\sum\limits_{i=1}^m {Check}\left(X, C_m\right), $

(2)

$ {Check}\left(X, C_m\right)= \begin{cases}1, & \text { 如果 }\left(x_1^i, x_2^i, \cdots, x_n^i\right) \in C_m, \\ 0, & \text { 如果 }\left(x_1^i, x_2^i, \cdots, x_n^i\right) \notin C_m, \end{cases} $

(3)

式中：m为约束个数。检查一组赋值X是否满足约束集C_m，目标是最小化求解过程中违反约束的数量，若一组赋值解X使得目标函数值F为0，则这组赋值为实例的解。

2.2 初始解空间离散

以图 1中RB模型实例为例，其中变量数n为6，值域个数d为4，值域D={0，1，2，3}，随机初始解如图 2所示。

P为一组随机初始解，令变量x_i取值为值域D中的随机值，每个变量x_i的取值都可重复地在对应的值域中取值。算法初始化时，每个变量x_i随机选取的值构成离散初始解。当变量数为n时，变量解空间大小为dⁿ，当种群数为q时，将在解空间内初始q组随机解。计算每组解的目标函数值，并更新顶级捕食者，当算法结束时，顶级捕食者为当前的最优解。

2.3 DMPA优化阶段

由于RB模型实例解的性质，MPA对每组变量进行随机性扰动并不适合求解，对于求解过程中若有部分解已满足个别约束，须保留该部分解，但同时又希望对其他变量进行扰动。所以在MPA的第一阶段基础上进行改进，对当前最优解P_b形成的精英矩阵E按照维度分为两部分，分别进行布朗随机扰动，前半部分扰动数学模型为

$ \begin{array}{l} \left\{ \begin{array}{l} \overrightarrow {{{\boldsymbol{Z}}_{i, j}}} = \overrightarrow {{{({\boldsymbol{A}} \oplus {\boldsymbol{B}})}_{i, j}}} \otimes \left( {\overrightarrow {{{\boldsymbol{E}}_{i, j}}} - \overrightarrow {{{({\boldsymbol{A}} \oplus {\boldsymbol{B}})}_{i, j}}} \otimes \overrightarrow {{{\boldsymbol{P}}_{i, j}}} } \right), \\ \overrightarrow {{{\boldsymbol{P}}_{i, j}}} = \left\lceil {\overrightarrow {{{\boldsymbol{P}}_{i, j}}} + U \cdot \overrightarrow {\boldsymbol{R}} \otimes \overrightarrow {{{\boldsymbol{Z}}_{i, j}}} } \right\rceil , \end{array} \right.\\ \;\;\;\;\;\;\;\;\;\;\;i = 1, 2, \cdots , n;j = 1, 2, \cdots , n/2, \end{array} $

(4)

后半部分扰动为

$ \begin{array}{l} \left\{ \begin{array}{l} \overrightarrow {{{\boldsymbol{Z}}_{i, j}}} = \overrightarrow {{{({\boldsymbol{B}} \oplus {\boldsymbol{A}})}_{i, j}}} \otimes \left( {\overrightarrow {{{\boldsymbol{E}}_{i, j}}} - \overrightarrow {{{({\boldsymbol{B}} \oplus {\boldsymbol{A}})}_{i, j}}} \otimes \overrightarrow {{{\boldsymbol{P}}_{i, j}}} } \right), \\ \overrightarrow {{{\boldsymbol{P}}_{i, j}}} = \left\lceil {\overrightarrow {{{\boldsymbol{P}}_{i, j}}} + U \cdot \overrightarrow {\boldsymbol{R}} \otimes \overrightarrow {{{\boldsymbol{Z}}_{i, j}}} } \right\rceil , \end{array} \right.\\ \;\;\;\;\;\;\;\;\;\;\;i = 1, 2, \cdots , n;j = 1, 2, \cdots , n/2, \end{array} $

(5)

其中：Z_{i, j}为计算的个体步长；S为移动步长；A为布朗随机游走矩阵，维度为n/2；B为空矩阵；⊕为逐项相加符号；⊗为逐项乘法符号；E为精英矩阵；P为猎物矩阵；U在上式中设为0.5；R为[0, 1]中均匀随机数的向量。初始时，海洋捕食者数量为n，维度为d，i∈[0, q], j∈[0, n]。由于实例的离散性质，故在每次扰动之后需要对其进行向上取整，计算个体适应度，并更新当前最优解P_b形成的精英矩阵。

第二阶段，逐渐从勘探向开发策略迈进，为处理离散型RB模型实例问题，以汉明距离重新定义步长因子Z为S，具体数学模型为

$ S=\sum\limits_{i=1}^n {Check}\left(\boldsymbol{P}_{b(j)}, \boldsymbol{P}_{(i, j)}\right), $

(6)

$ {Check}\left(\boldsymbol{P}_{b(j)}, \boldsymbol{P}_{(i, j)}\right)= \begin{cases}1, & \text { 如果 } \boldsymbol{P}_{b(j)}=\boldsymbol{P}_{(i, j)}, \\ 0, & \text { 如果 } \boldsymbol{P}_{b(j)} \neq \boldsymbol{P}_{(i, j)} \text { 。}\end{cases} $

(7)

假设给定两个解，分别为当前的最优解P_b与猎物矩阵中的解P，检查两组解中的每个变量是否相同，如图 3所示。

图 3中，P_b违反约束C₂，C₁₀，根据公式(2)，可得F(P_b)为2。P违反约束C₅，C₇，C₁₀，可得F(P)为3。根据步长因子公式(6)可得S为4。由原MPA算法可知CF控制算法的收敛速度，重新定义自适应参数CF，具体为

$ C F=\exp \left[\left(-s * S^a\right) / q^a\right], $

(8)

$ { Movement } \sum\limits_{i=1}^q\left(\boldsymbol{P}_b, \boldsymbol{P}\right)=\left(\boldsymbol{P}_{b(j)} \leftarrow \boldsymbol{P}_{(i, j)}\right), $

(9)

如果P_b≠ P & & rand>CF，

式(8)中：s和a是两个自定义参数，分别取0.25和2；q为种群数取20。参数s需认真定义，它影响着是否要将较优的P向当前最优解P_b移动。s分别取0.1, 0.2, 0.25，步长S取[0, 100]时，CF自适应概率如图 4所示。

在离散解空间的猎物矩阵中，P向当前最优解P_b移动，重新定义为公式(9)。rand取[0, 1]区间内随机数。从猎物矩阵P中依次检查每组解，并与当前最优解判断是否移动，当S值增大时，预示着猎物矩阵P与当前最优解P_b差距越大，图 4中对应的CF值越小，需要取rand>CF，以更大的概率去替换个别变量。相反，当CF值越大时，预示着P与当前最优解P_b更接近，反之需以小概率去替换个别变量。由上述公式(8)与图 4中各参数s概率图，再和公式(9)经实验共同得出：当s取0.25时，更适合每组变量规模的概率传递。若移动后的解优于当前最优解，则更新当前最优解P_b，直到循环结束。在离散空间中，猎物矩阵P向当前最优解P_b移动。如图 5所示。

在图 5的步骤中，首先检查P中变量与当前最优解P_b变量是否相同，若相同，则认为该变量为有效变量，保持不变；否则，其他变量为无效变量，以自适应参数CF概率来决定要替换的无效变量。假设，图 5中P违反约束C₅，C₇，C₁₀，当前最优解P_b违反约束C₄，CF值为0.25，以0.25的概率正好替换为P中的x₅，最后更新得到P_b={1, 3, 0, 2, 2, 0}，满足所有约束，即为找到的可行解。

第三阶段采用局部搜索策略，由于在大规模实例集中，MCH收敛性效果较差，故而在DMPA算法迭代后可融入快速模拟退火策略^[9]的SMCH算法再进行寻优。以1-(3/T)的概率进行选择扰动，利用最小约束数和当前退火温度T以exp[-(F_n-F_b)/T]概率来接受新赋值。具体SMCH算法伪代码如算法2。

算法2   SMCH算法

输入：RB模型实例I，当前最优解P_b，退火参数T, T₀, a, L。

输出：若找到解s，输出s，否则输出最小违反约束数V。

1) 初始解赋值s_n= P_b；

2) 主循环内，若rand < (3/T)，将最优解赋值s=s_n；

3) 否则将P_b传递至算法1，利用F(s)计算目标函数值F并判断；

4) 若F_n < F_b，更新目标函数F=F_n; 若(rand < (exp[-(F_n-F_b)/T]))，则接受新解s，更新最小违反约束数V=F；

5) 退火冷却表设置为T_i+1=a ·T_i；

6) 重复上述操作，若F=0，输出s；否则输出V。

2.4 变异算子

对于算法随机化的处理，文献[12]中离散化的萤火虫算法求解RB模型提出扰乱变异算子(scramble mutation, ScM)，具体是在更新的最优解P_b中随机取一个子集进行扰乱。并以0.2的概率在随机值域D中取值，具体如图 6所示。DMPA-SMCH算法伪代码如算法3。

算法3   DMPA-SMCH算法

输入：RB模型实例I，算法参数种群q, 影响概率H，最大迭代次数M。

输出：若找到解P_b，输出P_b，否则输出最小违反约束数V。

1) 初始化一组随机解，找出顶级捕食者；

2) 构建猎物矩阵P，确定P的上下界，赋值当前违反约束F_n=F，赋值最优解P_b= P，用公式(2)计算适应度值F；

3) 第一阶段优化：由P_b构建精英矩阵E，用公式(4)和(5)更新精英矩阵；

4) 第二阶段优化：若F(P) < F(P_b)，用公式(6)和(7)计算S=(P_b, P)，用公式(8)计算自适应因子CF，用公式(9)更新精英矩阵；

5) 第三阶段优化：执行算法1(实例I，最优解P_b，种群数q)；

6) 突变算子：若rand>H，则以ScM方式进行变异，边界处理，计算适应度，更新最小违反约束数；

7) 若V=0，输出违反约束数V；否则执行算法2，输出最优解P_b和违反约束数V。

3 算法时间复杂度分析

DMPA-SMCH算法的主要时间复杂度由DMPA和SMCH算法两部分组成，DMPA算法主要为目标函数计算适应度值的时间复杂度。即目标函数的时间复杂度为

$ T_1(m, t, k)=O(m * t * k), $

(10)

其中：$ m=r \cdot n \cdot \ln n ; t=p d_n^k, d_n=n^\alpha$；RB模型为二元约束，k为常数，取2；α取0.8；r和约束密度p为常数级。

MCH算法的时间复杂度可以表示为

$ T_2(n)=O(f(n, M, m)), $

(11)

其中：n为问题规模；M表示迭代次数；m表示局部搜索的复杂度主要为T₁。

SMCH算法在MCH算法上融入了退火策略，算法的主要时间复杂度为

$ T_3\left(T_0, T, L, T_1\right)=O\left(\log \left(T_0 / T\right) * L * T_1 * k\right), $

(12)

则DMPA-SMCH的主要时间复杂度为

$ T=T_1+T_2+T_3, $

整理后可得T=O(n*ln n)，即当确定了RB模型的参数后，整体的算法时间复杂度为

$ O(n * \ln n) \text { 。} $

4 实验分析 4.1 实验环境与参数设置

实验环境为Windows 11, 64B操作系统，CPU为AMD Ryzen 7 5800H，主频为3.20 GHz，内存16 GB。在DMPA-SMCH算法中，迭代次数设置为800次，种群数量q设置为20，突变概率H设置为0.2，第一阶段为50次迭代，第二阶段为350次迭代，第三阶段为400次迭代。SMCH算法中初始温度T₀为97，终止温度T为3，控制温度参数α为0.85。RB模型参数与文献[9]中保持一致，如表 2所示。计算可满足性相变点为p_cr=0.234。取n=[20∶20∶100]，取约束密度p=[0∶0.01∶0.2]，进行10组实验，数据取均值对比。

4.2 算法收敛概率分析

图 7(a)、(b)、(c)分别为RSA、GSA、DMPA-SMCH求解RB模型实例各组上的各个约束密度点的收敛概率。其中RSA在p≤0.16时、GSA在p≤0.16时两种算法在各组RB模型变量实例上以概率1收敛，DMPA-SMCH在p≤0.16时可有效收敛。RSA在p≤0.16时算法失效，GSA在除变量为20时、其他实例集上，p≤0.16时算法失效。而DMPA-SMCH在除变量n取80，100时，其余各组变量收敛概率有明显提升。随着约束密度增加，可以看到RSA算法收敛能力最差，GSA算法在n为20，接近可满足性相变点依然有概率求解，但当n增加到40时，算法在实例的高维区域已无法收敛。而DMPA-SMCH算法在n分别为20，40，60时，依然在p≥0.16时有概率收敛，在高维区域n分别为80、100时，对比RSA和GSA依然保持优秀的求解能力。

图 8(a)、(b)为DMPA-SMCH在n分别为40，60，约束密度为0.08~0.18中取随机值时的求解概率。其中完备性求解算法为回溯算法^[12](backtracking algorithm，BT)、基于度启发式的回溯算法^[12](heuristic backtracking algorithm，HBT)、动态度的回溯算法^[13](dynamic degree-maintaining arc consistency，ddeg-MAC)，而RSA、GSA都为启发式随机搜索算法。基础BT算法在高维难解实例中并不具有求解能力。可以观察到图 8(a)中，当n取40时，HBT与DMPA-SMCH在同组算法中，p为0.13~0.17时展现了高概率求解，均优于其余对比算法。在p>0.15时，HBT的求解性能稍优于DMPA-SMCH算法。而在p>0.18时，由于实例的复杂度增加，各算法均不能高效求解。当n取60，p=0.12时，如图 8(b)，DMPA-SMCH算法为对比算法中最优算法。在p>0.16时，HBT求解概率略高于DMPA-SMCH算法，但DMPA-SMCH算法与对比算法相比，将概率为1的解提升至p=0.13，其为所有对比算法中的最佳值。

4.3 求解效率对比分析

图 9(a)，(b)，(c)，(d)，(e)为三种算法在n分别取值为20，40，60，80，100时，每10组实验的平均时间对比图。由于启发式算法有一定的随机性，可以看出，除约束密度p < 0.11时，DMPA-SMCH算法的各组变量、各个约束密度求解的时间都要优于RSA和GSA算法。在p < 0.11时，DMPA-SMCH算法在各组变量低密度时，求解时间效率都劣于RSA与GSA算法。但当p>0.11时，可以看到第二阶段在难解的大规模实例上，可以有效地优化大量约束条件，在其他各组变量上都展现出了优秀的求解效率。图 9(f)是三种算法在每组变量约束密度为0~0.20范围内，运行10次的平均总时间对比图，在n分别取20，40，80时，DMPA-SMCH算法的总时间都优于RSA算法与GSA算法。虽然在n分别取60，80时，DMPA-SMCH算法在总平均时间上略差于RSA算法和GSA算法，但是DMPA-SMCH算法保持了优于RSA算法与GSA算法的高效求解精度。

表 3为三种算法在随机的RB模型实例上各组变量运行10次，取约束密度为0.16~0.20范围内接近可满足性相变的实验结果。可以看到，DMPA-SMCH算法无论在哪种实例集上，平均运行时间都要优于RSA算法与GSA算法，且随着变量规模与约束密度的增大，时间优化效率越明显。精确度在绝大多数实例集上优于RSA算法与GSA算法。

4.4 DMPA-SMCH参数选择与性能分析

DMPA-SMCH算法在第一阶段时，算法搜索处于一个全局扰动状态，改进的DMPA-SMCH第一阶段在第50次迭代时便能达到原MPA算法下100~200次的迭代效果。该阶段大于50次迭代时，并没有达到有效优化，经过多次实验取50为最优迭代次数。第二阶段是全局向局部搜索的一个转变过渡过程，该阶段由自适应因子CF控制，经过猎物矩阵P中的i个解进行概率传递优化，i对应为种群数q。由于种群q和自适应因子CF的大小控制着该阶段优化的随机性与多样性，理论上种群q越大，该阶段的随机性越好，但同样会进入盲目搜索，通过多次实验，选择出种群大小取20最合适。同样，如图 4所示，CF自适应因子的参数s在一定取值下与搜索效率成反比，参数s取值越大，对于小规模实例的搜索效率更高，但随机性与多样性也会大打折扣，陷入局部最优。参数s也不能过大，在大规模实例下会造成CF一直处于低概率，该阶段难以收敛。经多次实验，参数s取0.25，迭代次数取350次效率最好。这也是图 9中约束密度p < 0.11时，低维度下算法很难求解和时间效率低的原因。不过在大规模难解实例中，可以有效大幅降低违反约束，如图 9中可以看出，在p>0.11后，DMPA-SMCH算法在大规模难解实例上的时间效率反超RSA算法与GSA算法。第三阶段为MCH算法局部搜索。事实上，在小规模易解区域实例，第三阶段完全可以搜索得到可行解，理论上该阶段搜索时间越长，算法精度会越高。但当在大规模接近可满足相变点时，违反约束数越往后越复杂难解，算法类似穷举搜索，效率低下。故在第三阶段后再加上SMCH算法，在MCH算法中融入退火策略，当温度T变大时，以较大概率选择随机扰动，当温度T变小时进行针对性扰动。目的是以一定的概率选择赋值新解，避免搜索的过度随机性，在退火后期能进行针对性搜索，提高整体算法的搜索效率。经实验验证，当迭代次数达到400次时，MCH的优化效率最好。然而这并不意味着在DMPA算法的第三阶段的MCH搜索是多余的，相反MCH搜索同等重要，它为最后的SMCH搜索提供一个双层保障，保证了整个算法求解的稳定性。

5 结语

对于RB模型实例的求解探索提出了DMPA-SMCH算法，在整体规模接近可满足相变区域提升了元启发式的求解效率。但根据实验数据发现，在难解复杂实例集上，该离散化方法还有一定的优化空间，包括对CF自适应因子的定义公式，也并不是最理想的定义。在未来的工作中，可以进一步优化离散化方法，探索启发式优化算法对RB模型实例的求解研究。