0 引　言

随着船舶永磁同步电机在船舶推进系统中的广泛应用，对于船舶永磁同步电机的控制性能有了更高的要求^[1]。然而在船舶永磁同步电机运行中，转动惯量会受到无法预知的因素影响，如负载扰动、参数变化等，这些外部扰动的变化会影响船舶推进系统的控制性能^[2]。

对于船舶推进系统，需要在线调整控制器参数，使系统具有良好的运行状态。转动惯量是速度环PI控制器的基本设计参数，所以根据运行工况实时调整转动惯量，能够有效抑制系统出现超调和振荡的情况^[3]。负载转矩的突变会导致转速瞬间的跌落或上升，因此需要对负载转矩观测，并作为电流补偿，使得系统具备抗扰动性能。转动惯量辨识和负载转矩观测是提高船舶推进系统性能的重点研究问题。

惯量辨识的方法通常分为在线和离线两大类。常见的离线惯量辨识方法为加减法^[4]，先是规划三角波速度指令，根据转速反馈和转矩反馈计算得到惯量值，这种方法得到的辨识值精度高，但是在工况复杂的情况下，系统的诸多参数发生变化，采用离线惯量辨识不能反映惯量的动态变化。因此，诸多学者结合系统辨识理论，提出在线惯量辨识方法，主要有最小二乘法、模型参考自适应法、扩展卡尔曼滤波算法、人工智能算法、梯度算法等。

文献[5]提出一种基于惯量辨识的自适应扰动观测器，解决了辨识中参数耦合引起的波动问题。文献[6]采用递推最小二乘法对多参数同时辨识，但辨识参数稳定性有待提高。文献[7]提出一种变周期最小二乘法对转动惯量在线辨识，解决了在低加速度下的辨识问题，提高算法的适用范围。文献[8]提出在线辨识转动惯量的扰动观测器，并把辨识的转动惯量和负载转矩作用到速度控制器中,但辨识精度有待提升。文献[9]用蚁群算法同时辨识转动惯量和负载转矩，但算法运行量大，收敛速度不快。

上述文献对惯量辨识算法都有所优化，仍存在船舶推进系统在复杂工况下辨识的惯量误差较大和无法收敛，其难点在于惯量辨识需要观测实时的负载，负载的观测又需要实际的惯量^[10-11]。为此本文提出梯度校正法与扩展卡尔曼滤波器（Extended Kalman -Filter, EKF）结合的在线惯量辨识算法，采用后向差分法对电机的运动方程离散化，将梯度校正法用于此模型^[12]。在惯量辨识的同时需要获取负载转矩，利用扩展卡尔曼滤波器把负载转矩当成状态变量实现参数估计。对存在负载扰动的工况，可以有效解决转动惯量和负载转矩的耦合问题，具有收敛速度快、精度高、抗扰动性好等特点，仿真和实验验证了算法的有效性^[13]。

1 PMSM的数学模型

在同步旋转坐标系$ d - q $中，其定子电压方程可以表示为：

