论文学习笔记（2）：Multi-Robot Coordination for Estimation and Coverage of Unknown Spatial Fields

论文学习笔记的第二篇，具体篇目“Multi-Robot Coordination for Estimation and Coverage of Unknown Spatial Fields”，作者为Alessia Benevento，Maria Santos等

一 涉及的算法介绍

1 centroidal Voronoi tessellation (CVT)

1.1 Voronoi图及其定义

Voronoi图由一组由连接两邻点线段的垂直平分线组成的连续多边形组成。在论文中用于计算机器人的分布情况。

一个Voronoi多边形（如下公式的Vi）内的任一点到构成该多边形的控制点的距离小于到其他多边形控制点的距离：
$${\mathcal{V}}_i\triangleq \left\{q\in {\mathbb{R}}^2:\forall j\ne i{\Vert q-x_i\Vert }_2⩽{\Vert q-x_j\Vert }_2\right\}$$

Vi可以看作是机器人i的区域，它覆盖了此区域。

1.2 CVT

每个Voronoi多边形区域Vi的质心被定义为
$$c_i(\varphi )\triangleq \frac{1}{{\int }_{{\mathcal{V}}_i}\varphi (q)dq}{\int }_{{\mathcal{V}}_i}q\varphi (q)dq$$

关于密度函数Φ的位置
{ c i ( ϕ ) } i = 1 n \left\{c_{i}(\phi)\right\}_{i=1}^{n} { ci​(ϕ)}i=1n​

被称为CVT。

1.3 CVT的应用

CVT是为了便于设计一个分布式控制律还有一个通信方案，使机器人能最小化所谓的位置成本
$$\ell (\varphi ,\mathbf{x})\triangleq \sum _{i=1}^n{\int }_{q\in {\mathcal{V}}_i}{\Vert q-x_i\Vert }_2^2\varphi (q)dq$$位置成本体现了机器人集体移动到事件重要区域相对位置的需求，定义为Φ(·)。当密度函数Φ(·)已知的时候，位置成本的局部最小值（不是唯一的）由Voronoi多边形单元的质心组成，这可以使用连续时间的 Lloyd’s algorithm来实现。
关于Lloyd’s algorithm可以参考：S. Lloyd, “Least squares quantization in pcm,” IEEE transactions on information theory, vol. 28, no. 2, pp. 129–137, 1982.

2 贝叶斯估计

贝叶斯估计是概率论中学过的内容，这里对它和极大似然估计（MLE），最大后验概率估计（MAP）做一个区分。

2.1 最大似然估计（MLE）

