永磁同步电机模型的随机动力学分析

Stochastic Dynamics Analysis of Permanent Magnet Synchronous Motor Model

摘要: 首先基于同步旋转坐标的永磁同步电机模型引入随机噪声项，建立永磁同步电机系统的非线性微分方程，运用中心流形定理对随机系统进行降维处理，再运用随机平均法得到Ito ̂微分方程；其次根据Lyapunov指数法讨论了系统的局部稳定性、应用奇异边界理论讨论了全局稳定性并得到稳定性条件；根据拟不可积Hamilton系统随机平均法分析了系统的D-分岔和P-分岔；最后选取合适的分岔参数进行数值模拟，并验证了分岔位置。

Abstract: Firstly, the nonlinear differential equation of the PMSM systems was established by introducing random terms into the PMSM model based on synchronous rotating coordinates. The stochastic center manifold is used to reduce the dimension of the stochastic systems, and then the Ito differ-ential equations are obtained by the stochastic average method. Then, by using Lyapunov index and singular boundary theory, the local and global stability of the systems are discussed, and the conditions of stability are obtained. The D-bifurcation and P-bifurcation behavior of the systems is analyzed based on the stochastic average method for the quasi-integrable Hamiltonian systems. Finally, the appropriate bifurcation parameters are selected for numerical simulation, and the bi-furcation positions are verified.

1. 引言

永磁同步电机是一种新型同步电机，因其高精度、高功率密度、高转矩惯性比、小体积等优越性能广泛应用于航空航天、船舶电力、电动汽车等领域，但许多学者在研究永磁同步电机时会忽略或近似波动因素(如忽略温度和阻尼影响、如空间磁场近似视为正弦分布等)，讨论确定系统下的动力学行为，文献 [1] 在理想化条件下研究了车用永磁同步电机模型的抗扰和控制，其分析过程存在诸多限制，得不到更接近实际的结果；文献 [2] 在讨论永磁同步电机模型时指出了如磁链耦合等影响，并根据影响因素对系统进行了优化，但其优化是基于结构上进行的，并未对影响参数改进处理；文献 [3] 在确定的永磁同步电机模型下，引入了噪声干扰，但其动力学行为分析仍是基于规律性的正弦扰动。我们已知永磁同步电机工作时影响因素，电阻会受温度影响；摩擦系数会受阻尼影响；电感会受外部磁场影响。为了更好的提高永磁同步电机的鲁棒性需考虑上述影响参数，我们通过引入随机项对系统进行随机激励耗散的建立，便可获得一种更接近实际情况的描述思路。

目前，针对永磁同步电机模型，应用Hamilton理论进行分析的研究非常稀少，而随机动力系统的研究已有大量的理论方法，文献 [4] [5] [6] [7] 介绍了随机激励耗散Hamilton的知识，文献 [8] 介绍了调速器系统的随机分岔行为，这样便有机械系统引入随机噪声的参考可寻。本文在文献 [2] 的模型上引入了随机项，建立了随机参数激励下的永磁同步电机模型，首先运用中心流形定理对系统进行降维化简，再运用随机平均法得到Itô微分方程；接着讨论了系统稳定性并得出稳定性条件；最后进行分岔分析且选取合适的参数进行数值模拟。

2. 随机参数激励下的永磁同步电机模型建立

在交直轴电感近似相等的情形下，永磁同步电机基于同步旋转坐标下的数学模型 [2] 表示如下：

