Parameter Identification of Maritime Vessel Rudder PMSM Based on Extended Kalman Particle Filter Algorithm

Abstract

To address the issue of system parameter variations during the operation of a maritime light vessel rudder permanent magnet synchronous motor (PMSM), an extended Kalman particle filter (EKPF) algorithm that combines a particle filter (PF) with an extended Kalman filter (EKF) is proposed in this paper. This approach enables the online identification of motor resistance and inductance. For highly nonlinear problems that are challenging for traditional methods such as Kalman filtering, this algorithm is typically a statistical and effective estimation method that usually yields good results. Firstly, a standard linear discrete parameter identification model is established for a PMSM. Secondly, the PF algorithm based on Bayesian state estimation as a foundation for subsequent research is derived. Thirdly, the advantages and limitations of the PF algorithm are analyzed, addressing issues such as sample degeneracy, by integrating it with the Kalman filtering algorithm. Specifically, the EKPF algorithm for online parameter identification is employed. Finally, the identification model within MATLAB/Simulink is constructed and the simulation studies are executed to ascertain the viability of our suggested algorithm. The outcomes from these simulations indicate that the proposed EKPF algorithm identifies resistance and inductance values both swiftly and precisely, markedly boosting the robustness and enhancing the control efficacy of the PMSM.

1. Introduction

In modern maritime light vessel navigation systems, the accuracy and reliability of the autopilot system are crucial for navigational safety. As the core component of the autopilot system, the steering gear motor is required to possess high dynamic response, high efficiency, and precise control characteristics. Low-power permanent magnet synchronous motors (PMSMs) have been identified as the ideal choice for steering gear motors due to their excellent power density, efficiency, and control performance [1,2]. Over the years, with the advancement in PMSM drive control technology, scholars both domestically and internationally have proposed a plethora of control techniques to achieve the high-performance control of permanent magnet motors [3,4,5,6]. Whether it is speed control, position control, or torque control, achieving the high dynamic response and precise control goals requires accurate motor parameters. However, these parameters are significantly influenced by factors such as temperature, stator current, and magnetic flux saturation, leading to substantial variations. Consequently, methods such as vector control and direct torque control often fall short of achieving the results predicted by theoretical analysis. In the complex marine environment, the steering gear motor must maintain high-performance operation under various conditions. To realize efficient PMSM control, accurate parameter identification becomes crucial [7,8,9,10,11], leading to the emergence of various motor parameter identification methods. This is of great practical significance for enhancing the performance and reliability of the autopilot system. Depending on the operating conditions of the PMSM, the parameter identification can be categorized into offline parameter identification or online parameter identification. Offline parameter identification, due to its non-interference with the online operational status of the motor and conservation of the online control resources, has garnered widespread attention in the industrial sector. On the other hand, online parameter identification methods can be synchronized with the motor operation, enabling the real-time acquisition of the motor parameters. They meet the need for accurate parameter tracking under different motor operating conditions and reflect the changes in the motor status over time. Compared to offline parameter identification, the online methods offer better real-time performance and flexibility [12]. The traditional online identification methods include least squares estimation [13], model reference adaptive control [14,15], and Kalman filtering [16,17].

Considering that the inductance parameters are influenced by the PMSM operation status, a motor model utilizing the transient voltage equation is developed and a forgetting factor is incorporated to enhance the least squares identification method in [18], which can improve the algorithm’s tracking capabilities. When the forgetting factor is set to a certain value, the least squares method may struggle to ensure robustness. Ref. [8] treats the error between the theoretical and actual outputs as a variable and dynamically adjusts the forgetting factor during the identification process, thereby accelerating the convergence speed of the algorithm while ensuring good robustness. Assuming that the stator resistance and permanent magnet flux linkage are known and constant, Ref. [19] achieves the online identification of the d-axis and q-axis inductances through a model reference adaptive control algorithm. Ref. [20] proposes a disturbance compensation-based model reference adaptive system, designing a disturbance estimator to estimate the external disturbances in real time and updating the adaptation rate based on the disturbance, thus reducing the impact of parameter uncertainty and disturbances on the system and expanding the application scenarios of the model reference adaptive control algorithm. The Kalman filtering algorithm provides real-time updates of estimated values for linear systems, achieving optimal parameter estimation. However, in cases where the noise is unknown or the modeling error is significant, the Kalman filtering algorithm often fails to achieve the desired effect. Since the motor system is nonlinear, many scholars have proposed improvements to the Kalman filtering algorithm, such as extended Kalman filtering [21,22], unscented Kalman filtering [23], and cubature Kalman filtering [24], to extend its application scope. The prerequisite for using Kalman filtering is that the measurement noise follows a Gaussian distribution. However, in practice, the statistical properties of the system noise are often unknown or time-varying. In [25], an adaptive Kalman filtering algorithm based on variational Bayesian inference is presented, which selects appropriate conjugate distributions to estimate the covariance matrix of noise with minimal changes, to some extent addressing the aforementioned issues. However, it is limited by the linear Gaussian state model. The H_∞ filtering algorithm has good robustness to unknown noise but at the cost of accuracy. The H_∞ filtering algorithm is only suitable for linear systems, and, when applied to nonlinear systems, it still has many drawbacks. Therefore, modifications of the H_∞ filtering algorithm have been proposed. Xia proposes a fitted H_∞ filtering algorithm based on fitted transformation to address the problems in nonlinear uncertain systems, which shows good robustness [26,27]. Furthermore, outliers in the system can also reduce the estimation accuracy of the H_∞ filtering algorithm. Therefore, based on the H_∞ filtering algorithm, Zhao studies a series of derivative algorithms [28,29,30], among which the Krein space robust unscented Kalman filter based on the generalized maximum likelihood can effectively handle outliers [30].

