其实CPG网络原理非常简单，我们先介绍一下线性变换中的旋转变换：

旋转变换

在x，y平面的逆时针旋转矩阵如下，角度为${\theta }_0$的向量经过该旋转矩阵的变换则变为${\theta }_0+\theta$的角度
$$X({\theta }_0+\theta )=\left[\begin{array}{cc}cos(\theta )& -sin(\theta )\\ sin(\theta )& cos(\theta )\end{array}\right]\ast X({\theta }_0)$$
如果要求变化量，主要做差就找到了：
$$dX=X({\theta }_0+\theta )-X({\theta }_0)$$
令旋转矩阵为A，转化一下就是：

$$dX=A\ast X({\theta }_0)-I\ast X({\theta }_0)$$
得到X变换为dX的矩阵B为：
$$B=A-I=\left[\begin{array}{cc}cos(\theta )-1& -sin(\theta )\\ sin(\theta )& cos(\theta )-1\end{array}\right]$$
当然，实际上由于每次旋转角度$\theta$都是非常小的，所以$cos(\theta )$非常接近于1

旋转向量的收敛

如果只是上面那样旋转，那么效果就是一个半径$r=\sqrt{x^2+y^2}$的圆，x，y分别是X向量的坐标，由于只是旋转没有长度变化所以园是一个标准圆。现在我想让长度收敛到$L$怎么办，很显然我们会用上矩阵特征值：
$$R=\left[\begin{array}{cc}\lambda 1& 0\\ 0& \lambda 2\end{array}\right]$$
只要让特征值随着半径变化即可，这里需要用到负反馈知识：

如果向量长度r > 预期收敛长度L 则 特征值$\lambda <1$
如果向量长度r < 预期收敛长度L 则 特征值$\lambda >1$

所以我们可以设置$\lambda =k\ast (L-r)+1$，只不过得要求这个k非常小，必须让$\lambda$能保持1附近，如在 [0.5,1.5]区间

这里的也就是拉伸的比例系数，k越大，收敛的越快，当然太大也不行，会发生震荡，想必这很容易理解。
另外为什么$\lambda$不能小于0我就不说了，线性代数应该都熟悉。

接下来我们只要让R*A即可：
$$X2=\left[\begin{array}{cc}\lambda & 0\\ 0& \lambda \end{array}\right]\ast \left[\begin{array}{cc}cos(\theta )& -sin(\theta )\\ sin(\theta )& cos(\theta )\end{array}\right]\ast X$$
我们将这个矩阵记为C，则有：
$$C=\left[\begin{array}{cc}\lambda cos(\theta )& -\lambda sin(\theta )\\ \lambda sin(\theta )& \lambda cos(\theta )\end{array}\right]$$
得到X变换为dX的矩阵D为
$$D=C-I=\left[\begin{array}{cc}\lambda cos(\theta )-1& -\lambda sin(\theta )\\ \lambda sin(\theta )& \lambda cos(\theta )-1\end{array}\right]$$
我之前令$\lambda =k\ast (L-r)+1$，带入得：

$\lambda cos(\theta )-1=k\ast (L-r)cos(\theta )+[cos(\theta )-1]$

假设单位旋转角$\theta$非常小趋近于0，根据$cos(\theta )$的导数特性，$[cos(\theta )-1]$就基本为0，这一项就可以忽略

$li{m}_{x\to 0}co{s}^{\prime }(x)=0$

所以在旋转角非常小的情况下有：

$\lambda cos(\theta )-1=k\ast (L-r)cos(\theta )$

我推导的变换矩阵就变成了
$$D=\left[\begin{array}{cc}k\ast (L-r)cos(\theta )& -\lambda sin(\theta )\\ \lambda sin(\theta )& k\ast (L-r)cos(\theta )\end{array}\right]$$

HOPF振荡器

先列一下公式：
$\frac{dx}{dt}=\alpha (\mu -r^2)x-\omega y$

$\frac{dy}{dt}=\omega x+\alpha (\mu -r^2)y$

我们抽取一下形成矩阵：
$$dX=\left[\begin{array}{cc}\alpha (\mu -r^2)& -\omega \\ \omega & \alpha (\mu -r^2)\end{array}\right]\ast X$$
有了我之前的模型，这个模型理解起来也非常容易，对照一下有：

$k\ast (L-r)cos(\theta )=\alpha (\mu -r^2)$
$\lambda sin(\theta )=\omega$

很明显两者实际上异曲同工，只不过我用的是L1范数，而他用的是L2范数。为了方便计算，不用开根号，我也用L2范数好了：

$\lambda =k\times (\mu -r^2)+1$

$k\times cos(\theta )\times (\mu -r^2)=\alpha (\mu -r^2)$
$\lambda sin(\theta )=\omega$

我之前说过，k是特征值的拉伸系数，决定了长度L变化的快慢，实际上也就是收敛到圆环速度。我们再看HOPF的公式中的系数说明：

幅值：$A=\sqrt{\mu }$
​周期：$T=\frac{2\pi }{\omega }$
​收敛到圆环的速度：$\alpha 值$

第一个很好理解，第三个也很好理解，和我的模型都一样，但第二个还有点不一样，我们再来取极限分析：

$li{m}_{\theta \to 0}sin(\theta )=\theta$
这个高数的等价无穷小，为什么我就不过多解释了，因为sin函数的导数在0附近是等于1的。

由于假设k取值较小，则$\lambda$在1附近，我们按1算，所以我们近似得到：

$\lambda sin(\theta )=\theta$

单位旋转角为$\theta$，一周是2π，所以周期T 自然是$\frac{2\pi }{\theta }$

所以说在各种系数非常小时，我的模型就蜕变为HOPF振荡器了。

