本文参考Google OR-Tools官网文档介绍OR-Tools的使用方法。

1 约束满足问题

1.1 CSP定义

在前一篇文章中我提到Google OR-Tools中解决整数规划问题有MP Solver和CP Solver两种工具，但是只介绍了MP Solver，而这篇文章则会介绍CP Solver。之所以会有这两种工具，是因为虽然都是解决整数规划，各自面对的问题类型还是有区别的，CP Solver一般用于约束满足问题(Constraint Satisfaction Problem)，简称CSP。

一个约束满足问题包括有限变量集、每个变量的有限论域和有限约束集，每条约束限制了变量集赋值的组合，约束满足问题的解是为每个变量赋一个对应论域上的值，使之满足所有的约束。一般情况，解约束满足问题的目标是只要找到一个解就行，不过也可能要找出所有可行解中最优的那个。

举个简单的例子，比如一个人需要穿鞋、衬衫和裤子，其中每样东西都有一定的选择：鞋有运动鞋、皮鞋两双；衬衫有红色和白色两件；裤子有蓝色、灰色和牛仔裤三条。但是必须符合如下的规则：

运动鞋配牛仔裤；
皮鞋只能与灰裤子和白衬衫搭配；
白衬衫可以和牛仔裤或者蓝裤子搭配；
红衬衫只能与灰裤子搭配；

这个穿衣规则就是问题的约束，问题变量集合为{鞋，裤子，衬衫}，而每种物品可选择的样式就构成了变量的论域。

下面是CSP问题的数学定义：
用一个三元组 $(V,D,C)$ 来表示CSP，其中

$V$ 是由有限变量组成的集合 ${X_1,X_2,...,X_n}$ ；
$D$ 是一个函数，将每个变量映射到一个有限论域上，用 $D_{Xi}$ 表示变量 $X_i$ 在函数 $D$ 上的映射，即变量 $X_i$ 的论域；
$C$ 是一个约束的有限集合，而且在约束中出现的变量只能是 $V$ 中的变量。如果 $c\in C$ ，则约束 $c$ 表示的是其中出现的所有变量对应论域笛卡尔积的一个子集。

对于一个 $CSP(V,D,C)$ ：

如果一个 $CSP(V,D,C)$ 是 $k$ -可满足的，当且仅当 $V$ 中所有由 $k$ 个变量组成的子集都存在一组标识使得这组标识满足所有和它们相关的约束
一个有 $n$ 个变量的CSP是可满足的，如果它是 $n$ -可满足的

所以根据定义，如果一个CSP是可满足的，那么它至少存在一个解；如果一个CSP是不可满足的，则问题的解不存在。

乍看这个CSP的定义，会感觉很难和非线性规划的一般方程形式相关联，这个其实就是CSP问题和通常的整数规划问题的区别，还是那上面那个穿衣服的案例来分析，这个例子完全是用自然语言表示的各种逻辑关系构成的，没有人数、长度这类很显然的变量存在，所以这种问题很难直接用一般的非线性规划方程表达出来。当然只要变量定义得合适还是可以写成一般的非线性规划形式的，例如我们如果给每个衣物关联一个布尔变量，表示是否穿这个衣物，那么各种约束就变成了这些布尔变量的逻辑关系了。以我个人的理解看，CSP和一般的整数规划相比，主要是以下这些特点：

问题由很多带逻辑关系的约束构成
变量通常会是布尔变量
很难用线性表达式写成一般的方程形式
很多时候约束带有If-Than的条件信息

在现实中，很多问题(通常也会称为组合优化问题)都属于CSP，例如生产调度、物流规划、路由选择、资源分配等等，解空间规模一般比较大，而且需要符合很多现实规则约束，因此其求解难度一般都很大。

1.2 CSP相关算法