In addition to the aforementioned common traditional identification methods, there are some identification methods that are highly effective for handling highly nonlinear systems but have not been widely applied due to certain limitations. For example, the PF algorithm, which suffers from low computational efficiency [31], has only been applied with the continuous improvement in computer computational capabilities in recent years. While the PF algorithm is relatively simple and can accurately estimate target parameters, its limitations include high computational complexity, particle degeneracy, and potential delays in responding to dynamic system changes. To address these issues, this paper employs the EKPF algorithm, which combines the extended Kalman filtering algorithm with the particle filtering algorithm, for parameter identification. This approach mitigates problems such as particle degeneracy to some extent and improves the overall performance of the particle filter.

2. Establishing the Parameter Identification Model for PMSM

The mathematical model of the PMSM is typically formulated in the synchronous rotating reference frame d-q. The stator voltage equation is provided by

$$\left\{\begin{array}{l}{u}_{\mathrm{d}}=R_{\mathrm{s}}{i}_{\mathrm{d}}+\frac{\mathrm{d}{\psi }_{\mathrm{d}}}{\mathrm{d}t}-{\omega }_{\mathrm{e}}{\psi }_{\mathrm{q}}\\ u_{\mathrm{q}}=R_{\mathrm{s}}{i}_{\mathrm{q}}+\frac{\mathrm{d}{\psi }_{\mathrm{q}}}{\mathrm{d}t}+{\omega }_{\mathrm{e}}{\psi }_{\mathrm{d}}\end{array}\right.$$

(1)

The stator flux linkage equation is

$$\left\{\begin{array}{l}{\psi }_{\mathrm{d}}=L_{\mathrm{d}}{i}_{\mathrm{d}}+{\psi }_f\\ {\psi }_{\mathrm{q}}=L_{\mathrm{q}}{i}_{\mathrm{q}}\end{array}\right.$$

(2)

Substituting Equation (2) into Equation (1):

$$\left\{\begin{array}{l}{u}_{\mathrm{d}}=R_{\mathrm{s}}{i}_{\mathrm{d}}+L_{\mathrm{d}}\frac{\mathrm{d}}{\mathrm{d}t}{i}_{\mathrm{d}}-{\omega }_{\mathrm{e}}{L}_{\mathrm{q}}{i}_{\mathrm{q}}\hfill \\ u_{\mathrm{q}}=R_{\mathrm{s}}{i}_{\mathrm{q}}+L_{\mathrm{q}}\frac{\mathrm{d}}{\mathrm{d}t}{i}_{\mathrm{q}}+{\omega }_{\mathrm{e}}\left(L_{\mathrm{d}}{i}_{\mathrm{d}}+{\psi }_f\right)\hfill \end{array}\right.$$

(3)

In the equations, u_d and u_q are the voltages along the d and q axes, respectively; i_d and i_q are the currents along the d and q axes, respectively; R_s is the stator resistance; ψ_d and ψ_q are the components of stator flux linkage along the d and q axes, respectively; ω_e is the electrical angular speed of the rotor; L_d and L_q are the d- and q-axis inductances, respectively; ψ_f is the flux linkage of the permanent magnet.

Selecting the components of current along the d and q axes, i_d and i_q, and the parameters to be identified, L_d, L_q, and R_s, as state variables, we establish the state-space equations for the permanent magnet synchronous motor. This paper primarily focuses on surface-mounted permanent magnet synchronous motors, where L_d = L_q = L_s. Therefore, after rearrangement, Equation (3) becomes

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{i}_{\mathrm{d}}\\ i_{\mathrm{q}}\\ R_{\mathrm{s}}\\ L_{\mathrm{s}}\end{array}\right]=\left[\begin{array}{cccc}-\frac{R_{\mathrm{s}}}{L_{\mathrm{s}}}& {\omega }_{\mathrm{e}}& 0& 0\\ -{\omega }_{\mathrm{e}}& -\frac{R_{\mathrm{s}}}{L_{\mathrm{s}}}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{d}}\\ i_{\mathrm{q}}\\ R_{\mathrm{s}}\\ L_{\mathrm{s}}\end{array}\right]+\left[\begin{array}{cc}\frac{1}{L_{\mathrm{s}}}& 0\\ 0& \frac{1}{L_{\mathrm{s}}}\\ 0& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{u}_{\mathrm{d}}\\ u_{\mathrm{q}}-{\omega }_{\mathrm{e}}{\psi }_f\end{array}\right]$$

(4)

Due to the presence of coupling terms in the coefficient matrix of state Equation (4), direct identification of R_s and L_s is rather complex. Therefore, we introduce intermediate variables a and b to simplify the identification equation. Let a = R_s/L_s and b = 1/L_s, Equation (4) becomes

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{i}_{\mathrm{d}}\\ i_{\mathrm{q}}\\ a\\ b\end{array}\right]=\left[\begin{array}{cccc}0& {\omega }_{\mathrm{e}}& -i_{\mathrm{d}}& u_{\mathrm{d}}\\ -{\omega }_{\mathrm{e}}& 0& -i_{\mathrm{q}}& u_{\mathrm{q}}-{\omega }_{\mathrm{e}}{\psi }_f\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{d}}\\ i_{\mathrm{q}}\\ a\\ b\end{array}\right]+w$$