$\{\begin{array}{l}\frac{\text{d}w}{\text{d}t}=\frac{3P\phi }{2J}{i}_q-\frac{B}{J}w-\frac{T_L}{J}\\ \frac{\text{d}{i}_d}{\text{d}t}=-\frac{R}{L}{i}_d+Pw{i}_q+\frac{1}{L}{u}_d\\ \frac{\text{d}{i}_q}{\text{d}t}=-\frac{R}{L}{i}_q-Pw{i}_d-\frac{P\phi }{L}w+\frac{1}{L}{u}_q\end{array}$ (1)

系统变量及参数的物理意义： $i_d,i_q$ 为d-q轴定子电流， $u_d,u_q$ 为d-q轴定子电压，w为转子机械角速度，R为定子电阻，L为定子电感， $T_L$ 为负载转矩，J为转动惯量，B为粘滞摩擦系数，P为极对数， $\phi$ 为转子磁链，所涉及参数均为正数，求得非负平衡点为： $E_1\left(w_1^{\text{*}},i_{d_1}^{\text{*}},i_{q_1}^{\text{*}}\right)$。

$w_1^{*}=X_1=M-\frac{N}{M}-\frac{a}{3},i_{d_1}^{*}=\frac{2BL{X}_1^2+2T_{L}L{X}_1+3\phi u_q}{3\phi R},i_{q_1}^{*}=\frac{2B{X}_1+2T_L}{3P\phi }$

$M={\left(-8a^3+36ab-108c+12{\left(12a^{3}c-3a^2{b}^2-54abc+12b^3+81c^2\right)}^{\frac{1}{2}}\right)}^{\frac{1}{3}}/6$

$N=6\left(b/3-a^2/9\right)$

$a=\frac{T_L}{B},b=\frac{2B{R}^2+3\phi L{P}^2+3{\phi }^2{P}^{2}R}{2B{L}^2{P}^2},c=\frac{2T_L{R}^2-3\phi PR{U}_q}{2B{L}^2{P}^2}$

考虑系统(1)中参数电阻R，电感L，摩擦系数B，转子磁链 $\phi$ 会受内部温度外部磁场、阻尼、耦合等因素影响变得不稳定，把影响因素视为白噪声，对其引入随机项：

$R\to R-\sigma \xi \left(t\right),L\to L-\sigma \xi \left(t\right),B\to B-\sigma \xi \left(t\right)$，得到一个新的随机非线性微分方程：

$\{\begin{array}{l}\frac{\text{d}w}{\text{d}t}=\frac{3P\phi }{2J}{i}_q-\frac{B}{J}w-\frac{T_L}{J}+\sigma \left(w\left(t\right)-w_2^{*}\right)\xi \left(t\right)\\ \frac{\text{d}{i}_d}{\text{d}t}=-\frac{R}{L}{i}_d+Pw{i}_q+\frac{1}{L}{u}_d+\sigma \left(i_d\left(t\right)-i_{d_2}^{*}\right)\xi \left(t\right)\\ \frac{\text{d}{i}_q}{\text{d}t}=-\frac{R}{L}{i}_q-Pw{i}_d-\frac{P\phi }{L}w+\frac{u_q}{L}+\sigma \left(\left(i_q\left(t\right)-i_{q_2}^{*}\right)-\left(i_d\left(t\right)-i_{d_2}^{*}\right)\right)\xi \left(t\right)\end{array}$ (2)

其中 $\sigma$ 为白噪声强度， $\xi \left(t\right)$ 为零均值白噪声，其中 $\langle \xi \left(t\right)\rangle =0,\langle \xi \left(t\right)\xi \left(t+\tau \right)\rangle =\delta \left(\tau \right)$ 接着令 $x_1=w\left(t\right)-w_1^{\text{*}},x_2=i_d\left(t\right)-i_{d_1}^{\text{*}},x_3=i_q\left(t\right)-i_{q_1}^{\text{*}}$，系统(2)可化为：

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_1={\alpha }_1{x}_1+{\gamma }_1{x}_3+{\eta }_1+\sigma x_1\xi \left(t\right)\\ {\stackrel{\dot{}}{x}}_2={\alpha }_2{x}_1+{\beta }_2{x}_2+{\gamma }_2{x}_3+P{x}_1{x}_3+{\eta }_2+\sigma x_2\xi \left(t\right)\\ {\stackrel{\dot{}}{x}}_3={\alpha }_3{x}_1+{\beta }_3{x}_2+{\gamma }_3{x}_3+P{x}_1{x}_2+{\eta }_3+\sigma \left(x_3-x_2\right)\xi \left(t\right)\end{array}$ (3)

其中： ${\alpha }_1=-\frac{B}{J},{\beta }_1=0,{\gamma }_1=\frac{3P\phi }{2J},{\eta }_1=\frac{3P\phi }{2J}{i}_{q_2}^{*}-\frac{B}{L}{w}_2^{*}-\frac{T_L}{J}$

${\alpha }_2=P{i}_{q_2}^{*},{\beta }_2=-\frac{R}{L},{\gamma }_2=P{w}_2^{*},{\eta }_2=P{w}_2^{*}{i}_{q_2}^{*}-\frac{R}{L}{i}_{d_2}^{*}+\frac{u_d}{L}$

${\alpha }_3=-P{i}_{q_2}^{*}-\frac{P\phi }{L},{\beta }_3=-P{w}_2^{*},{\gamma }_3=-\frac{R}{L},{\eta }_3=-\frac{R}{L}{i}_{q_2}^{*}-P{w}_2^{*}{i}_{q_2}^{*}-\frac{P\phi }{L}{w}_2^{*}+\frac{u_d}{L}$

3. 永磁同步电机模型的处理

据上述推导关系知，系统(3)在 $E_0\left(0,0,0\right)$ 处的性质等价系统(2)在平衡点 $E_1\left(w_1^{\text{*}},i_{d_1}^{\text{*}},i_{q_1}^{\text{*}}\right)$ 的性质，我们可把系统(2)在 $E_1$ 的稳定性转化为系统(3)在 $E_0$ 处的来讨论。模型(3)的特征方程为：

${\lambda }^3+p{\lambda }^2+q\lambda +l=0$ (4)

其中： $p=-{\alpha }_1-{\beta }_2-{\gamma }_3,q={\alpha }_1{\beta }_2+{\alpha }_1{\gamma }_3-{\alpha }_3{\gamma }_1+{\beta }_2{\gamma }_3-{\beta }_3{\gamma }_2$

$l={\alpha }_1{\beta }_3{\gamma }_2-{\alpha }_2{\beta }_3{\gamma }_1-{\alpha }_1{\beta }_2{\gamma }_3+{\alpha }_3{\beta }_2{\gamma }_1$

求得方程(4)的特征值为： ${\lambda }_{1,2}=\frac{\epsilon }{2}\pm \sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}$, ${\lambda }_3=-p-\epsilon$

令 $q_0=-\frac{l}{p}$， $q=q_0+\epsilon \frac{{\left(p+\epsilon \right)}^2-\frac{l}{p}}{p+\epsilon }$，当参数 $\epsilon$ 趋近于零时，， ${\lim }_{\epsilon \to 0}{\lambda }_3=-p-\epsilon =-p$ 此时点 $E_0\left(0,0,0\right)$ 将退化为系统(3)的临界焦点，具有局部不变流行的性质。

令 $x_1,x_2,x_3$ 与 $\mu ,\eta ,\zeta$ 的关系式如下：

$x_1=\mu +\zeta$

$x_2=\frac{\frac{\epsilon }{2}{\beta }_2-\frac{l}{p+\epsilon }}{{\beta }_2^2+\epsilon {\beta }_2-\frac{l}{p+\epsilon }}\mu +\frac{{\beta }_2\sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}}{{\beta }_2^2+\epsilon {\beta }_2-\frac{l}{p+\epsilon }}\eta +\frac{p+\epsilon -{\alpha }_1}{-p-\epsilon +{\beta }_2}\zeta$