针对CSP的算法主要分成了两类，一类属于非完备性算法，另一类属于完备性算法。非完备性算法包括各种启发式或元启发式算法，例如模拟退火、遗传算法，它们的目的不在于找到精确的最优解，只要在合理的时间找到一个相对最优解，这些算法不会遍历整个解空间，因此如果原问题有解，那么通常非完备性算法可以快速找到一个合理的解，但是如果原问题无解，非完备性算法无法证明无解，换句话说即使跑了很长时间算法得到的解不合理，我们也不能认定该问题无解；而完备性算法则以各种搜索和推理技术为主，称其为完备是因为这些算法实际上是搜索了整个解空间的，因此如果这些算法得出的结论是无解，那么该问题就是无解的，同时如果问题有解，得到的最终解也肯定是最优的。从我个人的使用经验来猜测，CP Solver内部应该是采用的完备性算法，使用CP Solver求解CSP，如果约束有冲突，问题实际是无解的，CP Solver会第一时间判断出来。关于完备性算法，我参阅了一点资料，以下是我个人的一些理解。

目前通常的完备性CSP算法是将以回溯为主的搜索技术与以约束传播为主的推理技术相结合，通过搜索技术确保解空间被完整遍历，同时结合推理技术动态地省去冗余的子空间，从而加快求解速度。基础的回溯就类似于深度优先遍历，按照一定的次序依次实例化变量，如果进行到某一个变量时发现该变量的论域内没有任何值可以满足约束，则直接回退到上一个变量，取下一个值进行实例化。例如 $X$ 和 $Y$ 是实例化次序相邻的两个变量，它们的论域都是 $\{1,2,3\}$ ，并且需要满足约束 $X+Y=1$ ，假设 $X$ 实例化为了1，然后轮到 $Y$ 实例化，结果发现 $Y$ 取任何值都无法满足约束，则停止 $X$ 实例化为1的分支，回退到 $X$ ，开始 $X$ 实例化为2的分支。

CSP_1.png

回溯可以实现部分的剪枝效果，但是存在大量冗余的实例化操作，比如这个例子里，显然问题是无解的，把 $X$ 的论域都实例化一遍是没有意义的。针对这个问题，需要引入约束传播算法，就是在求解过程中使用约束信息删除变量中一些不参与解的值从而减小搜索空间。约束传播算法以相容性技术为主，而弧相容则是相容性技术最基础的概念。

对于一个二元CSP的约束图中的某一个弧 $(X_i, X_j)$ ，称它是弧相容的(Arc Consistency, AC)，当且仅当对于 $X_i$ 论域中能满足 $X_i$ 上一元约束的每一个值 $a$ ，都在 $X_j$ 的论域中存在一个值 $b$ ，使得 $b$ 满足 $X_j$ 上一元约束，并且 $(a,b)$ 满足 $X_i$ 和 $X_j$ 上的二元约束 $C(X_i,X_j)$
一个CSP是弧相容的，当且仅当它的约束图中每一条弧都是弧相容的
对于问题P在进行弧相容检查时，如果得到某一变量的论域为空，则此问题无解

我们还是用上面那个例子，不过把约束改成 $X+Y=3$ ，然后对 $X$ 和 $Y$ 进行相容性检查，变量 $X$ 中值 3 在变量 $Y$ 中没有值，能够满足约束 $X+Y=3$ ，因此被移走，同理变量 $Y$ 中值 3 也被移走，最后论域得到了缩减:

CSP_2.png

再把约束改回进行相容性检查，则的论域将为空，则该问题无解，可以直接退出计算了。所以通过约束传播，一方面可以迅速判断是否有解，另一方面可以大大缩小每一次搜索的解空间，从而达到加快计算的效果。

从上面对搜索和推理算法的简单介绍可以看到完备性算法的优点在于可以第一时间检查问题是否有解，并且一定能给出最优解，但缺点也很明显，虽然采用了各种减小搜索空间的技术，当问题的规模增大时，搜索空间仍然以指数级增长，这些技术对计算效率的提升也越来越小。

1.3 约束编程