(5)

The output y can be expressed as

$$y=\left[\begin{array}{llll}1\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 1\hfill & 0\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{d}}\\ i_{\mathrm{q}}\\ a\\ b\end{array}\right]+v$$

(6)

In the above two equations, w represents the process noise of the system, while v represents the measurement noise of the system. x = [i_d i_q a b]^T denotes the state variable matrix of the system, and y = [i_d i_q] represents the output variable matrix of the system. By discretizing Equations (5) and (6) with a sampling period T_s, we obtain the following equations:

$$\left\{\begin{array}{l}{x}_{\mathrm{k}+1}=F_{\mathrm{k}}{x}_{\mathrm{k}}+w_{\mathrm{k}}\\ y_{\mathrm{k}}=H_{\mathrm{k}}{x}_{\mathrm{k}}+v_{\mathrm{k}}\end{array}\right.$$

(7)

In the equations, w_k and v_k are mutually independent white noise processes. F_k and H_k are coefficient matrices defined as follows:

$$\left\{\begin{array}{l}{F}_{\mathrm{k}}=\left[\begin{array}{cccc}1& {\omega }_{\mathrm{e}}{T}_{\mathrm{s}}& -i_{\mathrm{d}}{T}_{\mathrm{s}}& u_{\mathrm{d}}{T}_{\mathrm{s}}\\ -{\omega }_{\mathrm{e}}{T}_{\mathrm{s}}& 1& -i_{\mathrm{q}}{T}_{\mathrm{s}}& \left(u_{\mathrm{q}}-{\omega }_{\mathrm{e}}{\psi }_f\right)T_{\mathrm{s}}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\\ H_{\mathrm{k}}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]\end{array}\right.$$

(8)

Using the filtering algorithm, parameters a and b can be identified through the above model. Then, by utilizing the relationship between a, b, R_s, and L_s, the identification values of resistance R_s and inductance L_s can be obtained. The specific identification process is illustrated in Figure 1, where the s-Function module represents the part where the filtering algorithm program is implemented.

3. Particle Filter

3.1. Bayesian State Estimation

Assume the equation describing a nonlinear system is as follows:

$$\begin{array}{l}{x}_{\mathrm{k}+1}=f_{\mathrm{k}}(x_{\mathrm{k}},w_{\mathrm{k}})\\ y_{\mathrm{k}}=h_{\mathrm{k}}(x_{\mathrm{k}},v_{\mathrm{k}})\end{array}$$

(9)

where k is the time index, x_k is the state variable, w_k is the process noise, y_k is the measurement value, and v_k is the measurement noise. The functions f_k(·) and h_k(·) are the process equation and measurement equation of the time-varying nonlinear system, respectively. The noises w_k and v_k are independent white noises with known probability density functions. The purpose of the Bayesian estimator is to estimate the probability density of the state x_k conditioned on the measurement values y₁, y₂, …, y_k. This conditional probability density is represented as

p(x_k|Y_k) = The probability density of the state x_k conditioned on the measurement values y₁, y₂, …, y_k.

(10)

At k = 1, the first measurement value is obtained. Therefore, the initial condition of the estimator is the probability density function of x₀ conditioned on the absence of measurement values, denoted as Y₀, expressed as

$$p\left(x_0\right)=p\left(x_0\right|Y_0)$$

(11)

Once p(x_k|Y_k) is computed, the most suitable estimate of x_k can be determined based on the specific problem. When the conditional probability density function p(x_k|Y_k) is multimodal, the mean of x_k cannot be used as an estimate.

To compute the conditional probability density function p(x_k|Y_k), a recursive method is required. Prior to this, it is essential to ascertain p(x_k|Y_k−1), representing the probability density of x_k conditioned on all measurements up to time k. This step is fundamental in deriving a recursive solution for evaluating p(x_k|Y_k). This can be expressed utilizing the Bayes’ theorem and the properties of joint probability density functions as follows:

$$p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)={\displaystyle \int p}\left[\left(x_{\mathrm{k}},x_{\mathrm{k}-1}\right)\mid Y_{\mathrm{k}-1}\right]d{x}_{\mathrm{k}-1}={\displaystyle \int p}\left[x_{\mathrm{k}}\mid \left(x_{\mathrm{k}-1},Y_{\mathrm{k}-1}\right)\right]p\left(x_{\mathrm{k}-1}\mid Y_{\mathrm{k}-1}\right)\mathrm{d}{x}_{\mathrm{k}-1}$$

(12)

From Equation (9), it can be observed that x_k is fully determined by x_k−1 and w_k−1. Therefore, p[x_k|(x_k−1,Y_k−1)] = p(x_k|x_k−1), and we have

$$p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)={\displaystyle \int p}\left(x_{\mathrm{k}}\mid x_{\mathrm{k}-1}\right)p\left(x_{\mathrm{k}-1}\mid Y_{\mathrm{k}-1}\right)\mathrm{d}{x}_{\mathrm{k}-1}$$