$\begin{array}{c}{x}_3=\frac{\frac{\epsilon }{2}{\beta }_2{\gamma }_3-\frac{l}{p+\epsilon }\left(\frac{\epsilon }{2}+{\beta }_2+{\gamma }_3\right)}{\left[{\beta }_2^2+\epsilon {\beta }_2-\frac{l}{p+\epsilon }\right]\left[{\gamma }_2^2+\epsilon {\gamma }_2-\frac{l}{p+\epsilon }\right]}\mu +\frac{\left({\beta }_2{\gamma }_3+\frac{l}{p+\epsilon }\right)\sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}}{\left[{\beta }_2^2+\epsilon {\beta }_2-\frac{l}{p+\epsilon }\right]\left[{\gamma }_2^2+\epsilon {\gamma }_2-\frac{l}{p+\epsilon }\right]}\eta \\ +\frac{{\beta }_3\left(p+\epsilon -{\alpha }_1\right)}{\left(-p-\epsilon +{\beta }_2\right)\left(-p-\epsilon +{\gamma }_3\right)}\zeta \end{array}$

则系统(3)可变为：

$\{\begin{array}{l}\frac{\text{d}\mu }{\text{d}t}=\frac{\epsilon }{2}\mu -\sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}\eta +f_1\left(\mu ,\eta ,\zeta \right)+\sigma \left[m_7\mu +m_8\eta +m_9\zeta \right]\xi \left(t\right)\\ \frac{\text{d}\eta }{\text{d}t}=\sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}\mu +\frac{\epsilon }{2}\eta +f_2\left(\mu ,\eta ,\zeta \right)+\sigma \left[n_7\mu +n_8\eta +n_9\zeta \right]\xi \left(t\right)\\ \frac{\text{d}\zeta }{\text{d}t}=-\left(p+\epsilon \right)\zeta +f_3\left(\mu ,\eta ,\zeta \right)+\sigma \left[r_7\mu +r_8\eta +r_9\zeta \right]\xi \left(t\right)\end{array}$ (5)

其中： $f_1\left(\mu ,\eta ,\zeta \right)=m_1{\mu }^2+m_2\mu \eta +m_3\mu \zeta +m_4{\eta }^2+m_5\eta \zeta +m_6{\zeta }^2$

$f_2\left(\mu ,\eta ,\zeta \right)=n_1{\mu }^2+n_2\mu \eta +n_3\mu \zeta +n_4{\eta }^2+n_5\eta \zeta +n_6{\zeta }^2$

$f_3\left(\mu ,\eta ,\zeta \right)=r_1{\mu }^2+r_2\mu \eta +r_3\mu \zeta +r_4{\eta }^2+r_5\eta \zeta +r_6{\zeta }^2$

此时系统(3)在原点附近具有局部不变流行：

$w_{loc}^c\left(o\right)=\left\{\left(\mu ,\eta ,\zeta \right)|\zeta =h\left(\mu ,\eta \right),\left|\mu \right|+\left|\eta \right|\le 1,h\left(0,0\right)=h_{\mu }\left(0,0\right)=h_{\eta }\left(0,0\right)=0\right\}$

对系统(5)降维可近似表示为下式：

$\{\begin{array}{l}\frac{\text{d}\mu }{\text{d}t}=\frac{\epsilon }{2}\mu -\sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}\eta +\psi \left(\mu ,\eta \right)+\sigma \left[m_7\mu +m_8\eta +m_9\zeta \right]\xi \left(t\right)\\ \frac{\text{d}\eta }{\text{d}t}=\sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}\mu +\frac{\epsilon }{2}\eta +\Phi \left(\mu ,\eta \right)+\sigma \left[n_7\mu +n_8\eta +n_9\zeta \right]\xi \left(t\right)\end{array}$ (6)

其中： $\psi \left(\mu ,\eta \right)=m_1{\mu }^2+m_2\mu \eta +m_3\mu h+m_4{\eta }^2+m_5\eta h+m_6{h}^2$

$\Phi \left(\mu ,\eta \right)=n_1{\mu }^2+n_2\mu \eta +n_3\mu h+n_4{\eta }^2+n_5\eta h+n_6{h}^2$

且h满足 $\zeta$ 的全微分 $\stackrel{\dot{}}{\zeta }=h_{\mu }\stackrel{\dot{}}{\mu }+h_{\eta }\stackrel{\dot{}}{\eta }$，所以将 $\zeta$ 的全微分代入(6)式得，

$\begin{array}{c}{h}_{\mu }\stackrel{\dot{}}{\mu }+h_{\eta }\stackrel{\dot{}}{\eta }=h_{\mu }\left(\frac{\epsilon }{2}\mu -\sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}\eta +\psi \left(\mu ,\eta \right)+\sigma \left[m_7\mu +m_8\eta +m_9\zeta \right]\xi \left(t\right)\right)\\ +h_{\eta }\left(\sqrt{-\frac{l}{p+\epsilon }-\frac{{\epsilon }^2}{4}}\mu +\frac{\epsilon }{2}\eta +\Phi \left(\mu ,\eta \right)+\sigma \left[n_7\mu +n_8\eta +n_9\zeta \right]\xi \left(t\right)\right)\\ =-\left(p+\epsilon \right)h+r_1{\mu }^2+r_2\mu \eta +r_3\mu \zeta +r_4{\eta }^2+r_5\eta \zeta +r_6{\zeta }^2+\sigma \left[r_7\mu +r_8\eta +r_9\zeta \right]\xi (t)\end{array}$

我们要求解的中心流形可用下式来逼近：

$h\left(\mu ,\eta \right)=h_1{\mu }^2+h_2\mu \eta +h_3\mu \sigma \xi \left(t\right)+h_4{\eta }^2+h_5\eta \sigma \xi \left(t\right)+h_6{\sigma }^2\xi {\left(t\right)}^2+\cdots$

代入上述方程两边相同幂系数合并得：

$h_1=\frac{p^2{r}_1+2q\left(r_1+r_4\right)-p\sqrt{q}{r}_2}{p\left(4q+p^2\right)},h_2=\frac{2\sqrt{q}\left(r_1-r_4\right)+p{r}_2}{4q+p^2},h_3=\frac{p{r}_7-\sqrt{q}{r}_8}{p^2+q}$

