0 引　言

随着现代交流技术的飞速发展，电力推进型船舶逐渐走入航运市场。由于永磁材料工艺、控制技术、电力电子技术的发展，永磁同步电机（Permanent Magnet Synchronous Motor，PMSM）在电力推进型船舶中的应用也愈发广泛^[1]。伴随着永磁同步电机应用范围的不断扩大，其调速控制策略也在研究和探索的过程中不断提出。其中，矢量控制是公认最有效的永磁同步电机控制策略之一^[2]。而就电力推进型船舶本身而言，其拥有传统柴油机推进型船舶所无法比拟的优点^[3]。因此为了更好地发展电力推进型船舶，如何建立有效的船舶推进电机运动控制系统的仿真模型具有重要意义^[4]。

为了建立合适的船舶推进电机运动控制系统仿真模型，国内外学者开展了相关研究。文献[5-6]对螺旋桨负载特性方面进行仿真研究，但缺乏对电机控制方面的研究；文献[7]基于直接转矩控制理论分析船舶交流电力推进系统，但并未给出整体调速系统模型；文献[8]采用永磁同步电机对螺旋桨负载特性进行仿真模拟，却未给出永磁同步电机的具体相关参数；文献[9]虽然设计了船舶电力推进仿真系统，但并未给出实验结果。除上述以外，国内外在对电力推进型船舶控制系统进行建模仿真时，绝大部分是基于Matlab/Simulink软件平台下搭建相关模型。而由于在Simulink中很多模块是现成封装好的，内部模型不公开，不便于程序调试，同时与其他算法程序接口匹配困难。

针对上述问题，本文采用Matlab脚本语言编程方式，基于牛顿拉夫逊法求解矢量控制下差分化的永磁同步电机方程组，通过搭建船机桨一体数学模型，分析船机桨在不同工况下的动态响应。本文采用全代码方式的仿真也使得整个过程可视易调，可以十分便捷地增加控制算法，仿真结果也验证了该模型和方法的正确性。

1 船舶电力推进系统建模

当选用低速电机作为船舶推进电机时，减速机构可被省略，螺旋桨直接通过联轴器与电机相连，因此两者转速一致。本文以永磁同步电机为推进电机，船桨运动产生的转矩作为推进电机的负载输入转矩，其整体运动控制系统框图如图1所示。

1.1 永磁同步电机控制系统建模

由于永磁同步电机具有非线性、强耦合的特点，为简化其模型，通常依据矢量控制的思想，利用坐标变换的方式将三相固定坐标系下电机数学模型转换为两相旋转坐标系下的电机数学模型，从而使得其中部分变量解耦，达到类似直流电机调速系统的控制效果。依据图1所示搭建静止坐标系下的永磁同步电机模型，由于2次坐标变换会导致计算繁琐。为了简化后期控制器的设计，本文选择在同步坐标系dq轴下搭建PMSM数学模型。PMSM数学模型应当包括电压方程、磁链方程、运动方程以及电磁转矩方程，其相关表达式如下^[10]：

${u_d} = {R_{\rm{s}}}{i_d} + \dfrac{{{\rm{d}}{\psi _d}}}{{{\rm{d}}t}} - {\omega _e}{\psi _q}{\text{，}}$

(1)

${u_q} = {R_{\rm{s}}}{i_q} + \dfrac{{{\rm{d}}{\psi _q}}}{{{\rm{d}}t}}{\rm{ + }}{\omega _{e}}{\psi _{d}}{\text{，}}$

(2)

${\psi_d} = {L_d}{i_d} + {\psi _f}{\text{，}}$

(3)

${\psi_q} = {L_q}{i_q}{\text{，}}$

(4)

$\dfrac{{{\rm d}{\omega _m}}}{{{\rm d}t}} = \dfrac{1}{J}({T_e} - {T_L}){\text{，}}$

(5)

${T_{\rm{e}}} = {\rm{1}}{\rm{.5}}{p_n}[{\psi _f}{i_q} + ({L_d} - {L_q}){i_q}{i_d}]{\text{。}}$

(6)

式中：u_d，u_q为定子电压的dq轴分量；R_s为定子电阻；i_d，i_q为定子电流的dq轴分量；ψ_d，ψ_q为磁链的dq轴分量；ω_e为电机电角速度；ψ_f为永磁体磁链；L_d，L_q为直轴和交轴电感；J为转动惯量；ω_m为机械转速；T_L为负载转矩；T_e为电磁转矩；p_n为磁极对数。

1.2 逆变器建模