(13)

In the equation, the second probability density function on the right-hand side, p(x_k|x_k−1), is initially unknown but becomes known at the start (refer to Equation (11)). This function p(x_k|x_k−1) denotes the probability density of the state at time k, given the state at time (k − 1), which is established. Given our knowledge of the system equation f_k(·) and the noise w_k, this probability density function is ascertainable.

Now, consider the posterior conditional probability density function of x_k. Express this probability density function utilizing Bayes’ theorem and the properties of joint probability density functions as follows:

$$\begin{array}{c}p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}}\right)=\frac{p\left(Y_{\mathrm{k}}\mid x_{\mathrm{k}}\right)}{p\left(y_{\mathrm{k}}\right)}p\left(x_{\mathrm{k}}\right)=\frac{p\left[\left(y_{\mathrm{k}},Y_{\mathrm{k}-1}\right)\mid x_{\mathrm{k}}\right]}{p\left(y_{\mathrm{k}},Y_{\mathrm{k}-1}\right)}\underset{p\left(x_{\mathrm{k}}\right)}{\underbrace{\frac{p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)p\left(Y_{\mathrm{k}-1}\right)}{p\left(Y_{\mathrm{k}-1}\mid x_{\mathrm{k}}\right)}}}\\ =\frac{p\left(x_{\mathrm{k}},y_{\mathrm{k}},Y_{\mathrm{k}-1}\right)}{p\left(x_{\mathrm{k}}\right)p\left(y_{\mathrm{k}},Y_{\mathrm{k}-1}\right)}\frac{p\left(x_{\mathrm{k}},Y_{\mathrm{k}-1}\right)p\left(Y_{\mathrm{k}-1}\right)}{p\left(Y_{\mathrm{k}-1}\right)p\left(Y_{\mathrm{k}-1}\mid x_{\mathrm{k}}\right)}\end{array}$$

(14)

Multiply the numerator and denominator simultaneously by p(x_k, y_k):

$$p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}}\right)=\frac{p\left(x_{\mathrm{k}},y_{\mathrm{k}},Y_{\mathrm{k}-1}\right)p\left(x_{\mathrm{k}},Y_{\mathrm{k}-1}\right)p\left(Y_{\mathrm{k}-1}\right)}{p\left(x_{\mathrm{k}}\right)p\left(y_{\mathrm{k}},Y_{\mathrm{k}-1}\right)p\left(Y_{\mathrm{k}-1}\right)p\left(Y_{\mathrm{k}-1}\mid x_{\mathrm{k}}\right)}\frac{p\left(x_{\mathrm{k}},y_{\mathrm{k}}\right)}{p\left(x_{\mathrm{k}},y_{\mathrm{k}}\right)}$$

(15)

In Equation (15), the conditional probability density function is derived several times by utilizing the ratio of the joint probability density function to the marginal probability density function. This method effectively captures the relationship between the conditioned and conditioning events; we obtain

$$p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}}\right)=\frac{p\left[Y_{\mathrm{k}-1}\mid \left(x_{\mathrm{k}},y_{\mathrm{k}}\right)\right]p\left(y_{\mathrm{k}}\mid x_{\mathrm{k}}\right)p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)}{p\left(y_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)p\left(Y_{\mathrm{k}-1}\mid x_{\mathrm{k}}\right)}$$

(16)

Note that y_k is a function of x_k, so p[Y_k−1|(x_k, y_k)] = p(Y_k−1|x_k). After rearranging, we obtain

$$p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}}\right)=\frac{p\left(y_{\mathrm{k}}\mid x_{\mathrm{k}}\right)p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)}{p\left(y_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)}$$

(17)

All probability density functions on the right-hand side of the equation are known. From the measurement equation h_k(·) and the probability density function of the measurement noise v_k, we can obtain p(y_k|x_k). From Equation (13), we know p(x_k|Y_k−1). The probability density function of p(y_k|Y_k−1) can be obtained.

$$\begin{array}{cc}\hfill p\left(y_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)& ={\displaystyle \int p}\left[\left(y_{\mathrm{k}},x_{\mathrm{k}}\right)\mid Y_{\mathrm{k}-1}\right]\mathrm{d}{x}_{\mathrm{k}}\hfill \\ & ={\displaystyle \int p}\left[y_{\mathrm{k}}\mid \left(x_{\mathrm{k}},Y_{\mathrm{k}-1}\right)\right]p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)\mathrm{d}{x}_{\mathrm{k}}\hfill \end{array}$$

(18)

y_k is completely determined by x_k and v_k, so p[y_k|(x_k,Y_k−1)] = p(y_k|x_k), and

$$p\left(y_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)={\displaystyle \int p}\left[y_{\mathrm{k}}\mid x_{\mathrm{k}}\right]p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)\mathrm{d}{x}_{\mathrm{k}}$$