按照定义，约束编程(Constraint Programming, CP)是围绕关系约束这一数学概念建立起来的方法论，是研究基于约束的组合优化问题的计算系统。实际生活中有很多复杂的组合优化问题，我们当然可以对不同的问题以合适的算法分别设计和开发计算系统，但是这样不具有通用性。约束编程或者说约束程序设计就是以高效的约束求解技术为核心，结合强说明性的问题描述方式建立的组合优化计算系统，对于使用者来说，可以用约束编程语言表达问题的任意约束，并且考虑新类型的约束时(特别是在面向对象框架中)，只需要对约束系统进行扩展。像Google OR-Tools的CP Solver就是一个CP系统，不仅内部提供了高效的求解CSP问题的算法，外部也提供了灵活的约束语言接口，用于构建适合自身项目的优化问题模型。

2 OR-Tools Demo

2.1 N皇后问题

N皇后问题可以描述为，是否可以在一个N X N的棋盘上放置N个皇后，使得它们之间互相不会攻击，所谓攻击指的是棋盘上同一行、同一列和同一对角线上不会有两个皇后。例如下面这张图就是4皇后的一个可行方案

sol_4x4_b.png

N皇后问题是一个典型的CSP问题，我们并不需要找到一个最优解，而是要找到满足约束的解。
我们来分析下如何对这个问题进行建模。首先需要确定决策变量是什么，对于CSP问题，确定决策变量是相当重要也是最具难度的步骤之一，好的决策变量可以方便约束的构建，反之不合适的决策变量会加大建立约束的难度，甚至影响计算的效率。对于N皇后问题，我们的一种方案是每个点位分配一个布尔变量，如果这个变量为True，表示第i行第j列放置了一个皇后。但这种方式并不是太好，对于N X N的棋盘，将会有N X N 个变量，无疑会加大搜索空间。另一种更好的方案是定义一个一维数组，表示第j列上的皇后在第i行，这样首先变量的数量只有N个，然后在变量阶段就已经隐含了一个约束：必须要有N个皇后被放置。
然后我们利用决策变量来定义约束：

同一列上不会有两个以上的皇后。这个约束已经隐含在变量定义里了，因为我们已经人为地用索引区分了
同一行上不会有两个以上的皇后。只要保证 $Q[j]$ 里没有两个以上的相同的值，或者说 $Q[j]$ 的值必须都不相同
同一对角线上不会有两个以上的皇后。这个稍微复杂点，假设两个皇后的坐标分别为 $(X_1,Y_1)$ 和 $(X_2,Y_2)$ ，则如果它们在对角线上，要么是 $\frac{Y_1-Y_2}{X_1-X_2}=1$ ，即 $X1-Y1=X2-Y2$ ，要么是 $\frac{Y_1-Y_2}{X_1-X_2}=-1$ ，即 $X1+Y1=X2+Y2$

CSP_3.png

因此我们只需要保证所有 $Q[j]$ 对应的两个数组 $Q[j]+j$ 和 $Q[j]-j$ 的值都不相同

2.2 代码

新建一个.Net Core应用，下载OR-Tools。在程序开头引用Google.OrTools.Sat库

using Google.OrTools.Sat;

首先我们定义一个CpModel对象。不同于之前线性规划，CP系统将模型model和求解器Solver在接口层分开了，这也符合CP系统的方法论。

            //Create CP Model
            var model = new CpModel();

然后定义决策变量数组

            //Create Variables
            List<IntVar> queens = new List<IntVar>();
            for(int j=0;j<N;j++)
            {
                queens.Add(model.NewIntVar(0, N - 1, $"queen_{j}"));
            }

接着定义约束一，不允许两个及以上的皇后出现在同一行。OR-Tools直接提供了AddAllDifferent接口实现一个数组元素互不相同的约束

            //All queens are in different rows
            model.AddAllDifferent(queens);