$ \left\{ \begin{gathered} {u_d} = R{i_d} + {L_d}\frac{{\rm{d}}}{{{\rm{d}}t}}{i_d} - {\omega _e}{L_q}{i_q}，\\ {u_q} = R{i_q} + {L_q}\frac{{\rm{d}}}{{{\rm{d}}t}}{i_q} + {\omega _e}({L_d}{i_d} + {\psi _f}) 。\\ \end{gathered} \right. $

(1)

式中：$ {u_d} $和$ {u_q} $分别为定子电压的$ d - q $轴分量；$ {i_d} $和$ {i_q} $分别为定子电流的$ d - q $轴分量；$ R $为定子电阻；$ {\omega _e} $是电角速度；$ {L_d} $和$ {L_q} $分别为$ d - q $轴电感分量；$ {\psi _f} $代表永磁体磁链。

PMSM的电磁转矩方程为：

$ {T_e} = \frac{3}{2}{P_n}{i_q}[{i_d}({L_d} - {L_q}) + {\psi _f}]。$

(2)

式中：$ {P_n} $为电机极对数；$ {T_e} $为电磁转矩。

对于表贴式三相PMSM，定子电感满足$ {L_d} = {L_q} = {L_s} $，电磁转矩方程可简化为：

$ {T_e} = \frac{3}{2}{P_n}{i_q}{\psi _f} = {K_t}{i_q}，$

(3)

式中，$ {K_t} $为转矩常数。

电机的机械运动方程为：

$ J\frac{{{\rm{d}}{\omega _m}}}{{{\rm{d}}t}} = {T_e} - {T_L} - B{\omega _m}。$

(4)

式中：$ {\omega _m} $为电机的机械角速度；$ J $为电机转子惯量和负载转动惯量之和；$ {T_L} $为负载转矩；$ B $为黏性阻尼系数。

2 基于梯度校正法与扩展卡尔曼滤波器的惯量辨识 2.1 惯量辨识算法原理

基于梯度校正法与扩展卡尔曼滤波器的惯量辨识框图分为两部分，如图1所示。第一部分是船舶永磁同步电机推进系统矢量控制模块，图中$ \omega _m^ * $为速度给定，速度环采用$ {\text{PI}} $控制器，输出为电流给定$ i_q^ * $，电流环同样采用$ {\text{PI}} $控制器，其输出电压$ u_q^ * $和$ u_d^ * $经过$ {\text{IPark}} $变换，经$ {\text{SVPWM}} $发波，三相电流$ {i_a} $、$ {i_b} $、$ {i_c} $经$ {\text{Clarke}} $变换和$ {\text{Park}} $变换，得到电流反馈$ {i_q} $和$ {i_d} $，速度信号$ {\omega _m} $由编码器反馈。

第二部分为惯量辨识模块，由电机运动方程可知，转动惯量和负载转矩存在强耦合关系。因此把梯度校正法与扩展卡尔曼滤波器结合对转动惯量进行辨识，速度反馈$ {\omega _m} $和电流反馈$ {i_q} $作$ {\text{EKF}} $的输入，同样也作为梯度校正法的输入。因$ {\text{EKF}} $的系数矩阵需要实际的转动惯量值，当转动惯量变化时，会导致$ {\text{EKF}} $观测的负载转矩波动，因此需要把梯度校正法辨识的惯量值实时输入到$ {\text{EKF}} $中。传统算法认为负载转矩的变化周期大于系统的采样周期，假定负载不变或变化缓慢，但是本文考虑存在负载转矩扰动的情况，因此需要把$ {\text{EKF}} $观测的负载转矩实时输入到梯度校正法中，从而在有负载扰动工况下能够辨识出准确的惯量值。

2.2 梯度校正法基本原理

梯度校正法的基本思想是：沿准则函数的负梯度方向，逐步修正模型参数估计值，直至准则函数达到最小值，对于$ {\text{PMSM}} $这样非线性、强耦合、变参数的对象进行参数辨识很适宜。

设梯度校正法中的准则函数为：

$ J(\theta ) = \frac{1}{2}{\left[ {y(k) - {\varphi ^{\rm{T}}}(k)\widehat \theta (k - 1)} \right]^2} 。$

(5)

沿$ J(\theta ) $的负梯度方向，逐步修正辨识值$ \widehat \theta (k) $，直到$ J(\widehat \theta ) $达到最小。梯度校正法的递推形式为：

$ \widehat \theta (k) = \widehat \theta (k - 1) - {\boldsymbol{R}}({\boldsymbol{k}})grad{\left. {[J(\theta )]} \right|_{\widehat \theta (k)}} 。$

(6)

式中，${\boldsymbol{R}}({\boldsymbol{k}})$为权矩阵。

$ {\rm{grad}}{\left. {[J(\theta )]} \right|_{\widehat \theta (k)}} = - [y(k) - {\varphi ^{\rm{T}}}(k)\widehat \theta (k - 1)]\varphi (k) ，$

(7)

得到梯度校正法的递推公式：

$ \widehat \theta (k) = \widehat \theta (k - 1) + R(k)\varphi (k)[y(k) - {\varphi ^{\rm{T}}}(k)\widehat \theta (k - 1)] 。$

(8)

式中：$ \varphi (k) $为系统输入；$ y(k) $为系统输出。

根据$ {\text{Lyapunov}} $稳定性定理，$ R(k) $一般选择为：

$ R(k)=\frac{1}{{\displaystyle \sum _{i=1}^{N}{\Lambda }_{i}(k){\phi }_{i}^{2}(k)}}\cdot {\rm{diag}}[{\Lambda }_{1}(k)\text{，}{\Lambda }_{2}(k)\text{，}\cdots \text{，}{\Lambda }_{N}(k)]，$

(9)

式中，一般选择$ {\Lambda _i}(k) $为单位矩阵$ I $。因此式（8）可表示为：

$ \widehat \theta (k) = \widehat \theta (k - 1) + \frac{{\varphi (k)}}{{{\varphi ^{\rm{T}}}(k)\varphi (k)}} \cdot [y(k) - {\varphi ^{\rm{T}}}(k)\widehat \theta (k - 1)] 。$

(10)

为了避免$ {\varphi ^T}(k)\varphi (k) = 0 $，对梯度校正法的递推公式修正为：

$ \widehat \theta (k) = \widehat \theta (k - 1) + \frac{{\alpha \varphi (k)}}{{\lambda + {\varphi ^T}(k)\varphi (k)}} \cdot [y(k) - {\varphi ^T}(k)\widehat \theta (k - 1)]。$

(11)

式中：$ \alpha $和$ \lambda $为可调参数，参数范围为$ 0 \leqslant \alpha \leqslant 2 $，$ \lambda > 0 $。

对式（4）进行离散化，并忽略黏性阻尼系数$ B $，可得差分方程为：

$ {\omega _m}(k) - {\omega _m}(k - 1) = \frac{{{T_s}}}{J} \cdot ({T_e}(k) - {T_L}(k))，$

(12)

根据式（12）将$ k $时刻和$ k - 1 $时刻的机械运动方程相减可得：

$ \begin{aligned} &{\omega _m}(k) - 2{\omega _m}(k - 1) + {\omega _m}(k - 2) = \frac{{{T_s}}}{J} \cdot \\ &[{T_e}(k) - {T_e}(k - 1) - {T_L}(k) + {T_L}(k - 1)]。\end{aligned}$

(13)

式中：$ {T_s} $为算法采样时间；定义$ {\text{y}}(k) = {\omega _m}(k) - 2{\omega _m}(k - 1) +$$ {\omega _m}(k - 2) $；${\phi }^{{\rm{T}}}（k）={T}_{e}(k)-$$ {T_e}(k - 1) - {T_L}(k) +{T_L}(k - 1) $；$ \theta = \dfrac{{{T_s}}}{J} $。

将式（13）代入式（11），可得：

$\begin{aligned} & \widehat \theta (k) = \widehat \theta (k - 1) + \frac{{\alpha \varphi (k)}}{{\lambda + {\varphi ^{\rm{T}}}(k)\varphi (k)}} \cdot \\ & [{\omega _m}(k) - 2{\omega _m}(k - 1) + {\omega _m}(k - 2) - {\varphi ^{\rm{T}}}(k)\widehat \theta (k - 1)] 。\end{aligned}$

(14)

式中：$ \widehat \theta (k) = \frac{{{T_s}}}{J} $，根据$ \widehat \theta (k) $可求出转动惯量$ J $。

2.3 扩展卡尔曼滤波器设计

负载转矩扰动影响速度的稳定性，因此需要对负载转矩实时观测，并把负载转矩作为电流前馈补偿。因船舶推进系统是非线性的，$ {\text{EKF}} $适用于非线性系统，从而选择$ {\text{EKF}} $进行观测负载转矩。

电机负载转矩的变化率对于船舶推进系统的采样率来说变化缓慢，可认为导数为0，即

$ \frac{{{\rm{d}}{T_L}}}{{{\rm{d}}t}} = 0。$

(15)

根据式（4）和式（15），构建PMSM的状态方程为：

$ \left[\begin{array}{c}\stackrel{·}{{\omega }_{m}}\\ \stackrel{·}{{T}_{L}}\end{array}\right]=\left[\begin{array}{cc}-\frac{B}{J}& -\frac{1}{J}\\ 0& 0\end{array}\right]\left[\begin{array}{c}\omega {}_{m}\\ {T}_{L}\end{array}\right]+\left[\begin{array}{c}\dfrac{{K}_{t}}{J}\\ 0\end{array}\right]{i}_{q} 。$

(16)

系统的状态方程和测量方程描述如下：

$ \Bigg\{\begin{array}{l}\stackrel{·}x{(t)}=f[x(t)]+Bu(t)+w(t)，\\ y(t)=h[x(t)]+v(t)。\end{array} $

(17)

式中：$f[x(t)] = \left[ \begin{gathered} {{ - }}\frac{{{B}}}{{{J}}}{{{\omega }}_{{m}}}{{ - }}\frac{{{{{T}}_{{L}}}}}{{{J}}} \\ {{0}} \\ \end{gathered} \right]$；$B = \left[ \begin{gathered} \frac{{{{{K}}_{{t}}}}}{{{J}}} \\ {{0}} \\ \end{gathered} \right]$；$h[x(t)] = \left[ {{{\mathbf{\omega }}_{{m}}}} \right]$。

对$ f[x(t)] $和$ h[x(t)] $线性化得到雅可比矩阵：

$ F[x(t)] = {\left. {\dfrac{{\delta f}}{{\delta x}}} \right|_{x = x(t)}} = \left[ {\begin{array}{*{20}{c}} {{{ - }}\dfrac{{{B}}}{{{J}}}}&{{{ - }}\dfrac{{{1}}}{{{J}}}} \\ {{0}}&{{0}} \end{array}} \right]，$

(18)

$ H[x(t)] = {\left. {\frac{{\delta h}}{{\delta x}}} \right|_{x = x(t)}} = \left[ {\begin{array}{*{20}{c}} {{1}}&{{0}} \end{array}} \right]。$

(19)

式（17）为非线性方程，设采样周期为$ {T_s} $，用前向差分对方程离散化为：

$ \left\{ \begin{gathered} x(k + 1) = {A_k}x(k) + {B_k}u(k) + w(k)，\\ y(k) = {H_k}x(k) + v(k) 。\\ \end{gathered} \right. $

(20)

式中：${A_k} = \left[ {\begin{array}{*{20}{c}} {{{1 - }}\frac{{{B}}}{{{J}}}{{{T}}_{{s}}}}&{{{ - }}\frac{{{1}}}{{{J}}}{{{T}}_{{s}}}} \\ {{0}}&{{1}} \end{array}} \right]$；${B_k} = \left[ {\begin{array}{*{20}{c}} \frac{{{{{K}}_{{t}}}}}{{{J}}}{{{T}}_{{s}}} \\ {{0}} \\ \end{array}} \right]$；$ {H_k} $$= \left[ {\mathbf{1}}{{0}} \right]$；$x(k) = {\left[ {\begin{array}{*{20}{c}} {{{{\omega }}_{{m}}}{{(k)}}}&{{{{T}}_{{L}}}{{(k)}}} \end{array}} \right]^{\rm{T}}} $；$u(k) = \left[ {{{{i}}_{{q}}}{{(k)}}} \right]$；$y(k) = $$\left[ {{{{\omega }}_{{m}}}{{(k)}}} \right]$；$ w(k) $为系统高斯噪声，$w(k) = {\left[ {\begin{array}{*{20}{c}} {{{{w}}_{{\omega }}}}&{{{{w}}_{{{{T}}_{{L}}}}}} \end{array}} \right]^{\rm{T}}}$，$ {w_\omega } $和$ {w_{{T_L}}} $之间没有相关性；$ v(k) $为测量高斯噪声，$v(k) = {\left[ {{{{v}}_{{w}}}} \right]^{\rm{T}}}$，定义$ w $的协方差矩阵${\boldsymbol{ Q}}$以及$ v $的协方差矩阵$ R $为：

$ \left\{ \begin{gathered} \cos (w) = E\left\{ {w{w^{\rm{T}}}} \right\} = Q = \left[ {\begin{array}{*{20}{c}} {{{{q}}_{{w}}}}&{{0}} \\ {{0}}&{{{{q}}_{{{{T}}_{{L}}}}}} \end{array}} \right]，\\ \cos (v) = E\left\{ {v{v^{\rm{T}}}} \right\} = R 。\qquad\qquad\qquad\quad\\ \end{gathered} \right. $

(21)

式中，$ {q_w} $和$ {q_{{T_L}}} $分别为$ {w_\omega } $和$ {w_{{T_L}}} $的方差。

利用$ {\text{EKF}} $对式（16）所描述的状态变量进行实时观测，具体如下：

步骤1　状态预测

$ \widetilde x(k + 1) = \widehat x(k) + {T_s}[f(\widehat x(k)) + B(k)u(k)]。$

(22)

步骤2　先验估计协方差矩阵

$ \widetilde p(k + 1) = F(k)\widehat p(k + 1)F{(k)^{\rm{T}}} + Q。$

(23)

步骤3　增益矩阵

$ K(k + 1) = \widetilde p(k + 1){H^T}{[H\widetilde p(k + 1){H^T} + R]^{ - 1}}。$

(24)

步骤4　后验估计

$ \widehat x(k + 1) = \widetilde x(k + 1) + K(k + 1)[y(k + 1) - \widetilde y(k + 1)]。$

(25)

步骤5　后验估计协方差矩阵

$ \widehat p(k + 1) = \widetilde p(k + 1) - K(k + 1)H\widetilde p(k + 1)。$

(26)

3 仿真分析

为分析梯度校正法和$ {\text{EKF}} $的有效性，在软件$ {\text{Matlab/Simulink}} $中搭建船舶推进系统的仿真模型，仿真中永磁同步电机的参数如表1所示。

由电机的机械运动方程可知，当电机的速度变化率${\rm{d}}{\omega _m}/{\rm{d}}t$发生变化时辨识惯量$ J $才有意义。在仿真中速度给定为幅值在$- 500 \sim 500\;\;{\rm{r}}/\min$之间交替变化的三角波，其上升周期为$0.02\;{\text{s}}$。

3.1 定惯量辨识仿真

为测试在定惯量的情况下可调参数$ \alpha $和$ \lambda $对辨识结果的影响，设$ \alpha $不变，改变$ \lambda $的值，或者$ \lambda $不变，改变$ \alpha $的值，比较辨识效果。电机空载运行时，系统总的转动惯量为电机转子惯量，即$J = 0.559 \times {10^{ - 4}}\;{\rm{kg}} \cdot {{\rm{m}}^2}$。

保持$ \alpha = 0.5 $不变，调节$ \lambda $的大小，仿真结果如图2所示。$ \lambda $在$ 0.1 \sim 0.5 $之间变化时，随着$ \lambda $越大，收敛的速度越慢，但辨识精度越高；$ \lambda $越小，收敛时间快，但波动较大，精度变低，在$ \lambda $过大或过小时，会出现结果发散的情况。

选择$ \lambda = 0.1 $不变，改变$ \alpha $的大小，仿真结果如图3所示。$ \alpha $在$ 0.1 \sim 1 $的范围内，$ \alpha $的值越大，辨识的波动越大，误差也越大，收敛速度随着$ \alpha $的增大越来越快。

通过分析曲线的辨识精度和收敛时间，当$ \alpha = 0.5 $、$ \lambda = 0.1 $时，稳定时间在$ 0.3\;{\text{s}} $，辨识稳定值为$5.331{\rm{e}} - 5$，稳态误差为−4.5%。

3.2 变惯量辨识仿真

在$ {\text{Matlab/Simulink}} $软件中，根据PMSM的数学方程自定义创建电机模型，目的是把转动惯量$ J $的端口引出，从而可以接任何信号，达到变惯量仿真的效果。为测试算法在转动惯量发生突变时的性能，开始电机空载运行，设定在$t = 0.5\;{\text{s}}$时，突然加入 $ 1 $ 倍和 $ 4 $ 倍电机转子惯量的负载惯量，则系统总的转动惯量${J'_1} = J + {J_{L1}} = 1.118 \times {10^{ - 4}}\;{\rm{kg \cdot {m^2}}}$和$ {J'_2} = J + {J_{L2}} = 2.795 \times $${10^{ - 4}}\;{\rm{kg \cdot {m^2}}}。$

从图4和图5可知，随着转动惯量变化越大，收敛时间变长，辨识精度变低，误差值也变大，可以看出负载惯量从$ J $增加到$ 4J $，误差增加了0.5%。

3.3 负载扰动仿真

对存在扰动的工况，即电机转动惯量和负载转矩同时变化，当电机空载运行，在$t = 0.5\;{\rm{s}}$突加$ 1 $倍转动惯量，负载转矩突加至$2\;{\rm{N \cdot m }}$。分别采用传统梯度校正法和本文算法对转动惯量和负载转矩进行辨识，辨识结果如图6所示。

可以看出，传统梯度校正法在负载突变时，会导致转动惯量明显突变上升到很大的值。相反，本文算法在突加负载时没有突变，逐渐跟随上理论辨识值，稳定后辨识误差也小于传统算法。本文算法改善了动态辨识效果，提高算法的自适应能力。

根据调试经验，将EKF的协方差矩阵设为$Q= {\rm{diag}}(0.1,0.01）$，$ R = 0.1 $，观测值如图7所示。

可知，在空载工况下，观测值很快就能收敛，误差较小；在带载工况下，观测值需经过短时间的波动才收敛到真实值，误差值在$ 0.5\% $左右。

4 实验分析

为验证算法的有效性，进行在线惯量辨识实验，其电机转子惯量约为$2.51{\rm{e}} - 5\;{\rm{kg}} \cdot {{\rm{m}}^2}$，通过离线惯量辨识的方法测得惯量盘的惯量为$3.13{\rm{e}} -$$4\;{\rm{kg}} \cdot {{\rm{m}}^2}$。

分别在空载和带负载的工况下测试算法的辨识性能，考虑速度和转矩电流的周期比辨识周期小，经过实验调试，辨识周期取$10\;{\text{ms}}$，为便于观察，将$ J/{T_s} $作为理论观测值。取速度给定信号为三角波，上升时间为$0.2\;{\rm{s}}$，幅值在$- 1000 \sim 1000\;\;{\rm{r}}/\min$之间交替变化。

在空载条件下，转动惯量理论辨识值为J/T_s=$0.00251\;\;{\rm{s}}/({\rm{kg}} \cdot {{\rm{m}}^2})$，实验中，可调参数$ \alpha $和$ \lambda $可参考仿真参数，但仿真环境比较理想化，不一定适用实验，需要细微调整，最终选取$ \alpha = 2 $，$ \lambda = 0.005 $。图8为空载条件下的辨识结果，转动惯量在$ 5 $ s内可以收敛，其辨识结果为$0.00239\text{ }{\rm{s/(kg\cdot {m}^{2}}})$，与理论值的误差为$ - 4.78\% $。

图9为带负载的工况，转动惯量理论辨识值为$J/{T_s} = 0.034{\text{ }}{\rm{s/(kg \cdot {m^2})}}$转动惯量经过4.5 s达到稳定值，此时转动惯量的稳态值为$0.032{\text{ }}{\rm{s/(kg \cdot {m^2})}}$,误差为−5.04%。

为验证$ {\text{EKF}} $的性能，通过对拖电机施加负载转矩。首先电机为空载运行，在$3\;{\rm{s}}$左右，突加$0.18\;{\rm{N \cdot m}}$的负载转矩。如图10所示，可以稳定的观测负载转矩，响应快速，在变负载转矩的瞬间会有较大的冲击，但会迅速跟随负载转矩的变化，具有良好的自适应能力。

5 结　语

本文针对船舶永磁同步电机转动惯量和负载转矩存在的扰动问题，提出梯度校正法与扩展卡尔曼滤波器结合的在线惯量辨识方法，很好解决了转动惯量和负载转矩之间因耦合问题导致惯量辨识误差大的问题。通过对梯度校正法选择合适的可调参数，对于定惯量和变惯量的工况，提高了算法的辨识效果。相比传统算法，本文算法在负载转矩阶跃工况下有较好的自适应能力。此外，利用扩展卡尔曼滤波器对负载转矩实时观测，并作为电流前馈补偿。提高了船舶推进系统对负载转矩的抗扰动性能。仿真和实验验证了算法的有效性。