(19)

From the above discussion, we can conclude that P(y_k|x_k) can be obtained from the probability density function of the measurement equation h_k(·) and the measurement noise v_k, and p(x_k|Y_k−1) can be obtained from Equation (13).

In summary, the recursive equations for the Bayesian state estimator can be summarized as follows:

1. The system equation and measurement equation are as follows:

$$\begin{array}{l}{x}_{\mathrm{k}+1}=f_{\mathrm{k}}(x_{\mathrm{k}},w_{\mathrm{k}})\\ y_{\mathrm{k}}=h_{\mathrm{k}}(x_{\mathrm{k}},v_{\mathrm{k}})\end{array}$$

(20)

where w_k and v_k are independent white noise processes with known probability density functions.

2. Assuming the probability density function p(x₀) of the initial state is known, the estimator is initialized as follows:

$$p(x_0|Y_0)=p(x_0)$$

(21)

3. For k = 1, 2, …, execute the following equations:

(a) Obtain the prior probability density function from Equation (13) as follows:

$$p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)={\displaystyle \int p}\left(x_{\mathrm{k}}\mid x_{\mathrm{k}-1}\right)p\left(x_{\mathrm{k}-1}\mid Y_{\mathrm{k}-1}\right)\mathrm{d}{x}_{\mathrm{k}-1}$$

(22)

(b) Obtain the posterior probability density function from Equations (17) and (19) as follows:

$$p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}}\right)=\frac{p\left(y_{\mathrm{k}}\mid x_{\mathrm{k}}\right)p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)}{{\displaystyle \int p}\left(y_{\mathrm{k}}\mid x_{\mathrm{k}}\right)p\left(x_{\mathrm{k}}\mid Y_{\mathrm{k}-1}\right)\mathrm{d}{x}_{\mathrm{k}}}$$

(23)

3.2. Establishment of the Particle Filter

The PF is derived from Bayesian state estimation and is used for numerically implementing Bayesian filters. The concept is straightforward and intuitive. At the outset of the estimation process, we randomly generate a specific number N of state vectors according to the initial probability density function p(x₀). These state vectors are termed particles and are labeled as x⁺_0,i (where i = 1, …, N). For each subsequent time step k = 1, 2, …, the next batch of particles is produced using state equation f(·).

$$x_{\mathrm{k},i}^{-}=f_{\mathrm{k}-1}\left(x_{\mathrm{k}-1,i}^{+},w_{\mathrm{k}-1}^i\right)\hspace{1em}(i=1,\cdots ,N)$$

(24)

where each noise vector wⁱ_k−1 is randomly generated based on the known probability density function of w_k−1. After obtaining the measurement at time k, the likelihood probability density of each particle x⁻_k,i, denoted as p(y_k|x⁻_k,i), is computed. For example, if given an m-dimensional measurement equation y_k = h(x_k) + v_k, and v_k ~ N(0, R), the likelihood probability density q_i for the measurement y* under the condition that x_k equals the particle x⁻_k,i can be computed as follows:

$$\begin{array}{cc}\hfill q_i& =p\left[\left(y_{\mathrm{k}}=y^{*}\right)\mid \left(x_{\mathrm{k}}=x_{\mathrm{k},i}^{-}\right)\right]=p\left[v_{\mathrm{k}}=y^{*}-h_{\mathrm{k}}\left(x_{\mathrm{k},i}^{-}\right)\right]\hfill \\ & ~\frac{1}{{(2\pi )}^{m/2}|R{|}^{1/2}}\exp \left(\frac{-{\left[y^{*}-h\left(x_{\mathrm{k},i}^{-}\right)\right]}^T{\mathrm{R}}^{-1}\left[y^{*}-h\left(x_{\mathrm{k},i}^{-}\right)\right]}{2}\right)\hfill \end{array}$$

(25)

In the above equation, “~” indicates that this probability density function is proportional to the expression on the right-hand side. If this equation is applied to all particles x⁻_0,i (i = 1, …, N), then the “relative” likelihood probability density of the state being equal to each particle will be correct. Normalize the likelihood probability density obtained from Equation (25).

$$q_i=\frac{q_i}{{\displaystyle {\sum }_{j=1}^N{q}_j}}$$

(26)

This guarantees that the total of all likelihood probability densities sums to 1.

Next, we perform resampling. The resampling process entails discarding particles with minimal weights and concentrating on those with larger weights. By resampling, we amplify the presence of particles with substantial weights, thereby better representing the posterior distribution using particles and their associated weights. Resample particles from the computed likelihood probability density functions; i.e., randomly generate a new set of particles x⁺_k,i based on the likelihood probability density function q_i. There can be several different resampling methods; one of the most direct methods is as follows. For 1, …, N, perform the following two steps:

1. Generate a random number r from a uniform distribution in the interval [0, 1].

2. Accumulate the likelihood function q_i until the cumulative sum is greater than r. That is, for ${\displaystyle {\sum }_{m=1}^{j-1}{q}_m}<r$, but ${\displaystyle {\sum }_{m=1}^j{q}_m}\ge r$, the new particle x⁺_k,i is set equal to the old particle x⁻_k,j.