在电机矢量控制系统的电流环中，经过Park逆变换得到的u_α，u_β被送到SVPWM发生器，由此产生的PWM波经过驱动电路控制IGBT开关管的导通以及关断，从而控制电机的旋转。为了简化系统，方便分析，本文将逆变器部分看成增益为常数的电流环放大器。

1.3 船桨模型建模

船舶在航行过程中，依靠推进电机带动螺旋桨旋转而产生的推力，在克服海水对螺旋桨产生阻力矩的基础上，若剩余推力方向与船舶航行方向一致，则推动船舶加速前进，当剩余推力为0时，船速稳定。依据螺旋桨的基本原理^[11]可得直线航行时四象限船桨数学模型为：

$P = {K_p}\rho {n^2}D_P^4{\text{，}}$

(7)

${T_L} = {K_m}\rho {n^2}D_P^5{\text{，}}$

(8)

${J_P} = \dfrac{{{h_p}}}{{{D_p}}} = \dfrac{{{v_p}}}{{n{D_p}}}{\text{，}}$

(9)

$w = \dfrac{{{v_s} - {v_p}}}{{{v_s}}} = 1 - \dfrac{{{v_p}}}{{{v_s}}}{\text{，}}$

(10)

${v_p} = (1 - w){v_s}{\text{，}}$

(11)

$P = R + \Delta P{\text{，}}$

(12)

$R = P(1 - {t_s}){\text{，}}$

(13)