$h_4=\frac{p^2{r}_4+2q\left(r_1+r_4\right)+p\sqrt{q}{r}_2}{p\left(4q+p^2\right)},h_5=\frac{p{r}_8-\sqrt{q}{r}_7}{p^2+q},h_6=0$

忽略高阶项，可进一步得到系统(7)

$\{\begin{array}{l}\frac{\text{d}\mu }{\text{d}t}=-\sqrt{p}\eta +m_1{\mu }^2+m_2\mu \eta +m_3{h}_1{\mu }^3+\left(m_3{h}_2+m_5{h}_1\right){\mu }^2\eta \\ +\left(m_3{h}_4+m_5{h}_2\right)\mu {\eta }^2+m_5{h}_4{\eta }^3+\sigma \left(m_7\mu +m_8\eta \right)\xi \left(t\right)\\ \frac{\text{d}\eta }{\text{d}t}=\sqrt{p}\mu +n_1{\mu }^2+n_2\mu \eta +n_3{h}_1{\mu }^3+\left(n_3{h}_2+n_5{h}_1\right){\mu }^2\eta \\ +\left(n_3{h}_4+n_5{h}_2\right)\mu {\eta }^2+n_5{h}_4{\eta }^3+\sigma \left(n_7\mu +n_8\eta \right)\xi \left(t\right)\end{array}$ (7)

令 $\mu =r\cos \theta ,\eta =r\sin \theta$，进行极坐标变换系统(7)可转化为：

$\{\begin{array}{l}{g}_1\left(r,\theta ,{\xi }_t\right)=\stackrel{\dot{}}{r}\left(t\right)=\frac{\text{d}\mu }{\text{d}t}\cos \theta +\frac{\text{d}\eta }{\text{d}t}\sin \theta \\ g_2\left(r,\theta ,{\xi }_t\right)=\stackrel{\dot{}}{\theta }\left(t\right)=-\frac{\text{d}\mu }{\text{d}t}\frac{\sin \theta }{r}+\frac{\text{d}\eta }{\text{d}t}\frac{\cos \theta }{r}\end{array}$ (8)

其中：

$\begin{array}{c}{g}_1\left(r,\theta ,{\xi }_t\right)=m_1{r}^2{\cos }^3\theta +m_3{h}_1{r}^3{\cos }^4\theta +\left(m_2+n_1\right)r^2{\cos }^2\theta \sin \theta +m_4{n}_5{r}^3\sin \theta \\ +m_1{n}_3{r}^3{\cos }^3\theta +\left(m_3{h}_2+m_5{h}_1\right)r^3{\cos }^3\theta \sin \theta +m_5{h}_4{r}^3\cos \theta \sin \theta \\ +n_3{h}_4{r}^3\cos \theta \sin \theta +\left(m_3{h}_4+m_5{h}_2+n_3{h}_2+n_5{h}_1\right)r^3{\cos }^2\theta \sin \theta \\ +\sigma \left(m_7{\cos }^2\theta +\left(m_7+n_8\right)\cos \theta \sin \theta +n_8{\sin }^2\theta \right)\xi (t)\end{array}$

$\begin{array}{c}{g}_2\left(r,\theta ,{\xi }_t\right)=\sqrt{q}+n_1{r}^2{\cos }^3\theta +n_3{h}_1{r}^3{\cos }^4\theta -\left(n_1-m_2\right)r{\cos }^2\theta \sin \theta \\ +n_3{h}_2{r}^3{\cos }^3\theta \sin \theta -n_5{h}_4{r}^2{\sin }^4\theta +\left(n_5{h}_1-m_3{h}_1\right)r^2{\cos }^2\theta \sin \theta \\ +\left(n_5{h}_4-m_3{h}_4-m_5{h}_2\right)r^2\cos \theta {\sin }^3\theta -m_{2}r\cos \theta {\sin }^2\theta \\ +\left(n_3{h}_4+n_5{h}_2-m_3{h}_2+m_5{h}_1\right)r^2{\cos }^2\theta {\sin }^2\theta \\ +\sigma \left(n_7{\cos }^2\theta -\left(m_7-n_8\right)\cos \theta \sin \theta -m_8{\sin }^2\theta \right)\xi (t)\end{array}$

在解决系统(8)时，精确解不易求出，但可基于Khasminskii极限定理 [4] [9] [10] [11]，当白噪声强度充分小时，其反应 $\left\{r\left(t\right),\theta \left(t\right)\right\}$ 弱收敛于一个二维的马尔可夫扩散过程，运用随机平均法，系统(8)可得到Itô随机微分方程：

$\{\begin{array}{l}\text{d}r=m_r\text{d}t+{\sigma }_{11}\text{d}w\\ \text{d}\theta =m_{\theta }\text{d}t+{\sigma }_{21}\text{d}w\end{array}$ (9)

对于上述二维马尔可夫扩散过程，计算它的转换概率密度，方程(9)表达为：

$\{\begin{array}{l}\text{d}r=\left(\frac{v_1}{8}r+\frac{v_2}{8}{r}^3\right)\text{d}t+\sqrt{\frac{v_3}{8}{r}^2}\text{d}w\\ \text{d}\theta =\left(\sqrt{q}+\frac{v_4}{8}{r}^2\right)\text{d}t+\sqrt{v_{5}r}\text{d}w\end{array}$ (10)

其中：

$v_1=5{\sigma }^2{m}_7^2+5{\sigma }^2{n}_8^2+3{\sigma }^2{m}_8^2+3{\sigma }^2{n}_7^2+6{\sigma }^2{m}_8^2{n}_7^2-2{\sigma }^2{m}_7{n}_8$

$v_2=3m_3{h}_1+m_3{h}_4+m_5{h}_5+n_3{h}_2+n_5{h}_1+3n_3{h}_4$

$v_3=3{\sigma }^2{m}_7^2+3{\sigma }^2{n}_8^2+{\sigma }^2{n}_7^2+{\sigma }^2{m}_8^2+2{\sigma }^2{m}_8^2{n}_7^2+2{\sigma }^2{m}_7^2{n}_8^2$

$v_4=6n_3{h}_1+2n_5{h}_2-2m_3{h}_2+2n_3{h}_4-2m_5{h}_1-6h_4{h}_2$

$v_5=\frac{1}{4}{\sigma }^2\left(m_7+n_8\right)\left(n_7-m_8\right)$

根据扩散矩阵，可以得到平均振幅是一维的马尔可夫过程，如下式(11)：

$\text{d}r=\left(\frac{v_1}{8}r+\frac{v_2}{8}{r}^3\right)\text{d}t+{\left(\frac{v_3}{8}{r}^2\right)}^{\frac{1}{2}}\text{d}w$ (11)

这样，在讨论模型的稳定性和随机分岔时获得一个思路，即根据参数的随机概率密度分布，分析在概率密度意义下平均振幅的稳定性和临界分岔。

4. 稳定性分析

4.1. 局部稳定性

讨论局部稳定时，基于Oseledec乘性遍历理论 [4] 和最大Lyapunov指数法，随机方程平凡解在概率意义下稳定的充要条件为：最大Lyapunov指数 $\lambda <0$。

我们讨论线性Itô随机微分方程的稳定性，首先令 $v_2=0$，即讨论 $v_2=0$ 时方程(11)的稳定性，此时方程变为：

$\text{d}r=\frac{v_1}{8}r\text{d}t+{\left(\frac{v_3}{8}{r}^2\right)}^{\frac{1}{2}}\text{d}w$ (12)

对微分方程(12)求解得，

$r\left(t\right)=r\left(0\right)\exp ({\displaystyle {\int }_0^t\left(m^{\prime }\left(0\right)-{\left({\sigma }^{\prime }\left(0\right)\right)}^2/2\right)\text{d}s}+{\displaystyle {\int }_0^t{\sigma }^{\prime }\left(0\right)\text{d}w(s))}$

其中 $m\left(0\right)=\frac{v_1}{8}$, $\sigma \left(0\right)={\left(\frac{v_3}{8}\right)}^{\frac{1}{2}}$

根据拟不可积Hamilton系统理论 [5]，线性Itô随机微分方程的Lyapunov指数近似表示为：

$\lambda ={\lim }_{t\to +\infty }\frac{1}{t}\ln \left(\Vert x\left(t,x_0\right)\Vert \right)$。

这样可得方程(12)的Lyapunov指数：

$\lambda =\underset{t\to +\infty }{\lim }\frac{1}{t}\ln {\left(r\left(t\right)\right)}^{\frac{1}{2}}=\frac{m^{\prime }\left(0\right)-{\left({\sigma }^{\prime }\left(0\right)\right)}^2/2}{2}=\frac{1}{2}\left(\frac{v_1}{8}-\frac{v_3}{16}\right)$

① 当 $\frac{v_1}{8}-\frac{v_3}{16}<0$，即 $\lambda <0$ 时，线性Itô随机微分方程(11)是局部稳定的，即原系统(2)在平衡点处是稳定的；

② 当 $\frac{v_1}{8}-\frac{v_3}{16}>0$，即 $\lambda >0$ 时，线性Itô随机微分方程(11)是局部不稳定的，即原系统(2)在平衡点处是不稳定的。

4.2. 全局稳定性

在分析局部稳定性时，应用Oseledec乘性遍历理论和最大Lyapunov指数判定方法，但该理论方法只适用于做局部稳定性判断，接下来分析全局稳定性，我们可用奇异边界理论 [4] [5] 对其全局稳定性做分析判断。

首先介绍奇异边界理论，我们已知r，此外还有扩散指数 ${\alpha }_r$ ；漂移指数 ${\beta }_r$ ；特征指数 $c_r$。

I) 如果有 $r=0$，即 $\sigma \left(r\right)=0$，系统为第一类奇异边界，第一类边界下：

当 ${\alpha }_r={\beta }_r+1,{\beta }_r=1$ 时， $c_r>1$ 是排斥自然边界， $c_r<1$ 是吸引自然边界， $c_r=1$ 是严格自然边界；

II) 如果有 $r=+\infty$，即 $m\left(r\right)=+\infty$，系统为第二类奇异边界，第二类边界下：

当 ${\beta }_r={\alpha }_r-1$， ${\beta }_r\le 1$ 时， $c_r>-{\beta }_r$ 是排斥自然边界， $c_r<-1$ 是吸引自然边界， $c_r\le -{\beta }_r$ 且 $c_r\ge -1$ 是严格自然边界；

当 ${\beta }_r>{\alpha }_r-1$ 时， ${\beta }_r>1$ 是进入边界。

上述为奇异边界理论介绍，接下来我们分两种情形进行分析，即 $v_2=0,v_2\ne 0$ 情形。

1) 当 $v_2=0$ 时，系统如上述(12)， $r=0$ 时系统(12)为第一类奇异边界，经过计算可以得， ${\alpha }_r=2$， ${\beta }_r=1$，

$c_r=\underset{r\to 0^{+}}{\lim }\frac{2m\left(r\right)\cdot r^{{\alpha }_r-{\beta }_r}}{{\left(\sigma \left(r\right)\right)}^2}=\underset{r\to 0^{+}}{\lim }\frac{2\frac{v_1}{8}{r}^2}{\frac{v_3}{8}{r}^2}=\frac{2v_1}{v_3}$，则有： $\frac{v_1}{v_3}>\frac{1}{2}$ 是排斥自然边界， $\frac{v_1}{v_3}<\frac{1}{2}$ 是吸引自然边界， $\frac{v_1}{v_3}=\frac{1}{2}$

是严格自然边界；

$r=+\infty$ 时系统(12)为第二类奇异边界，经过计算可以得， ${\alpha }_r=2$， ${\beta }_r=1$，

$c_r=-\underset{r\to +\infty }{\lim }\frac{2m\left(r\right)\cdot r^{{\alpha }_r-{\beta }_r}}{{\left(\sigma \left(r\right)\right)}^2}=-\underset{r\to +\infty }{\lim }\frac{2\frac{v_1}{8}{r}^2}{\frac{v_3}{8}{r}^2}=-\frac{2v_1}{v_3}$ 则有， $\frac{v_1}{v_3}<\frac{1}{2}$ 是排斥自然边界， $\frac{v_1}{v_3}>\frac{1}{2}$ 是吸引自然边界， $\frac{v_1}{v_3}=\frac{1}{2}$ 是严格自然边界。

综上可得出结论：

① 当 $\frac{v_1}{v_3}<\frac{1}{2}$， $r=0$ 为吸引自然边界， $r=+\infty$ 是排斥自然边界，则Itô微分方程(11)有概率意义下稳定的平凡解，说明系统(2)在白噪声激励下平衡点处是稳定的；

② 当 $\frac{v_1}{v_3}>\frac{1}{2}$， $r=0$ 为排斥自然边界， $r=+\infty$ 是吸引自然边界，则Itô微分方程(11)有概率意义下不稳定的平凡解，说明系统(2)在白噪声激励下平衡点处是不稳定的；

③ 当 $\frac{v_1}{v_3}=\frac{1}{2}$，和 $r=+\infty$ 为严格自然边界。

2) 当 $v_2\ne 0$ 时，系统为方程(11)， $r=0$ 系统(11)是第一类奇异边界，经过计算可得到， ${\alpha }_r=2$， ${\beta }_r=1$，