As the sample size N tends to ∞, the set probability density function of the new particles x⁺_k,i approximates p(x_k|y_k).

The steps of resampling can be summarized as follows:

$$x_{\mathrm{k},i}^{+}= \mathrm{Probability} \mathrm{of} x_{\mathrm{k},j}^{-} \mathrm{in}\mathrm{terms}\mathrm{of}{q}_j$$

(27)

The resampling process diagram is as Figure 2.

The primary bottleneck of the PF is often computational complexity. To address this, more efficient resampling techniques can be used. With particles x⁺_k,i now distributed according to p(x_k|y_k), arbitrary statistical properties of this function can be computed.

In summary, the steps of particle filter can be summarized as follows:

1. The system equation and measurement equation are as follows:

$$\begin{array}{l}{x}_{\mathrm{k}+1}=f_{\mathrm{k}}(x_{\mathrm{k}},w_{\mathrm{k}})\\ y_{\mathrm{k}}=h_{\mathrm{k}}(x_{\mathrm{k}},v_{\mathrm{k}})\end{array}$$

(28)

where w_k and v_k are independent white noise processes with known probability density functions.

2. Assuming the probability density function p(x₀) of the initial state is known, N initial particles x⁺_0,i(i = 1, …, N) are randomly generated based on p(x₀). The parameter N serves as a trade-off between computational load and estimation accuracy.

3. For k = 1, 2, …, execute the following steps:

(a) Execute time update using the known process equation and the probability density function of the process noise to obtain the prior particles x⁻_k,i:

$$x_{\mathrm{k},i}^{-}=f_{\mathrm{k}-1}\left(x_{\mathrm{k}-1,i}^{+},w_{\mathrm{k}-1}^i\right)\hspace{1em}(i=1,\cdots ,N)$$

(29)

where each noise vector wⁱ_k−1 is randomly generated according to the known probability density function of w_k−1.

(b) Compute the likelihood probability density q_i for each particle x⁻_k,i given the measurement y_k. This can be obtained by estimating p(y_k|x⁻_k,i) using the probability density function of the nonlinear measurement equation and the measurement noise.

(c) Normalize the obtained likelihood probability density as follows:

$$q_i=\frac{q_i}{{\displaystyle {\sum }_{j=1}^N{q}_j}}$$

(30)

Now, the sum of all likelihood probability densities equals 1.

(d) Perform resampling, i.e., randomly generate a set of posterior particles x⁺_k,i based on the likelihood probability densities q_i.

(e) Now that we have a set of particles x⁺_k,i distributed according to the probability density function p(x_k|y_k), arbitrary statistical properties of this probability density function can be computed.

4. Extended Kalman Particle Filter

PF has many advantages, such as being a statistical and effective estimation method for highly nonlinear systems, often yielding good results. However, achieving good performance with PF comes with significantly increased computational complexity. Additionally, when the regions of high probability density for the state-space probability density functions p(y_k|x_k) and p(x_k|y_k−1) do not overlap, sample degeneracy occurs, causing all particles to converge to a single value. This issue can be mitigated by increasing the number of particles N, but this rapidly leads to a large computational burden and typically only delays sample degeneration. Sample degeneration and other factors may prevent PF from promptly reflecting system dynamics.

Improving the performance of the particle filter can be achieved through combination with other filters. In this study, we adopt the EKPF, which integrates the EKF with the PF to enhance its performance and mitigate the issue of sample degeneration. In this approach, at each measurement instant, the EKF iterates over the particles, and these measurements are used to resample the particles. This is akin to running N Kalman filters concurrently and performing one resampling step after each measurement. After obtaining x⁻_k,i as per Equation (24), x⁺_k,i can be iteratively obtained from x⁻_k,i using the iteration equations of the EKF [32].

$$\begin{array}{l}{P}_{\mathrm{k},i}^{-}=F_{\mathrm{k}-1,i}{P}_{\mathrm{k}-1,i}^{+}{F}_{\mathrm{k}-1,i}^T+Q_{\mathrm{k}-1}\\ K_{\mathrm{k},i}=P_{\mathrm{k},i}^{-}{H}_{\mathrm{k},i}^T{\left(H_{\mathrm{k},i}{P}_{\mathrm{k},i}^{-}{H}_{\mathrm{k},i}^T+R_{\mathrm{k}}\right)}^{-1}\\ x_{\mathrm{k},i}^{+}=x_{\mathrm{k},i}^{-}+K_{\mathrm{k},i}\left[y_{\mathrm{k}}-h\left(x_{\mathrm{k},i}^{-}\right)\right]\\ P_{\mathrm{k},i}^{+}=\left(I-K_{\mathrm{k},i}{H}_{\mathrm{k},i}\right)P_{\mathrm{k},i}^{-}\end{array}$$

(31)

K_k,i represents the Kalman gain for the i-th particle, and P⁻_k,i denotes the a priori estimate error covariance for the i-th particle. The matrices F and H, representing the partial derivatives, are defined as follows:

$$\begin{array}{l}{F}_{\mathrm{k}-1,i}={\left.\frac{\partial f}{\partial x}\right|}_{x=x_{\mathrm{k}-1,i}^{+}}\\ H_{\mathrm{k},i}={\left.\frac{\partial h}{\partial x}\right|}_{x=x_{\mathrm{k},i}^{-}}\end{array}$$