${P^{'}} = (1 - {t_s}){K_p}\rho {n^2}D_p^4 {\text{，}}$

(14)

$k{M_s}\dfrac{{{\rm d}{v_s}}}{{{\rm d}t}} = {P^{'}} - R{\text{。}}$

(15)

式中：P，T_L分别为螺旋桨旋转时产生的推力以及阻力矩；ρ为海水密度；n为螺旋桨转速；D_P为螺旋桨直径；K_p为螺旋桨推力的无因次系数；K_m为螺旋桨阻力矩的无因次系数；v_p为桨相对水的轴向速度；h_p为桨回转一周的的轴向进程，J_p为定义轴向进程与轴向速度的比值为进速比；w为伴流系数；v_p为螺旋桨进速；P为螺旋桨总推力；R为船舶所受的阻力；ΔP为推力减额；其中，若推力减额系数t_s已知的话，可根据式（13）计算阻力。

综上所述，简化以后的船舶推进电机运动控制系统如图2所示^[12]。

2 永磁同步电机电压方程数值求解 2.1 数值解法分析

由上述式（1）和式（2）可知，电机电压方程为一组常微分方程组。虽然常微分方程的解法有多种，但就电机这种具有非线性、强耦合特点的特殊系统来说，由于经常遇到操作及故障，因此一般不采用无法自起步、需要多于一步的状态量历史信息的多步法进行数值求解工作，而是选择可以自起步，即可以通过稳态值出发进行下一步数值求解工作的单步法。此外，在求解常微分方程组的过程当中，数值稳定性将是决定能否采用这种方法的重要考量因素。而所谓数值稳定性，是指在计算过程中前一步形成的误差在后续时步计算中是否收敛，若收敛，则数值稳定。

对设特征根为λ，计算步长为h的一阶微分方程进行不同单步法下的数值稳定性分析。计算真值与实际计算机得到的值之间有误差，将相邻两次误差之间比值称之为此方法的误差传递规律。在误差比值的绝对值小于1的判定条件下，得到每种方法的数值稳定域。

由表1计算可得知，梯形法的数值稳定域为整个左半平面。当λ<0时，即系统稳定，无论步长为何值，hλ均小于0，位于数值稳定域内，方法的计算误差都能收敛，结果不会发生畸变，。而当λ>0时，无论步长为何值，hλ均大于0，即系统无法稳定。这使得在实际计算中，不会误得“系统稳定”的结果^[13]。同理可知，其他几种单步法在数值稳定性方面均没有梯形法好。

2.2 数值求解

利用梯形法^[14]离散dq轴电机电压微分方程：

$\left\{\begin{aligned} & {i_d}(t) = \dfrac{h}{2}({f_1}({i_d}(t - 1),{i_q}(t - 1)) + {f_1}({i_d}(t),{i_q}(t))) + {i_d}(t - 1), \\ & {i_q}(t) = \dfrac{h}{2}({f_2}({i_d}(t - 1),{i_q}(t - 1)) + {f_2}({i_d}(t),{i_q}(t))) + {i_q}(t - 1); \end{aligned}\right.$

(16)

$\left\{\begin{aligned} & {f_1} = \dfrac{1}{{{L_d}}}{u_d}(t) + \dfrac{1}{{{L_d}}}{\omega _e}(t){i_q}(t){L_q} - \dfrac{1}{{{L_d}}}{R_s}{i_d}(t), \\ & {f_2} = \dfrac{1}{{{L_q}}}{u_q}(t) - \dfrac{1}{{{L_q}}}{\omega _e}(t){i_d}(t){L_d} - \dfrac{1}{{{L_q}}}{R_s}{i_q}(t) - \dfrac{1}{{{L_q}}}{\omega _e}(t){\psi _r}; \end{aligned}\right. $

(17)

$\small \left\{\begin{aligned} & {F_1}({i_d}(t),{i_q}(t)) = \left(\dfrac{{h{R_s}}}{{2{L_d}}} - 1\right){i_d}(t - 1) - \dfrac{{h{u_d}(t)}}{{{L_d}}} + \left(1 + \dfrac{{h{R_s}}}{{2{L_d}}}\right){i_d}(t) -\\ &\quad \dfrac{{h{L_q}{\omega _e}(t)({i_q}(t) + {i_q}(t - 1))}}{{2{L_d}}} = 0, \\ & {F_2}({i_d}(t),{i_q}(t)) = \dfrac{{h{L_d}{\omega _e}(t)({i_d}(t) + {i_d}(t - 1))}}{{2{L_q}}} + \dfrac{{h({\psi _r}{\omega _e}(t) - {u_q}(t))}}{{{L_q}}}+ \\ &\quad \left(1 + \dfrac{{h{R_s}}}{{2{L_q}}}\right){i_q}(t) + \left(\dfrac{{h{R_s}}}{{2{L_q}}} - 1\right){i_q}(t - 1) = 0\text{。} \end{aligned}\right. $

(18)

式（16）表示梯形法离散的表达式，式（17）表示式（16）中的f₁，f₂，式（18）则表示整体化简后的表达式。

利用牛顿-拉夫逊法差分后的方程为一组非线性方程，其求解步骤如下：

步骤1 　设定迭代次数变量k，而i_d（t）和i_q（t）为t时刻的未知变量。当牛顿迭代程序运行时，此时，迭代次数k=1，i_d（t-1）和i_q（t-1）为i^k_d（t）和i^k_q（t）的初值。

步骤2 　根据牛顿迭代法的公式计算求解非线性方程组Jacobean矩阵：

${{J}} = \left[ {\begin{array}{*{20}{c}} {{{\left. {\dfrac{{\partial {F_1}}}{{\partial {i_d}(t)}}} \right|}_k}}&{{{\left. {\dfrac{{\partial {F_1}}}{{\partial {i_q}(t)}}} \right|}_k}} \\ {{{\left. {\dfrac{{\partial {F_2}}}{{\partial {i_d}(t)}}} \right|}_k}}&{{{\left. {\dfrac{{\partial {F_2}}}{{\partial {i_q}(t)}}} \right|}_k}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {1 + \dfrac{{h{R_s}}}{{2{L_d}}}}&{ - \dfrac{{h{L_q}{\omega _e}(t)}}{{2{L_d}}}} \\ {\dfrac{{h{L_d}{\omega _e}(t)}}{{2{L_d}}}}&{1 + \dfrac{{h{R_s}}}{{2{L_q}}}} \end{array}} \right]{\text{。}}$

(19)

步骤3 　计算两次迭代电流之间的误差值。

步骤4 　若某次迭代电流的误差值小于一个设置的精度值ε，则中止迭代，输出上次迭代的电流值。

步骤5 　输出的电流计算公式如下：

$\left[ {\begin{array}{*{20}{c}} {i_d^{k + 1}} \\ {i_q^{k + 1}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\delta i_d^{k + 1}} \\ {\delta i_q^{k + 1}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {i_d^k} \\ {i_q^k} \end{array}} \right]{\text{。}}$

(20)

步骤6 　若迭代次数k大于预设的最大迭代值，迭代程序结束并输出i^k_d（t）和i^k_q（t），否则回到步骤2进行新的一轮迭代。

3 仿真实验分析 3.1 电机数值仿真

本文以中铁渤海铁路渡轮1号船舶^[15]为仿真对象，基于Matlab R2014a平台，搭建该船电机矢量控制系统仿真模型，对其进行不同数值解法的恒负载运行仿真实验,以此来模拟出电机从静止启动直到给定转速运行的状态。其中，仿真时间为5 s，给定恒负载为1000 N，电机步长为5*10^–5s给定转速为120 rad/s。实验所用电机相关参数如表2所示。

由图3可以看出，在dq轴下搭建的永磁同步电机模型中，q轴的电压电流在经历过PI调节后迅速达到恒定。与此同时，图3（d）中的转矩也达到恒定。由此可见，dq轴下电机搭建的模型与实际契合。图3（c）为各种算法下的电机转速曲线。其中，在求解电机差分方程时均为牛顿法，区别在于微分方程的离散方法不同。可以看出，虽然每种方法都能很快达到转速稳定，但是相对而言，前进欧拉法相较于其他几种方法波动较大。在起步时刻，梯形法响应速度较其他方法更快。虽然梯形法的超调相较于其他2种方法更大，但是仍属于正常范围内。在达到转速稳定期间，梯形法相较于其他方法优势很明显。因此，梯形法离散加牛顿法求解更适用于电机差分化后方程组。

3.2 船机桨负载仿真

由于K_p和K_m是在螺旋桨敞水试验中，通过曲线拟合得到其与进速比J_p之间的关系，表3为螺旋桨敞水特性曲线上不同的离散点。

利用Matlab内置函数polyfit进行2阶多项式曲线拟合，所得函数中的每项系数取4位有效数字。

${K_P} = 0. 389\;5 - 0.271\;2 J- 0.102\;6{J^2}{\text{，}}$

(21)

$10{K_m} = 0.495\;4 - 0.218\;3J - 0.209\;8{J^2}{\text{，}}$

(22)

文献[16]给出了中铁渤海铁路轮渡1号船舶在某几个船速下船舶受到阻力的详细计算结果，如表4所示。

与上述计算过程同理，得到拟合后的函数表达式如下：

$R = 0.201\;5 - 0.059\;05{V_s} + 0.017\;69V_s^2{\text{。}}$

(23)

至此，可利用上述算法以及船机桨数学模型来分析船舶在不同工况下的动态响应。

船舶启航操作一般可以分为直接启动和分级启动2种方式。

1）直接启动

船舶从零速状态下起车到达额定转速的的操作叫做正车。其正车时的船舶航速v_s、螺旋桨负载转矩T_L、螺旋桨转速v_p如图4所示。当电机开始工作时，螺旋桨转速和船舶航速迅速达到稳定，其中螺旋桨转速在大约5 s时开始达到稳定状态，如图4（c）所示。而螺旋桨转矩则经历了达到峰值后随即降低的过程，在500 s左右时达到稳定。图4（a）和图4（b）对比可知，当螺旋桨转矩上升时，船速也随之上升，而在螺旋桨转矩下降初期，船速仍处于上升期，直到螺旋桨转矩稳定时，此时船速恒定。由此可以得出结论，当船舶以额定转矩直接启动，在船舶加速的过程中，推进电机的转矩过载。

2）分级启动

为了避免推进电机转矩过载而导致的设备损坏，在实际船舶正车当中，一般选择分级启动。本文仿真实验，分级启动分为3个阶段，即通过阶段性地提高船舶转矩，使得船舶航速阶段性地稳定，从而达到预定航速，其船舶航速v_s、螺旋桨转矩T_L、螺旋桨转速v_p如图5所示。当船舶进行第1级启动时，此时船舶处于加速阶段，而船舶负载转矩处于较小阶段，船舶约在800 s时达到此级的稳定转速，随后进行下一级启动。在图5（b）可以看出，在每次级别启动的时刻螺旋桨转矩都有一次突变，随后又稳定在与螺旋桨转速相对应的转矩值上。在整个启动过程中，螺旋桨转矩的峰值约在170 kN，与直接启动中图4（b）相比，大大减小了推进电机的转矩，也减小了事故发生的概率或维护成本。同时，从图5（a）和图5（b）可以看出，分级启动的所用时间比直接启动所用时间更长。

4 结　语

本文在dq轴下搭建永磁同步电机矢量控制简化模型，通过船机桨之间的基本特性分析，结合螺旋桨敞水特性曲线以及实例船舶阻力计算，利用全代码的方式实现了船舶推进电机运动控制系统的仿真。电机数值仿真实验结果验证了牛顿加梯形法相较于其他方法，更适用于求解差分化的永磁同步电机电压微分方程组。在船舶推进电机运动控制系统仿真当中，不同工况下船机桨的动态响应也证明了该模型具有可行性。本文建立的数学模型及相应数值仿真解法可为后续研究提供理论基础。