$c_r=\underset{r\to 0^{+}}{\lim }\frac{2m\left(r\right)r^{{\alpha }_r-{\beta }_r}}{{\left(\sigma \left(r\right)\right)}^2}=\underset{r\to 0^{+}}{\lim }\frac{2\left(\frac{v_1}{8}r+\frac{v_3}{8}{r}^3\right)r}{\frac{v_3}{8}{r}^2}=\frac{2v_1}{v_3}$，则有： $\frac{v_1}{v_3}>\frac{1}{2}$ 是排斥自然边界， $\frac{v_1}{v_3}<\frac{1}{2}$ 是吸引自然边界， $\frac{v_1}{v_3}=\frac{1}{2}$ 是严格自然边界； $r=+\infty$ 系统(11)为第二类奇异边界，经计算 ${\alpha }_r=2$， ${\beta }_r=3$，则 $r=+\infty$ 是二类奇异边界中的进入边界。

综上可以得出结论：

① 当 $\frac{v_1}{v_3}<\frac{1}{2}$， $r=0$ 是吸引自然边界，且 $r=+\infty$ 是进入边界，则Itô微分方程(11)有概率意义下稳定的平凡解，说明系统(2)在白噪声激励下平衡点处是稳定的；

② 当 $\frac{v_1}{v_3}>\frac{1}{2}$， $r=0$ 时为排斥自然边界，而 $r=+\infty$ 是进入自然边界，则Itô微分方程(11)有概率意义下不稳定的平凡解，说明系统(2)在白噪声激励下平衡点处是不稳定的；