(32)

Next, following the procedure in Section 3.2, resample x⁺_k,i and its corresponding covariance matrix P⁺_k,i. This method, which updates prior particles x⁻_k,i based on the measurement at time k before resampling, helps to mitigate particle degeneracy.

In summary, the EKPF is as follows.

1. The system equation and measurement equation are as follows:

$$\begin{array}{l}{x}_{\mathrm{k}+1}=f_{\mathrm{k}}(x_{\mathrm{k}},w_{\mathrm{k}})\\ y_{\mathrm{k}}=h_{\mathrm{k}}(x_{\mathrm{k}},v_{\mathrm{k}})\end{array}$$

(33)

where w_k and v_k are independent white noise processes with known probability density functions.

2. Assuming the probability density function p(x₀) of the initial state is known, N initial particles x⁺_0,i are randomly generated based on p(x₀), with corresponding covariance matrices P⁺_0,i = P⁺₀(i = 1, …, N). The choice of N should balance computational complexity and estimation accuracy.

3. For sampling time k = 1, 2, …, follow these steps:

(a) Based on the process equation and the probability density function of the process noise, obtain the prior particles x⁻_k,i and covariance P⁻_k,i through the time update equation:

$$\begin{array}{l}{x}_{\mathrm{k},i}^{-}=f_{\mathrm{k}-1}\left(x_{\mathrm{k}-1,i}^{+},w_{\mathrm{k}-1}^i\right)\\ P_{\mathrm{k},i}^{-}=F_{\mathrm{k}-1,i}{P}_{\mathrm{k}-1,i}^{+}{F}_{\mathrm{k}-1,i}^T+Q_{\mathrm{k}-1}\\ F_{\mathrm{k}-1,i}={\left.\frac{\partial f}{\partial x}\right|}_{x=x_{\mathrm{k}-1,i}^{+}}\end{array}$$

(34)

Here, each noise vector wⁱ_k−1 is randomly generated based on the known probability density function of w_k−1.

(b) Update the posterior particles and their covariance matrices based on the prior particles and their covariance matrices:

$$\begin{array}{l}{H}_{\mathrm{k},i}={\left.\frac{\partial h}{\partial x}\right|}_{x=x_{\mathrm{k},i}^{-}}\\ K_{\mathrm{k},i}=P_{\mathrm{k},i}^{-}{H}_{\mathrm{k},i}^T{\left(H_{\mathrm{k},i}{P}_{\mathrm{k},i}^{-}{H}_{\mathrm{k},i}^T+R_{\mathrm{k}}\right)}^{-1}\\ x_{\mathrm{k},i}^{+}=x_{\mathrm{k},i}^{-}+K_{\mathrm{k},i}\left[y_{\mathrm{k}}-h\left(x_{\mathrm{k},i}^{-}\right)\right]\\ P_{\mathrm{k},i}^{+}=\left(I-K_{\mathrm{k},i}{H}_{\mathrm{k},i}\right)P_{\mathrm{k},i}^{-}\end{array}$$

(35)

(c) Compute the likelihood probability density q_i for each particle x⁺_k,i conditioned on the measurement y_k. The likelihood probability density p(y_k|x⁺_k,i) can be estimated based on the nonlinear measurement equation and the probability density function of the measurement noise.

(d) Normalize the probabilities obtained in the previous step:

$$q_i=\frac{q_i}{{\displaystyle {\sum }_{j=1}^N{q}_j}}$$

(36)

At this point, all likelihood probabilities sum to 1.

(e) Correct the posterior particles x⁺_k,i and their corresponding covariance matrices P⁺_k,i based on the probabilities q_i. This constitutes the resampling process.

(f) Now that we have a set of posterior particles x⁺_k,i and covariance matrices P⁺_k,i, we can compute statistical quantities based on these.

5. Simulation Analysis and Comparison

To confirm the feasibility of the proposed parameter identification algorithm, this chapter carries out relevant simulation verification. In the simulation process, a simulation model is built using MATLAB/Simulink, and the two methods mentioned above are implemented in the s-Function module of Simulink. The flowchart of the s-Function module is shown in Figure 3.

The s-Function module is implemented using a non-graphical method to create a dynamic system. In this module, the input variables needed for the motor control system are obtained in real time, the state variables are updated in real time, and, thus, the real-time identification of the motor parameters in the control model is achieved. The selected input variables are u = [i_d i_q ω_e u_d u_q], state variables are x = [i_d i_q a b] (where a = R_s/L_s and b = 1/L_s), and output variables are y = [i_d i_q]. The chosen number of particles N is 2000. The particle cloud is generated by providing an initial particle and standard deviation of the initial particle cloud using the “randn” function. The motor parameters used are shown in

Table 1. Motor parameters used in the present study.

Electromagnetic Parameters	Value	Unit
DC voltage	24	V
Stator resistance	0.84	Ω
d-axis inductance	3	mH
q-axis inductance	3	mH
Flux linkage	0.01	Wb
Number of pole pairs	4	-

.