最大似然估计的思想是使得观测数据（样本）发生概率最大的参数就是最好的参数。
如对于正态分布求最大似然，似然函数为：
$$L\left(\mu ,{\sigma }^2\right)=\prod _{i=0}^n\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(x_i-\mu \right)}^2}{2{\sigma }^2}}$$对其取对数得
$$\ln L\left(\mu ,{\sigma }^2\right)=-\frac{n}{2}\ln (2\pi )-\frac{n}{2}\ln \left({\sigma }^2\right)-\frac{1}{2{\sigma }^2}\sum _{i=0}^n{\left(x_i-\mu \right)}^2$$分别对$\mu ,{\sigma }^2$取偏导，就可以解得
$$\{\begin{array}{l}\hat{\mu }=\frac{1}{n}{\sum }_{i=0}^n{x}_i=\bar{x}\\ \widehat{{\sigma }^2}=\frac{1}{n}{\sum }_{i=0}^n{\left(x_i-\bar{x}\right)}^2\end{array}$$所以最大似然估计的求解步骤：
1.确定似然函数
2.将似然函数转换为对数似然函数
3.求对数似然函数的最大值（求导，解似然方程）

2.2 最大后验概率估计(MAP)

MAP与MLE的区别：
最大似然估计认为使P(X | θ)最大的θ就是最好的，是因为它把θ看做一个常值，只是其值未知。而最大后验概率认为θ是一个随机变量，符合某种概率分布，即先验分布，即为P（θ），则使P（θ）P（X | θ）取最大时的对应的θ就是最好的。

由于P（X）是已知的（其实我们对X分布毫无兴趣），我们需要最大化的函数：
$$\frac{P(X\mid \theta )P(\theta )}{P(X)}=P(\theta \mid X)$$其中P(θ | X）即为后验分布。
最大后验概率估计的公式表示：
$$\underset{\theta }{\mathrm{argmax}}P(\theta \mid X)=\underset{\theta }{\mathrm{argmax}}\frac{P(X\mid \theta )P(\theta )}{P(X)}\propto \underset{\theta }{\mathrm{argmax}}P(X\mid \theta )P(\theta )$$举一个抛硬币的例子，连抛十次硬币，其中6次正面朝上，4次反面朝上。在例子中，由于均值为0.5的概率最大，故用均值0.5，方差0.1的高斯分布来描述θ的先验概率。计算出来的最佳θ值为0.52。如果用均值0.6去实验，最佳θ值为0.6附近。这就说明了合理的先验分布假设非常重要，一个错误的先验分布会得到完全错误的结果。

最大后验概率估计的求解步骤：
1.确定参数的先验分布以及似然函数
2.确定参数的后验分布函数
3.将后验分布函数转换为对数函数
4.求对数函数的最大值（求导，解方程）

2.3 贝叶斯估计

贝叶斯估计与MAP的区别：
贝叶斯估计是最大后验估计的进一步扩展，贝叶斯估计同样假定θ是一个随机变量，但贝叶斯估计并不是直接估计出θ的某个特定值，而是估计θ的分布，这是贝叶斯估计与最大后验概率估计不同的地方。

贝叶斯估计的公式在概率论中已经学过：
$$P(\theta \mid X)=\frac{P(X\mid \theta )P(\theta )}{{\int }_{\Theta }P(X\mid \theta )P(\theta )d\theta }$$

一般来说，我们无法直接对公式分母中的积分进行计算，所以确定一个合适的先验概率就非常重要。这里的先验概率应该是后验概率的共轭先验。
（具体到模型中：

• 二项分布参数的共轭先验是Beta分布
• 多项式分布参数的共轭先验是Dirichlet分布
• 指数分布参数的共轭先验是Gamma分布
• ⾼斯分布均值的共轭先验是另⼀个⾼斯分布
• 泊松分布的共轭先验是Gamma分布。）

注意，贝叶斯估计要解决的不是如何估计参数，而是用来估计新测量数据出现的概率，对于新出现的数据x~：
$$P(\tilde{x}\mid X)={\int }_{\Theta }P(\tilde{x}\mid \theta )P(\theta \mid X)d\theta ={\int }_{\Theta }P(\tilde{x}\mid \theta )\frac{P(X\mid \theta )P(\theta )}{P(X)}d\theta$$

贝叶斯估计的求解步骤：
1.确定参数的似然函数
2.确定参数的先验分布，应是后验分布的共轭先验
3.确定参数的后验分布函数
4.根据贝叶斯公式求解参数的后验分布

3 贝叶斯优化（Bayesian Optimization，BO）

贝叶斯优化的过程中会涉及到高斯过程（Gaussian Process，GP），推到了再介绍。

3.1 适用范围

贝叶斯优化适用两类情况：

• 需要优化的function计算起来非常费时费力
• 需要优化的function没有导数信息

3.2 基本步骤

比如我们要优化的function是
$$f:\mathcal{X}\to \mathbb{R}$$

为了方便可以假设这里的X为离散的。然后假设我们要解决的优化问题是
$$x^{\ast }=\arg \underset{x\in \mathcal{X}}{\max }f(x)$$BO需要不断进行迭代得到结果。在每一次迭代t=1,2,3,…,T中，我们选一个输入$x_t\in X$,然后观察对应的f的值$f\left(x_t\right)$。
多数情况下我们观察的结果都含有噪声，即
$$y_t=f\left(x_t\right)+ϵ,噪声\in \sim N\left(0,{\sigma }^2\right)$$我们把新观测到的$\left(x_t,y_t\right)$加入到观测数据中，然后进行下一次迭代。
由此可以发现，BO问题的核心在于如何选择每一次迭代里面要观测的$x_t$。这需要通过acquisition function 采集函数来决定，即
x t = arg ⁡ max ⁡ x ∈ X α t (