③ 当 $\frac{v_1}{v_3}=\frac{1}{2}$， $r=0$ 时为严格自然边界。

5. 随机分岔分析

随机分岔理论 [4] [5] 是通过研究参数的变化，对随机动态系统进行性态分析，我们将应用拟不可积Hamilton系统随机平均法来解决随机分岔动力学行为。

5.1. 随机D-分岔分析

当Itô随机微分方程(11)中 $v_2=0$ 时，方程即为：

$\text{d}r=\frac{v_1}{8}r\text{d}t+{\left(\frac{v_3}{8}{r}^2\right)}^{\frac{1}{2}}\text{d}w$ (13)

由于 $v_3=0$ 方程变为不发生分岔的确定系统，我们讨论 $v_3\ne 0$ 时的情形，令 $m\left(0\right)=\left(\frac{v_1}{8}-\frac{v_3}{16}\right)r$， $\sigma \left(0\right)={\left(\frac{v_3}{8}\right)}^{\frac{1}{2}}r$，方程(13)生成的连续动态系统为：

$\mu \left(t\right)=x+{\displaystyle {\int }_0^{t}m\left(\mu \left(s\right)x\right)\text{d}s}+{\displaystyle {\int }_0^t\sigma \left(\mu \left(s\right)x\right)\circ \text{d}w}$ (14)

式(14)为方程(13)以x为初值的唯一强解。当 $m\left(0\right)=0,\sigma \left(0\right)=0$ 时，0是 $\mu$ 的固定点，对此固定点， $\circ \text{d}w$ 是Stratonovich的随机意义下的激参，设 $m\left(r\right)$ 有界，则对 $r\ne 0$ 满足椭圆性条件的 $\sigma \left(0\right)\ne 0$，保证最多有一个平稳概率，则对应的FPK方程为：

$\frac{\partial p}{\partial t}=-\frac{\partial }{\partial r}\left\{\left(\frac{v_1}{8}r\right)p\right\}-\frac{{\partial }^2}{\partial r^2}\left\{\left(\frac{v_3}{8}{r}^2\right)p\right\}$ (15)

其中，令 $\frac{\partial p}{\partial t}=0$ 得到与方程(15)相应的平稳概率密度：

$p\left(r\right)=c\left|{\sigma }^{-1}\left(r\right)\right|\exp \left({\displaystyle {\int }_0^r\frac{2m\left(v\right)}{{\sigma }^2\left(v\right)}\text{d}v}\right)$ (16)

这样，则该动态系统可能有不动点、非平凡平稳两种平稳状态，且二者不变测度 ${\delta }_0$ 和 $\epsilon$ 的密度分别为 $\delta \left(x\right)$ 和(16)式，此时还需确定两不变测度的Lyapunov指数，解线性Itô随机微分方程(13)为：