In the first part of this section, simulation analysis is conducted on the PF algorithm to verify its performance and identify any limitations. In the second part, simulation analysis is performed on the EKPF algorithm. Firstly, steady-state condition parameter identification simulation is conducted to validate the effectiveness of the method. Subsequently, to assess the robustness of the EKPF algorithm, simulation analyses are performed under three different operating conditions: variations in motor load, alterations in stator resistance, and changes in stator inductance.

5.1. PF Simulation Analysis

In this section, the PF algorithm is simulated with the following parameter settings: initial particles are set to [0; 5; 208.696; 434.78], initial standard deviation of particle cloud is set to [2; 8; 30; 60], process noise covariance matrix Q is diag([0 0 0.9 1.18]), measurement noise covariance matrix R is diag([1 1]), motor load is set to 0.3 N·m, and motor speed is set to 900 rpm, as shown in Figure 4. Three scenarios are simulated: (1) resistance and inductance remain stable; (2) resistance changes abruptly from 0.84 Ω to 1.2 Ω; (3) inductance changes abruptly from 3 mH to 4 mH.

(1) Resistance and inductance remain stable. The simulation results are shown in Figure 5.

From Figure 5, it can be observed that, for stable parameters, the PF algorithm achieves fast and accurate identification of R_s and L_s.

(2) When the resistance changes abruptly from 0.84 Ω to 1.2 Ω, the simulation results are shown in Figure 6.

(3) The inductance abruptly changes from 3 mH to 4 mH, as shown in the simulation results in Figure 7.

From Figure 6 and Figure 7, it can be observed that, when the stator resistance or stator inductance undergoes a sudden change, the identification process fails to promptly reflect the dynamic changes in the system.

In summary, the PF algorithm can perform identification relatively quickly and accurately to some extent. However, due to the issue of particle degeneracy in the PF algorithm, it may fail to respond promptly to dynamic changes in the system.

5.2. EKPF Simulation Analysis

In this section, simulations are conducted to analyze the EKPF algorithm. The parameter settings for the identification algorithm are as follows: the initial particles are set to [0; 5; 208.696; 434.78], the standard deviation of the initial particle cloud is [2; 8; 30; 60], the process noise covariance matrix Q is diag([0 0 0.9 1.18]), the measurement noise covariance matrix R is diag([1 1]), and the estimation error covariance matrix P₀ is diag([0.01 0.1 3 3]). The performance and robustness of the algorithm are validated separately for steady-state and dynamic conditions.

5.2.1. Steady-State Performance

The motor operates under steady-state conditions with a load torque of 0.3 N·m and a speed of 900 rpm, as shown in Figure 8, while the stator resistance and stator inductance remain constant at 0.84 Ω and 3 mH, respectively. The simulation results are shown in Figure 9.

From Figure 9, it can be observed that the EKPF algorithm achieves parameter identification of resistance and inductance in a short time, with the final identification results deviating by around 3% from the actual values. This demonstrates the effectiveness of the parameter identification algorithm.

5.2.2. Robustness Verification

(a) Load Torque

The motor speed is maintained at 900 rpm, while the torque is adjusted from 0.2 N·m to 0.4 N·m within a span of 0.5 s. The outcomes of this simulation are displayed in Figure 10.

From Figure 11, it can be observed that the EKPF algorithm ensures parameter identification even when the torque changes (doubles). The difference between the identification results before and after the change is within 2%, demonstrating the robustness of the EKPF parameter identification algorithm to torque variations.

(b) Stator Resistance

In this section, simulations were performed under two conditions: abrupt changes in resistance due to motor faults and gradual increases in resistance from factors such as temperature rise. The stator resistance shifts suddenly from 0.84 Ω to 1.2 Ω and then increases gradually to 1.2 Ω. The motor speed is fixed at 900 rpm, with the motor load set at 0.3 N·m. The outcomes of these simulations are presented in Figure 12.

From Figure 12, it can be observed that the EKPF algorithm ensures parameter identification during resistance changes, and it maintains a satisfactory response speed during sudden resistance changes, demonstrating the robustness of the EKPF parameter identification algorithm to resistance variations.

(c) Stator inductance

In this section, simulations were conducted for the inductance change conditions corresponding to the previous section: sudden inductance change due to motor failure and gradual increase in inductance due to improper design. The stator inductance changed abruptly from 3 mH to 4 mH, and gradually increased to 4 mH. The motor speed was established at 900 rpm, and the motor load was configured at 0.3 N·m. The results of the simulation are depicted in Figure 13.

From Figure 13, it can be observed that the EKPF algorithm ensures parameter identification during inductance changes, and it maintains a satisfactory response speed during sudden inductance changes, demonstrating the robustness of the EKPF parameter identification algorithm to inductance variations.

6. Conclusions

PMSMs are widely used in various fields, such as vessel autopilots. To achieve efficient PMSM control, the accurate identification of motor parameters is particularly important. This paper combines the PF algorithm with the EKF algorithm to form the EKPF algorithm for the online identification of motor resistance and inductance. The theoretical analysis and simulation results show that, compared to the PF algorithm, the EKPF algorithm can identify results more quickly and accurately, and its robustness is significantly improved. When the parameters to be identified change, whether suddenly or gradually, the EKPF algorithm maintains good identification accuracy and fast convergence speed.