注意：使用我的模型请注明出处，商业用途先征得我的同意。

相变

由于CPG通常用于机器人上，假设机器人的脚有一半时间是用来支撑地面的，另一半时间用来摆动，就是从地上提起来，换到另一个地方。这两段时间我们称为摆动相和支撑相。

其实这个在我的模型里非常好解决，就是角速度$\theta$在不同区间，发生改变，这个改变可以参考绝对角度，比如角在[-90,90]区间，$\theta$变大一些，在[90，270]区间，$\theta$变小一些。

这个在CPG模型里，就不是特别好解决，引入了一个$\beta$项和a参数，支撑相和摆动相的频率分别为${\omega }_{st}$和${\omega }_{sw}$：

实际上也就是角速度是一个关于y的函数，我们来看一下系数$c=\frac{1}{e^{ay}+1}$

当y>0时，$e^{ay}变得非常大$，系数c非常小，区间为[ 0 , 0.5]
当y<0时，$e^{ay}变得非常小$，系数c非常大，区间为[ 0.5 ,1]

另一个系数$c1=\frac{1}{e^{-ay}+1}$，其实和c一样，相当于$c1=\frac{1}{e^{a(-y)}+1}$。这两个系数的作用，就是人y在变化的时候，使得角速度$\omega$趋向于${\omega }_{st}$或者${\omega }_{sw}$，原来式子改写为

$\omega =c1\ast {\omega }_{st}+c\ast {\omega }_{sw}$

当y>0时，系数c<0.5，c1>0.5，$\omega$趋向于${\omega }_{st}$
当y<0时，系数c>0.5，c1<0.5，$\omega$趋向于${\omega }_{sw}$
当y=0时，系数c=0.5，c1=0.5，$\omega$为两者取平均

分析完系数的作用，我们来分析${\omega }_{st}=\frac{1-\beta }{\beta }{\omega }_{sw}$，我们移项一下得到：

$\beta \ast {\omega }_{st}=(1-\beta )\ast {\omega }_{sw}$

我们假设有$1={\omega }_{st}+{\omega }_{sw}$，实际上我们换一个方式来理解：

${\omega }_{st}=1-{\omega }_{sw}$

也就是说$\beta$和${\omega }_{st}$是对偶形式，对原式变换有

$(1-{\omega }_{sw})\ast \beta =(1-\beta )\ast {\omega }_{sw}$

两边各占半个周期，总周期为两边之和

调节$\beta$越大，相当于调节支撑相的${\omega }_{st}$越大(即 y<0时变化速度越快)，如图所示：

由于向量是逆时针旋转，所以黄色部分是y>0时的状态，绿色部分是y<0时的状态。

黄色部分角速度趋向于${\omega }_{sw}$（摆动相），绿色部分趋向于${\omega }_{st}$(支撑相)

耦合

其中${\theta }_j=2\pi ({\phi }_i-{\phi }_j)$

这个是相位差，比如${\theta }_j=45度$，则说明向量Xi比Xj多转了45度，也就说Xj逆时针再旋转45度，将会和Xi的角度重合

假设我有两个旋转向量X1,X2，相位差 45度，那么上式就变为：

$\frac{dx1}{dt}=\alpha (\mu -r^2)x1-\omega y1+[cos(45)x2-sin(45)y2]$

$\frac{dy1}{dt}=\omega x1+\alpha (\mu -r^2)y1+[sin(45)x2+cos(45)y2]$

我们把尾部的耦合项分离一下有：
$d{X}_1=D{X}_1+M{X}_2$

我们来看一下后面的耦合项的含义：
$$\delta X=M{X}_2=\left[\begin{array}{cc}cos(45)& -sin(45)\\ sin(45)& cos(45)\end{array}\right]\ast X_2$$
我们便抽取出了逆时针旋转45度的旋转矩阵，我们问为什么能这样做呢？

我们之前提到过，X2旋转45度，将会和X1的角度重合

不过这种耦合方式会引起联合反馈震荡，很容易会出相位震荡问题，如果要我来设计耦合模型，起码会换成如下形式：

$d{X}_1=D(a_1\ast X_1+a_2\ast ||X_1||\ast \frac{M{X}_2}{||X_2||})$

其中a1取较大值，后面的a2,a3都取较小值，并且有a1+a2+a3… =1

注意：使用我的模型请注明出处，商业用途先征得我的同意。

我的模型大集会

X变换为dX的矩阵D为
$$d{X}_1=D{X}_1=\left[\begin{array}{cc}\lambda cos(d\theta )-1& -\lambda sin(d\theta )\\ \lambda sin(d\theta )& \lambda cos(d\theta )-1\end{array}\right]X_1$$

令$\lambda =k\ast (L-r)+1$，其中$r=||X_1||$

若带旋转耦合项，则如下：
$d{X}_1=D(a_1\ast X_1+a_2\ast ||X_1||\ast \frac{M{X}_2}{||X_2||})$

其中a1取较大值，后面的a2,a3都取较小值，并且有a1+a2+a3… =1

当然，我们要首先明确使用耦合的目的。

一个稳态系统必然有稳态解$[X(t_0)，X(t_1)......]$，这里的$X(t)$是状态空间表示法，$X(t)=[x_1，y_1，x_2，y_2,......]$，由当前时刻所有耦合的向量构成。

耦合的目的就是当系统受到干扰脱离稳态空间，会随时间自适应的回到稳态空间。

假设某个分量$x_1(t_0)$由于外界干扰变成了$x_1(t_1)$，其他项如$y_1(t_0)$也会自适应的变化到$y_1(t_1)$。

注意：使用我的模型请注明出处，商业用途先征得我的同意。