约束二，对角线上不允许两个及以上的皇后

            //No two queens can be on the same diagonal.
            List<IntVar> diag1 = new List<IntVar>();
            List<IntVar> diag2 = new List<IntVar>();
            for (int j = 0; j < N; j++)
            {
                var q1 = model.NewIntVar(0, 2 * N, $"diag1");
                diag1.Add(q1);
                model.Add(q1 == queens[j] + j);
                var q2 = model.NewIntVar(-N, N, $"diag2");
                diag2.Add(q2);
                model.Add(q2 == queens[j] - j);
            }
            model.AddAllDifferent(diag1);
            model.AddAllDifferent(diag2);

然后就可以定义Solver了

            //Create Solver
            var solver = new CpSolver();

但是在调用计算之前，我们可以先定义一个回调接口供Solver使用，方便获取所有的解

        public class QueenSolutionCallback:CpSolverSolutionCallback
        {
            private int solutionCount=0;
            private List<IntVar> queens;
            public QueenSolutionCallback(List<IntVar> queens)
            {
                this.queens = queens;
            }

            public override void OnSolutionCallback()
            {
                Console.WriteLine($"#Solution {solutionCount}:");
                for(int row=0;row<queens.Count;row++)
                {
                    string s = "";
                    for(int column=0; column < queens.Count;column++)
                    {
                        if(Value(queens[column])==row)
                        {
                            s += "Q ";
                        }
                        else
                        {
                            s += "_ ";
                        }
                    }
                    Console.WriteLine(s);
                }
                Console.WriteLine();
                solutionCount++;
            }
        }

最后计算结果

            //Get the result
            QueenSolutionCallback cb = new QueenSolutionCallback(queens);
            var status = solver.SearchAllSolutions(model, cb);

我们看一下4皇后问题的结果：

1.PNG

4个皇后时有2个解。再看下8皇后的情况：

2.PNG

这时将有92个解。
完整的程序

using System;
using System.Collections.Generic;
using Google.OrTools.Sat;

namespace Demo4
{
    class Program
    {
        static void Main(string[] args)
        {
            Nqueen(8);
        }

        static void Nqueen(int N)
        {
            //Create CP Model
            var model = new CpModel();

            //Create Variables
            List<IntVar> queens = new List<IntVar>();
            for(int j=0;j<N;j++)
            {
                queens.Add(model.NewIntVar(0, N - 1, $"queen_{j}"));
            }

            //All queens are in different rows
            model.AddAllDifferent(queens);

            //No two queens can be on the same diagonal.
            List<IntVar> diag1 = new List<IntVar>();
            List<IntVar> diag2 = new List<IntVar>();
            for (int j = 0; j < N; j++)
            {
                var q1 = model.NewIntVar(0, 2 * N, $"diag1");
                diag1.Add(q1);
                model.Add(q1 == queens[j] + j);
                var q2 = model.NewIntVar(-N, N, $"diag2");
                diag2.Add(q2);
                model.Add(q2 == queens[j] - j);
            }
            model.AddAllDifferent(diag1);
            model.AddAllDifferent(diag2);

            //Create Solver
            var solver = new CpSolver();


            //Get the result
            QueenSolutionCallback cb = new QueenSolutionCallback(queens);
            var status = solver.SearchAllSolutions(model, cb);

            Console.WriteLine(status);
        }

        public class QueenSolutionCallback:CpSolverSolutionCallback
        {
            private int solutionCount=0;
            private List<IntVar> queens;
            public QueenSolutionCallback(List<IntVar> queens)
            {
                this.queens = queens;
            }

            public override void OnSolutionCallback()
            {
                Console.WriteLine($"#Solution {solutionCount}:");
                for(int row=0;row<queens.Count;row++)
                {
                    string s = "";
                    for(int column=0; column < queens.Count;column++)
                    {
                        if(Value(queens[column])==row)
                        {
                            s += "Q ";
                        }
                        else
                        {
                            s += "_ ";
                        }
                    }
                    Console.WriteLine(s);
                }
                Console.WriteLine();
                solutionCount++;
            }
        }
    }
}

Google OR-Tools(四) 约束编程 Constraint Programming