$r\left(t\right)=r\left(0\right)\exp \left[{\displaystyle {\int }_0^t\left[\stackrel{\dot{}}{m}\left(r\right)+\frac{\sigma \left(r\right)\ddot{\sigma }\left(r\right)}{2}\right]\text{d}s}+{\displaystyle {\int }_0^t\stackrel{\dot{}}{\sigma }\left(r\right)\text{d}w}\right]$ (17)

已知系统测度的Lyapunov指数：

${\lambda }_{\mu }={\lim }_{t\to +\infty }\frac{1}{t}\ln \Vert r\left(t\right)\Vert$ (18)

根据公式 $\sigma \left(0\right)=0,\stackrel{\dot{}}{\sigma }\left(0\right)=0$ 得，不动点状态的Lyapunov指数：

$\begin{array}{c}{\lambda }_{\mu }\left({\delta }_0\right)=\underset{t\to +\infty }{\lim }\frac{1}{t}\ln \Vert r\left(0\right)\Vert +\stackrel{\dot{}}{m}\left(0\right){\displaystyle {\int }_0^t\text{d}s}+\stackrel{\dot{}}{\sigma }\left(0\right){\displaystyle {\int }_0^t\text{d}w\left(s\right)}\\ =\stackrel{\dot{}}{m}\left(0\right)+\stackrel{\dot{}}{\sigma }\left(0\right)\underset{t\to +\infty }{\lim }\frac{w_r\left(t\right)}{t}=\stackrel{\dot{}}{m}\left(0\right)=\frac{v_1}{8}-\frac{v_3}{16}\end{array}$ (19)

设 $\stackrel{\dot{}}{\sigma }$ 有界， $\stackrel{\dot{}}{m}\left(r\right)+\sigma \left(r\right)\ddot{\sigma }\left(r\right)$ 可积，将式(17)代入方程(18)，非平凡平稳状态的Lyapunov指数为：

$\begin{array}{c}{\lambda }_{\mu }\left(\epsilon \right)=\underset{t\to +\infty }{\lim }\frac{1}{t}\left[{\displaystyle {\int }_0^t\left[\stackrel{\dot{}}{m}\left(r\right)+\sigma \left(r\right)\ddot{\sigma }\left(r\right)\right]\text{d}s}\right]={\displaystyle \int \left[\stackrel{\dot{}}{m}\left(r\right)+\frac{\sigma \left(r\right)\ddot{\sigma }\left(r\right)}{2}\right]p\left(r\right)\text{d}r}\\ =-2{\displaystyle \int {\left(\frac{m\left(r\right)}{\sigma \left(r\right)}\right)}^{2}p\left(r\right)\text{d}r}\stackrel{\dot{}}{m}\left(0\right)=-2v_3^{\frac{3}{2}}\left(\frac{v_1}{8}-\frac{v_3}{16}\right)\exp \left[\left(\frac{16}{v_3}\right)\left(\frac{v_1}{8}-\frac{v_3}{16}\right)\right]\end{array}$ (20)

令 $\kappa =\frac{v_1}{8}-\frac{v_3}{16}$，则根据二者不变测度的Lyapunov指数可判定，当 $\kappa \le 0$ 时前者不变测度稳定，当 $\kappa >0$

时后者不变测度稳定，即系统(2)的一个D-分岔在 $\kappa ={\kappa }_D=0$ 处发生。

我们进一步处理系统(16)可以得到：

$p_{st}\left(r\right)=c{r}^{\frac{2\left(v_1-v_3\right)}{v_3}}=o\left(r^{\chi }\right)r\to 0$ (21)

其中：c为归一化常数， $\chi =\frac{2\left(v_1-v_3\right)}{v_3}$，当 $\chi <-1$，即 $\frac{v_1}{v_3}<\frac{1}{2}$ 时， $p_{st}\left(r\right)$ 是一个 $\delta$ 函数，当 $-1<\chi <0$，即 $\frac{v_1}{v_3}>\frac{1}{2}$ 时， $r=0$ 是 $p_{st}\left(r\right)$ 在该区间的最大值点，所以随机系统(2)发生D-分岔的临界条件为， $\chi =-1$，即 $\frac{v_1}{v_3}=\frac{1}{2}$。

当Itô随机微分方程(11)中 $v_2\ne 0$ 时，令 $\Phi ={\left(-\frac{v_3}{8}r\right)}^{\frac{1}{2}}$, $v_2<0$，则方程(11)变为：

$\text{d}\Phi =\left(\kappa \cdot \Phi -{\Phi }^3\right)\text{d}t+{\left(\frac{v_3}{8}\right)}^{\frac{1}{2}}\Phi \circ \text{d}w$ (22)

当 $v_3=0$ 时，系统为确定分岔-叉形分岔，当 $v_3\ne 0$ 时，讨论 $\kappa <0$ 时，有唯一密度为 $\delta \left(x\right)$ 不变测度 ${\delta }_0$ ； $\kappa >0$ 时，动态系统的平衡状态相应的不变测度分别为 ${\delta }_{\kappa }$ 和非平凡平稳状态测度，密度为：

$P_{\kappa }^{+}\left(\Phi \right)=\{\begin{array}{l}{G}_{\kappa }{\Phi }^{\frac{2\kappa }{v_3/8}-1}\exp \left(-\frac{8{\Phi }^2}{v_3}\right),\Phi >0\\ 0,\Phi \le 0\end{array}$ (23)

$F_{\kappa }^{-}\left(\Phi \right)=P_{\kappa }^{+}\left(-\Phi \right)$ (24)

其中 $G_{\kappa }^{-1}=\Gamma \left(\frac{8\kappa }{v_3}\right)\cdot {\left(\frac{v_3}{8}\right)}^{-\frac{8\kappa }{v_3}}$，则根据式(19)得 ${\delta }_{\kappa }$ 的Lyapunov指数为：

${\lambda }_{\mu }\left({\delta }_{\kappa }\right)=\kappa =\frac{v_1}{8}-\frac{v_3}{16}$ (25)

根据方程(20)，(23)，(25)可得两个非平凡平稳状态测度 ${\epsilon }_{\chi }^{\pm }$ 的Lyapunov指数为：

${\lambda }_{\mu }\left({\epsilon }_{\chi }^{\pm }\right)=-2\kappa =-\frac{v_1}{4}+\frac{v_3}{8}$ (26)

所以 $\kappa ={\kappa }_D=0$ 时，随机系统(2)的发生D-分岔。

5.2. 随机P-分岔分析

令 $\Phi ={\left(-\frac{v_3}{8}r\right)}^{\frac{1}{2}}$, $v_2<0$，则方程(11)变为：

$\text{d}\Phi =\left(\left(\frac{v_1}{8}-\frac{v_3}{16}\right)\Phi -{\Phi }^3\right)\text{d}t+{\left(\frac{v_3}{8}\right)}^{\frac{1}{2}}\Phi \circ \text{d}w$ (27)

可得到方程(27)的FKP方程为：

$\frac{\partial p}{\partial t}=-\frac{\partial }{\partial r}\left[\left(\frac{v_1}{8}r-\frac{v_3}{2v_2}r-r^3\right)p\left(r\right)\right]-\frac{v_3}{2v_2}\frac{{\partial }^2}{\partial r^2}\left(r^{2}p\left(r\right)\right)$ (28)

其初值条件为： $v_2=0$, $p\left(r,t|r_0,t_0\right)\to \delta \left(r-r_0\right)$, $t\to t_0$, $p\left(r,t|r_0,t_0\right)$ 为扩散过程的转移概率密度， $r\left(t\right)$ 的不变测度是平稳概率密度，若令 $\frac{\partial p}{\partial t}=0$，则退化系统的解为：

$p\left(r\right)=\frac{\exp \left(\frac{r^2{v}_2}{v_3}\right)r^{1-\frac{v_1{v}_2}{4v_3}}}{\Gamma \left[-\frac{v_1{v}_2}{8v_3}\right]{\left(-\frac{v_2}{v_3}\right)}^{\frac{v_1{v}_2}{8v_3}}}$

对系统进行P-分岔分析时，我们的思路是由线性Itô随机微分方程的概率密度 $p\left(r\right)$ 计算不变测度极值，进而分析分岔。根据Namachivay’s理论 [4] [10] [12]，当噪声 $\sigma$ 趋于0时，的极值趋表现稳态特征，当系统响应 $r\left(t\right)$ 是遍历的，由乘积遍历性定理得 $p\left(r\right)$ 是样本轨线在 的邻域内访问时间的度量。

若 $p\left(r\right)$ 在 $r_0$ 处为极大值，轨线在 $r_0$ 的邻域内停留时间较长，说明 $r_0$ 在概率意义下是稳定点；若 $p\left(r\right)$ 在 $r_0$ 处为极小值，轨线在 $r_0$ 的邻域内停留时间较短，说明 $r_0$ 在概率意义下是不稳定点。

6. 数值模拟

综上可知原系统的响应过程为 $\left(\mu \left(t\right),\eta \left(t\right)\right)$，且 $r\left(t\right)=\sqrt{{\mu }^2\left(t\right)+{\eta }^2\left(t\right)}$，所以我们主要分析 $\left(\mu \left(t\right),\eta \left(t\right)\right)$ 的联合概率密度 $p\left(\mu ,\eta \right)$ 即可，在实际意义下，赋值参数 $m_i,n_i\left(i=1,2,3,\cdots ,8\right)$，使得 $\frac{v_2}{v_3}=-4$，在此情况

下，给出平稳概率密度函数和联合概率密度函数模拟图像。

若设 $\rho =-1+v_1$，图1中左图表示，取不同的负值 $\rho$，平稳概率密度函数是单调递减函数，且随 $\rho$ 的增大图像减势变弱，右图表示 $\rho =-0.02$ 时联合概率密度图像呈现单峰尖状；图2中左图表示在 $\rho =0$ 附近，随 $\rho$ 值增大平稳概率密度函数从递减变为有极值状态，右图表示 $\rho =0$ 时联合概率密度图像呈现单峰状，说明此时并未发生分岔；图3中左图表示，取不同的正值 $\rho$，平稳概率密度函数是极值函数，且随 $\rho$ 的增大平稳概率密度函数极值峰谷趋于平缓，右图表示 $\rho >0$ 时联合概率密度函数呈现火山口形状，则系统在临界点 $\rho =0$ 后发生了分岔。

Figure 1. The left figure shows stationary probability density function image of $\rho =-0.22,-0.12,-0.02$ ; The right figure shows joint probability density function image of $\rho =-0.02$

图1. 左图表示 $\rho$ 分取−0.22，−0.12，−0.02时平稳概率密度函数图像；右图表示 $\rho =-0.02$ 联合概率密度函数图像

Figure 2. The left figure shows stationary probability density function image of $\rho =-0.02,0,0.02$ ; The right figure shows joint probability density function image of $\rho =0$

图2. 左图表示 $\rho$ 分取−0.02，0，0.02时平稳概率密度函数图像；右图表示 $\rho =0.32$ 联合概率密度函数图像

Figure 3. The left figure shows stationary probability density function image of $\rho =0.12,0.22,0.32$ ; The right figure shows joint probability density function image of $\rho =0.32$

图3. 左图表示 $\rho$ 分取0.12，0.22，0.32时平稳概率密度函数图像；右图表示 $\rho =0.32$ 联合概率密度函数图像

7. 结论

本文讨论了在随机参数激励下永磁同步电机的稳定性和分岔，把随机参数视为白噪声，运用随机激耗散的Hamilton理论，对系统的动力学行为进行分析，结合实际意义选取分岔参数 $\rho$，根据数值模拟结果得出当 $\rho >0$ 时系统发生分岔行为，说明此时永磁同步电机系统有可能会变得不稳定。