*Article* **Design and Implementation of Position Sensorless Field-Excited Flux-Switching Motor Drive Systems**

## **Tian-Hua Liu \*, Muhammad Syahril Mubarok and Yu-Hao Xu**

Department of Electrical Engineering, National Taiwan University of Science and Technology, Taipei 106, Taiwan; syahril.elmubarok@gmail.com (M.S.M.); eddiechesterbn@gmail.com (Y.-H.X.)

**\*** Correspondence: Liu@mail.ntust.edu

Received: 20 June 2020; Accepted: 15 July 2020; Published: 16 July 2020

**Abstract:** Field-excited flux-switching motor drive systems have become more and more popular due to their robustness and lack of need for a permanent magnet. Three different types of predictive controllers, including a single-step predictive speed controller, a multi-step predictive speed controller, and a predictive current controller are proposed for sensorless flux-switching motor drive systems in this paper. By using a 1 kHz high-frequency sinusoidal voltage injected into the field winding and by measuring the a-b-c armature currents in the stator, an estimated rotor position that is near ±2 electrical degrees is developed. To improve the dynamic responses of the field-excited flux-switching motor drive system, predictive controllers are employed. Experimental results demonstrate the proposed predictive controllers have better performance than PI controllers, including transient, load disturbance, and tracking responses. In addition, the adjustable speed range of the proposed drive system is from 4 r/min to 1500 r/min. A digital signal processor, TMS-320F-2808, is used as a control center to carry out the rotor position estimation and the predictive control algorithms. Measured results can validate the theoretical analysis to illustrate the practicability and correctness of the proposed method.

**Keywords:** field-excited flux-switching motor; high-frequency injection; predictive controller; digital signal processor

## **1. Introduction**

The flux-switching motor is a type of double-salient structured motor with two windings in the stator, including an armature winding and a field winding, which could be replaced by a permanent magnet. The flux-switching motor has several advantages, including a robust structure, a sinusoidal back-electromotive force (back-EMF) waveform, and a reasonable torque density [1,2]. Several researchers have investigated the design of different types of flux-switching motors [3–5]. Recently, field-excited flux-switching motors have become more and more popular due to there being no need for a permanent magnet and the good flux-weakening operational characteristics [6]. Many papers have investigated field-excited flux-switching motors. For example, Ullah et al. proposed a field-excited linear flux-switching motor [7]. Gaussens proposed an analytical method of air-gap modeling for a field-excited flux-switching motor, by which the flux linkage and torque were derived [8].

Several papers have studied the control of flux-switching motor drive systems. For instance, Zhao et al. employed a model predictive controller for a flux-switching drive system to decrease its torque and flux ripples and to enhance the quality of the drive system [9]. Moreover, Zhao et al. investigated the low copper-loss control of a field-excited flux-switching drive system to improve its adjustable speed range and efficiency [10]. Zhao et al. implemented vector control of a field-excited flux-switching drive to maximize its output torque by using particle swarm optimization [11]. Yang et al. studied flux-weakening control of a hybrid flux-switching motor to increase its output torque and

high-speed operating range [12]. Wu et al. investigated field-oriented control and direct torque control for a five-phase flux-switching motor drive system with a fault-tolerant capability to improve its dynamics during faulty conditions [13]. Nguyen et al. proposed rotor position sensorless control of a field-excited flux-switching motor drive system using a high-frequency square-wave voltage injecting method to replace an encoder [14]. Zhang proposed model reference adaptive control to improve the dynamics of a sensorless flux-switching motor drive system [15].

Several researchers have proposed sensorless methods for motor drives. For example, Fan et al. proposed sensorless control of a five-phase interior permanent magnet synchronous motor (IPMSM) based on a high-frequency sinusoidal voltage injection [16]. Wang et al. investigated a position-sensorless control method at low speed for PMSM based on high-frequency signal injection into a rotating reference frame [17]. Compared to previously published papers [9–17], this paper is the first to propose predictive controllers for a sensorless field-excited flux-switching motor drive system. To the authors' best knowledge, this idea is an original idea. This is the first time that a multi-step predictive speed controller is developed and compared to a single-step predictive speed controller for a flux-switching motor drive system. This is another main contribution of the paper. The implemented drive system could be applied for grass cutters and vacuum cleaners due to its robustness, lack of need for a permanent magnet, and high torque. In addition, the high torque ripple, large volume, and serious acoustic noise may not be the main issues for the applications of vacuum cleaners and grass cutters.

## **2. Flux-Switching Motor**

This paper investigates a 3-phase, 6-slot armature winding stator, 7-salient-tooth rotor, field-excited, flux-switching motor, which is shown in Figure 1. The armature winding is located inside the motor, and the salient-tooth rotor is located outside of the motor. In Figure 1, the A1, A2, B1, B2, C1, and C2 are armature windings, and F1, F2, F3, F4, F5, and F6 are field windings. All of them are located on the teeth of the stator. Between the stator and rotor, there is a non-uniform air gap.

**Figure 1.** The structure of a flux-switching motor.

To explain the basic principle of the flux-switching motor, Figure 2a–d illustrate the different rotor positions between the rotor and the stator, in which *S*1, *S*2, and *S*<sup>3</sup> are the slot numbers of the stator. Figure 3 shows the induced flux linkage and its related back-EMF. Figure 2a illustrates that the flux linkage of the F1 field winding rises, but the flux linkage of the F2 field winding decreases. As a result, the total flux linkage in the A2 armature winding is zero. This situation is also indicated at the a-point in Figure 3. Here the flux linkage of the A2 armature winding is zero, but the back-EMF of the A2 armature winding reaches its maximum value. Figure 2b illustrates that both the F1 and F2 field windings provide rising flux at the A2 armature winding. As a result, the total flux of the A2 armature winding reaches its maximum value, and the back-EMF of the A2 armature winding is zero, which is shown at the b-point in Figure 3. Figure 2c illustrates that the flux linkage of the field winding F1 also rises, but the flux linkage of the field winding F2 also decreases. As a result, the total flux linkage of the A2 armature winding is zero, and the back-EMF of the A2 armature winding reaches its negative minimum value, which is shown at the c-point of Figure 3. Figure 2d illustrates that both the flux linkage of the F1 and F2 field windings decrease. As a result, the total flux linkage of the A2 armature winding reaches its negative minimum value, and the back-EMF of the A2 armature winding is zero. This is illustrated as the d-point in Figure 3. In Figure 3, ω*re* is the electric rotor speed of the rotor and λ*<sup>m</sup>* is the flux linkage of the stator when rotor rotates. According to the above analysis, it is feasible to generate a flux linkage waveform and a back-EMF waveform, which are sinusoidal and cosine waveforms, respectively, and are illustrated in Figure 3. By using the back-EMF from the armature winding, torque is generated and the motor rotates smoothly.

**Figure 2.** Rotor flux at different rotor positions. (**a**) R is aligned with the A2 axis; (**b**) R is between the F1 and A2 axes; (**c**) R is aligned with the F1 axis; (**d**) R is on the left of the F1 axis.

**Figure 3.** Flux linkage and back-EMF waveforms.

### **3. Mathematical Model of a Flux-Switching Motor**

The synchronous frame d-q axis stator voltage of a flux-switching motor is expressed as follows:

$$
\begin{bmatrix} v\_d \\ v\_q \end{bmatrix} = \begin{bmatrix} r\_s + L\_d \frac{d}{dt} & -\omega\_{re} L\_q & L\_{df} \frac{d}{dt} \\\\ \omega\_{re} L\_d & r\_s + L\_q \frac{d}{dt} & \omega\_{re} L\_{qf} \end{bmatrix} \begin{bmatrix} i\_d \\ i\_q \\ i\_f \end{bmatrix} \tag{1}
$$

where *vd* and *vq* are the d- and q-axis voltages, *rs* is the rotor resistance, *Ld* and *Lq* are the d- and q-axis self-inductances, *Ld f* is the mutual inductance between the d-axis of the armature winding and the field winding, *Lq f* is the mutual inductance between the q-axis of the armature and the field winding, <sup>ω</sup>*re* is the electrical rotor speed, *<sup>d</sup> dt* is the differential operator, *id* and *iq* are the d- and q-axis currents, and *if* is the field current.

The flux linkage of the field winding, λ*<sup>f</sup>* , can be expressed as follows:

$$
\lambda\_f = \mathbf{L}\_{ff}\mathbf{i}\_f + \frac{\mathfrak{Z}}{2}\mathbf{L}\_{df}\mathbf{i}\_d \tag{2}
$$

where λ*<sup>f</sup>* is the total flux linkage of the field winding, *Lf f* is the self-inductance of the field winding, and *Ld f* is the mutual inductance between the field winding and the d-axis winding. The voltage of the field winding is expressed as follows:

$$
\mathbf{v}\_f = \mathbf{r}\_f \mathbf{i}\_f + \mathbf{L}\_{ff} \frac{d}{dt} \mathbf{i}\_f + \frac{3}{2} \mathbf{L}\_{df} \frac{d}{dt} \mathbf{i}\_d \tag{3}
$$

where *vf* is the voltage of the field winding, and *rf* is the resistance of the field winding. The total torque is shown as follows:

$$T\_e = \frac{3}{2} P\_m \left[ L\_{dqf} i\_f i\_q + (L\_d - L\_q) i\_d i\_q \right] \tag{4}$$

The mechanical speed of the motor is

$$\frac{d}{dt}\omega\_{rm} = \frac{1}{f\_{st}}(T\_{\varepsilon} - T\_{L} - B\_{\text{sf}}\omega\_{rm})\tag{5}$$

where ω*rm* is the mechanical speed of the rotor, *Jst* is the inertia of the motor and load, *Te* is the total output torque, *TL* is external load, and *Bst* is the total viscous friction coefficient. The differential of the mechanical position of the motor is

$$\frac{d}{dt}\theta\_{rm} = \omega\_{rm} \tag{6}$$

## **4. Rotor Position Estimator Design**

In this paper, assuming *Ldh* = *Lqh* = *Lsh*, the d-q-f axis high-frequency voltages and high- frequency currents can be described as follows:

$$
\begin{bmatrix} \upsilon\_{dh} \\ \upsilon\_{qh} \\ \upsilon\_{fh} \end{bmatrix} = jf\_h \begin{bmatrix} L\_{sh} & 0 & L\_{msf} \\ 0 & L\_{sh} & 0 \\ \frac{3}{2}L\_{msf} & 0 & L\_{ffh} \end{bmatrix} \begin{bmatrix} i\_{dh} \\ i\_{qh} \\ i\_{fh} \end{bmatrix} \tag{7}
$$

where *vdh*, *vqh*, and *vf h* are the high-frequency d-axis, q-axis, and field winding voltages, *Lsh* is the self-inductance of the d- and q-axis, *Lms f* is the mutual inductance between the d-axis or q-axis and the field winding inductance, and *Lf fh* is the high-frequency self-inductance of the field winding. From Figure 4, the coordinate transformation between d-q-f and α-β-f can be expressed as

$$
\begin{bmatrix} f\_d \\ f\_q \\ f\_f \end{bmatrix} = \begin{bmatrix} \cos\theta\_{\nu\varepsilon} & \sin\theta\_{\nu\varepsilon} & 0 \\ -\sin\theta\_{\nu\varepsilon} & \cos\theta\_{\nu\varepsilon} & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} f\_a \\ f\_\beta \\ f\_f \end{bmatrix} = T(\theta\_{\nu i}) \begin{bmatrix} f\_a \\ f\_\beta \\ f\_f \end{bmatrix} \tag{8}
$$

where *fd* and *fq* are the d- and q-axis voltages or currents, *f*<sup>α</sup> and *f*<sup>β</sup> are the α- and β-axis voltages or currents, and *T*(θ*re*) is the coordinate transformation matrix. Substituting Equation (8) into (7), one can obtain

$$
\begin{bmatrix} v\_{ah} \\ v\_{fh} \\ v\_{fh} \end{bmatrix} = j f\_h \begin{bmatrix} \theta\_{rh} \\ \theta\_{rh} \end{bmatrix}^{-1} \begin{bmatrix} L\_{sh} & 0 & L\_{dbf} \\ 0 & L\_{sh} & 0 \\ \frac{3}{2} L\_{dbf} & 0 & L\_{ffh} \end{bmatrix} T(\theta\_{lv}) \begin{bmatrix} i\_{ah} \\ i\_{fh} \\ i\_{fh} \end{bmatrix} \tag{9}
$$

$$
= j f\_h \begin{bmatrix} \cos \theta\_{lv} & \sin \theta\_{lv} & 0 \\ -\sin \theta\_{lv} & \cos \theta\_{lv} & 0 \\ 0 & 0 & 1 \end{bmatrix}^{-1} \begin{bmatrix} L\_{sh} & 0 & L\_{dbf} \\ 0 & L\_{sh} & 0 \\ \frac{3}{2} L\_{dbf} & 0 & L\_{ffh} \end{bmatrix} \begin{bmatrix} \cos \theta\_{lv} & \sin \theta\_{lv} & 0 \\ -\sin \theta\_{lv} & \cos \theta\_{lv} & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} i\_{ah} \\ i\_{fh} \\ i\_{fh} \end{bmatrix} \tag{9}
$$

where *v*α*h*, *v*β*h*, and *vf h* are the α − *axis*, β − *axis*, and field winding voltages, and *i*α*h*, *i*β*h*, and *if h* are the α − *axis*, β − *axis*, and field winding currents.

From Equation (9), after doing some mathematical processes, one can obtain

$$
\begin{bmatrix} v\_{ah} \\ v\_{\rho h} \\ v\_{fh} \end{bmatrix} = jf\_h \begin{bmatrix} L\_{\rm sh} & 0 & L\_{dqf} \cos \theta\_{re} \\ 0 & L\_{\rm sh} & L\_{dqf} \sin \theta\_{re} \\ \frac{3}{2} L\_{dqf} \cos \theta\_{re} & \frac{3}{2} L\_{dqf} \sin \theta\_{re} & L\_{ffl} \end{bmatrix} \begin{bmatrix} i\_{ah} \\ i\_{\rho h} \\ i\_{fh} \end{bmatrix} \tag{10}
$$

From Equation (10), one can derive the following equation:

$$\begin{bmatrix} i\_{\rm ah} \\ i\_{\rm fh} \\ i\_{\rm fh} \end{bmatrix} = \frac{1}{jf\_{\rm h}(L\_{\rm sh}^{2}L\_{\rm ffh} - \frac{3}{2}L\_{\rm dff}^{2}\sin\theta\_{\rm r}\tau^{2})} \begin{bmatrix} L\_{\rm ck}L\_{\rm ffh} - \frac{3}{2}L\_{\rm dff}^{2}\sin\theta\_{\rm r}\cos\theta\_{\rm r} & -\frac{3}{2}L\_{\rm ck}L\_{\rm dff}\cos\theta\_{\rm r} \\ 3L\_{\rm dff}\sin\theta\_{\rm r}\cos\theta\_{\rm r} & L\_{\rm sk}L\_{\rm ffh} - \frac{3}{2}L\_{\rm df}^{2}\cos\theta\_{\rm r}\tau^{2} - \frac{3}{2}L\_{\rm dk}L\_{\rm dff}\sin\theta\_{\rm r} \\ -\frac{3}{2}L\_{\rm sh}L\_{\rm dff}\cos\theta\_{\rm r} & -\frac{3}{2}L\_{\rm sh}L\_{\rm dff}\sin\theta\_{\rm r} & L\_{\rm sh}^{2} \end{bmatrix} \begin{bmatrix} v\_{\rm ah} \\ v\_{\rm fh} \\ v\_{\rm fh} \end{bmatrix} \tag{11}$$

In this paper, a high-frequency voltage *vf h*(*t*), which is injected into the field winding, is expressed as follows:

$$
\begin{bmatrix} v\_{ah} \\ v\_{fh} \\ v\_{fh} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ v\_{fh}(t) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ V\_{\text{inj}} \sin f\_{\text{ft}} t \end{bmatrix} \tag{12}
$$

Substituting Equation (12) into (11) and doing some mathematical processes, one can obtain the high-frequency α-axis and β-axis currents as follows:

$$
\begin{bmatrix} i\_{\text{alt}} \\ i\_{\text{fl}} \\ i\_{\text{fl}} \end{bmatrix} = \frac{V\_{\text{inj}}}{j f\_h (L\_{\text{sh}}^2 L\_{fft} - \frac{3}{2} L\_{\text{dq}f}^2)} \begin{bmatrix} -\frac{3}{2} L\_{\text{sl}} L\_{\text{dq}f} \cos \theta\_{\text{re}} + \frac{3}{2} L\_{\text{sl}} L\_{\text{dq}f} \cos \theta\_{\text{re}} \cos (2f\_{\text{fl}}t) \\\\ -\frac{3}{2} L\_{\text{sl}} L\_{\text{dq}f} \sin \theta\_{\text{re}} - \frac{3}{2} L\_{\text{sl}} L\_{\text{dq}f} \cos (2f\_{\text{fl}}t) \\\\ L\_{\text{sl}}^2 - L\_{\text{sl}}^2 \cos (2f\_{\text{fl}}t) \end{bmatrix} \tag{13}
$$

After using a low-pass filter to remove the high-frequency components, one can obtain the following equation:

$$
\begin{bmatrix} i\_{\text{sh}\\_p} \\ i\_{\text{fl}\\_p} \\ i\_{\text{fl}\\_p} \end{bmatrix} = \frac{V\_{\text{inj}}}{j f\_h (L\_{\text{sh}}^2 L\_{ffh} - \frac{3}{2} L\_{\text{dq}f}^2)} \begin{bmatrix} -\frac{3}{2} L\_{\text{sl}} L\_{\text{dq}f} \cos \theta\_{re} \\ -\frac{3}{2} L\_{\text{sl}} L\_{\text{dq}f} \sin \theta\_{re} \\ L\_{\text{sl}}^2 \end{bmatrix} \tag{14}$$

where *i*α*h*\_*p*, *i*β*h*\_*p*, and *if h*\_*<sup>p</sup>* are the high-frequency currents after a band-pass filtering. From Equation (14), the estimated rotor position, θ*re*, can be derived as follows:

$$\widehat{\boldsymbol{\theta}}\_{n\varepsilon} \cong \tan^{-1} \left| \frac{\boldsymbol{i}\_{\varepsilon h - \boldsymbol{p}} \boldsymbol{p}}{\boldsymbol{i}\_{\alpha h - \boldsymbol{p}} \boldsymbol{p}} \right| \cong \tan^{-1} \frac{\boldsymbol{(} \boldsymbol{V}\_{\mathrm{inj}} \frac{\boldsymbol{3}}{2} \boldsymbol{\mathcal{L}}\_{\mathrm{sl}f} \boldsymbol{\mathcal{L}}\_{\mathrm{d}qf} \boldsymbol{\hat{\theta}} \boldsymbol{\mathcal{L}}\_{\mathrm{d}qf} \boldsymbol{\mathcal{I}})}{\boldsymbol{(} \boldsymbol{V}\_{\mathrm{inj}} \frac{\boldsymbol{3}}{2} \boldsymbol{\mathcal{L}}\_{\mathrm{sl}\boldsymbol{h}} \boldsymbol{\mathcal{L}}\_{\mathrm{d}qf} \boldsymbol{\cos\theta} \boldsymbol{\mathcal{e}}\_{\mathrm{rc}}} \tag{15}$$

The estimated speed can be described as follows:

$$
\widehat{\omega}\_{\text{re}}(k) = \frac{\widehat{\theta}\_{\text{re}}(k) - \widehat{\theta}\_{\text{re}}(k-1)}{\Delta T} \tag{16}
$$

where ω*re*(*k*) is the estimated speed and Δ*T* is the sampling interval. In this paper, a speed state estimator is constructed [18]. A high-order speed estimator is employed to avoid the high-frequency noise caused by the difference operator. Figure 4 shows the relationship between the estimated rotor position θˆ*re* and the real position θ*re*. Figure 5 demonstrates the control block diagram of the sensorless drive system, which includes the current control of the armature winding, the current control of the field winding, the rotor position estimation, and the rotor speed control.

**Figure 4.** Estimated and real rotor positions.

**Figure 5.** Control block diagram of proposed sensorless drive system.

## **5. Predictive Controller Design**

Predictive control began in the late 1970s and has been developed significantly since then. It has been widely used in chemical processes, robotics, drying towers, and steam generators [19,20]. Recently, it has been used for power converters and motor drives [21–23]. The predictive speed-loop controllers include single-step and multi-step control inputs. The design methods are discussed as follows.

#### *5.1. One-Step Predictive Speed Controller*

By assuming that the external load *TL* is zero, the discrete dynamics of the speed can be expressed as

$$\begin{split} \omega\_{rm}(n+1) &= e^{-\frac{B\_{\rm st}}{\int\_{\rm st}}T\_{\rm s}} \omega\_{rm}(n) + \frac{1}{B\_{\rm st}} (1 - e^{-\frac{B\_{\rm st}}{\int\_{\rm st}}T\_{\rm s}}) \mathcal{K}\_{\rm st} i\_{\rm q}(n) \\ &= a\_{\rm rm} \omega\_{rm}(n) + b\_{\rm rm} i\_{\rm q}(n) \end{split} \tag{17}$$

and

$$a\_{rm} = e^{-\frac{B\_{st}}{J\_{st}}r\_s} \tag{18}$$

$$b\_{rm} = \frac{1}{B\_{st}} \left( 1 - e^{-\frac{B\_{st}}{J\_{st}}T\_s} \right) \mathbf{K}\_{st} \tag{19}$$

where *Kst* is the torque constant of the motor, *arm*, and *brm* are the discrete time parameters of the motor, and n is the step number. By defining *xrm*(*n*) = ω*rm*(*n*), and *u*(*n*) = *iq*(*n*), one can obtain

$$\mathbf{x}\_{rm}(n+1) = a\_{rm}\mathbf{x}\_{rm}(n) + b\_{rm}u(n) \tag{20}$$

The output *yrm*(*n*) is defined as

$$
\lambda y\_{rm}(n) = \mathfrak{x}\_{rm}(n) \tag{21}
$$

By using Equation (20) and taking one step back, one can obtain

$$
\alpha\_{rm}(n) = a\_{rm}\chi\_{rm}(n-1) + b\_{rm}u(n-1) \tag{22}
$$

By subtracting (20) from (22), one can obtain

$$\begin{aligned} \Delta \mathbf{x}\_{rm}(n+1) &= \mathbf{x}\_{rm}(n+1) - \mathbf{x}\_{rm}(n) \\ &= a\_{rm} \Delta \mathbf{x}\_{rm}(n) + b\_{rm} \Delta u(n) \end{aligned} \tag{23}$$

and

$$
\Delta \mathbf{x}\_{rm}(n) = \mathbf{x}\_{rm}(n) - \mathbf{x}\_{rm}(n-1) \tag{24}
$$

The input difference Δ*u*(*n*) is defined as

$$
\Delta u(n) = u(n) - u(n-1) \tag{25}
$$

The output Δ*yrm*(*n* + 1) is defined as follows:

$$
\Delta y\_{rm}(n+1) = y\_{rm}(n+1) - y\_{rm}(n) = \mathbf{x}\_{rm}(n+1) - \mathbf{x}\_{rm}(n) \tag{26}
$$

Then we can obtain

$$\begin{aligned} y\_{rm}(n+1) &= \Delta y\_{rm}(n+1) + y\_{rm}(n) \\ &= \Delta \mathbf{x}\_{rm}(n+1) + y\_{rm}(n) \\ &= a\_{rm} \Delta \mathbf{x}\_{rm}(n) + b\_{rm} \Delta u(n) + y\_{rm}(n) \end{aligned} \tag{27}$$

After that, we can define a cost function as follows [19,20]:

$$\begin{split} \left(\Omega\_{\text{sp}}(n)\right) &= \left(\omega\_{\text{rm}}^{\*}(n+1) - y\_{\text{sm}}(n+1)\right)^{2} + q\Delta u^{2}(n) \\ &= \left(\Delta \omega\_{\text{rm}}(n+1)\right)^{2} + q\left(\Delta u\_{\text{n}}(n)\right)^{2} \end{split} \tag{28}$$

where q is the weighting factor between the (Δ*un*(*n*))<sup>2</sup> and the (Δω*rm*(*n* + 1))<sup>2</sup> . Then, by taking the differential of the Ω*sp*(*n*) to the Δ*u*(*n*) and setting its result to be zero, one can obtain

$$
\Delta u(n) = \frac{b\_{rm}\omega\_{rm}^\*(n+1) - a\_{rm}b\_{rm}\Delta\omega\_{rm}(n) - b\_{rm}\omega\_{rm}(n)}{b\_{rm}^2 + q} \tag{29}
$$

The Δ*iq*(*n*) is used to replace the Δ*u*(*n*), and then Equation (29) can be rewritten as follows:

$$\begin{split} \Delta i\_{\boldsymbol{q}}(\boldsymbol{n}) &= \frac{b\_{rm}(\boldsymbol{\omega}\_{\boldsymbol{rm}}^{\*}(\boldsymbol{n}+1) - \boldsymbol{\omega}\_{\boldsymbol{rm}}(\boldsymbol{n})) - a\_{rm}b\_{rm} \Delta \boldsymbol{\omega}\_{\boldsymbol{rm}}(\boldsymbol{n})}{b\_{rm}^{2} + q} \\ &= k\_{1}(\boldsymbol{\omega}\_{\boldsymbol{rm}}^{\*}(\boldsymbol{n}+1) - \boldsymbol{\omega}\_{\boldsymbol{rm}}(\boldsymbol{n})) - k\_{2} \Delta \boldsymbol{\omega}\_{\boldsymbol{rm}}(\boldsymbol{n}) \end{split} \tag{30}$$

and

$$k\_1 = \frac{b\_{rm}}{b\_{rm}^2 + q} \tag{31}$$

$$k\_2 = \frac{a\_{rm}b\_{rm}}{b\_{rm}^2 + q} \tag{32}$$

Finally, the q-axis current command is

$$
\dot{i}\_q^\*(n) = \dot{i}\_q^\*(n-1) + \Delta \dot{i}\_q^\*(n) \tag{33}
$$

From Equations (30)–(33), the control block diagram of the single-step predictive speed-loop controller is demonstrated in Figure 6.

**Figure 6.** Single-step predictive speed controller.

#### *5.2. Multi-Step Predictive Speed Controller*

In this paper, to implement the multi-step predictive speed controller, a two-step control horizon and a two-step prediction horizon are developed. The details are discussed as follows. The (*n* + 1) augmented model is

$$X\_{\rm sm}(n+1) = A\_{\rm sm} X\_{\rm sm}(n) + B\_{\rm sm} \Delta u(n) \tag{34}$$

and

$$A\_{sm} = \begin{bmatrix} a\_{rm} & 0 \\ a\_{rm} & 1 \end{bmatrix} \tag{35}$$

$$B\_{sm} = \begin{bmatrix} b\_{rm} \\ b\_{rm} \end{bmatrix} \tag{36}$$

The augmented state variable for the (*n* + 2) step is

$$\begin{aligned} \left(X\_{\rm sm}(n+2)\right) &= A\_{\rm sm}X\_{\rm sm}(n+1) + B\_{\rm sm}\Delta u(n+1) \\ &= A\_{\rm sm}^2 X\_{\rm sm}(n) + A\_{\rm sm}B\_{\rm sm}\Delta u(n) + \Delta u(n+1) \end{aligned} \tag{37}$$

In addition, the (*n* + 1) step predictive output is

$$y\_{\rm sm}(n+1) = \mathcal{C}\_{\rm sm} A\_{\rm sm} X\_{\rm sm}(n) + \mathcal{C}\_{\rm sm} B\_{\rm sm} \Delta u(n) \tag{38}$$

The (*n* + 2) step predictive output is

$$y\_{\rm sm}(\mathfrak{n}+\mathfrak{2}) = \mathbb{C}\_{\rm sm}A\_{\rm sm}^2X\_{\rm sm}(\mathfrak{n}) + \mathbb{C}\_{\rm sm}A\_{\rm sm}B\_{\rm sm}\Delta\mathfrak{u}(\mathfrak{n}) + \mathbb{C}\_{\rm sm}B\_{\rm sm}\Delta\mathfrak{u}(\mathfrak{n}+\mathfrak{1})\tag{39}$$

Then, a new predictive output vector can be defined as follows:

$$Y\_{sm} = \begin{bmatrix} y\_{sm}(n+1) \\ y\_{sm}(n+2) \end{bmatrix} \tag{40}$$

and

$$X\_{\rm sm} = \begin{bmatrix} \varkappa\_{\rm sm}(n+1) \\ \varkappa\_{\rm sm}(n+2) \end{bmatrix} \tag{41}$$

After that, one can obtain

$$Y\_{sm} = F\_{sm}X\_{sm} + \Theta\_{sm} \Delta U\_{sm} \tag{42}$$

and

$$F\_{sm} = \begin{bmatrix} \mathbb{C}\_{sm} A\_{sm} \\ \mathbb{C}\_{sm} A\_{sm}^2 \end{bmatrix} \tag{43}$$

$$
\Theta\_{sm} = \begin{bmatrix}
\mathcal{C}\_{sm} B\_{sm} & 0 \\
\mathcal{C}\_{sm} A\_{sm} B\_{sm} & \mathcal{C}\_{sm} B\_{sm}
\end{bmatrix} \tag{44}
$$

$$
\Delta l I\_{\rm sm} = \begin{bmatrix} \Delta u(n) \\ \Delta u(n+1) \end{bmatrix} \tag{45}
$$

Now, one can define the performance index as follows [17–19]:

$$\Psi\_{sp} = \left(R\_{sm}^\* - \Upsilon\_{sm}\right)^T \left(R\_{sm}^\* - \Upsilon\_{sm}\right) + \Delta t I\_{sm}^T Q \Delta t I\_{sm} \tag{46}$$

and

$$R\_{sm}^{\*} = \begin{bmatrix} \omega\_{rm}^{\*}(n+1) \\ \omega\_{rm}^{\*}(n+2) \end{bmatrix} \tag{47}$$

$$Q = \begin{bmatrix} q & 0 \\ 0 & q \end{bmatrix} \tag{48}$$

By substituting (42) into (46) and doing the ∂Ψ*sp* ∂Δ*Usm* processes, one can obtain

$$\frac{\partial \Psi\_{sp}}{\partial \Delta \mathcal{U}\_{\rm sm}} = -2\Theta\_{\rm sm}^{T}(\mathcal{R}\_{\rm sm}^{\*} - \mathcal{F}\_{\rm sm}X\_{\rm sm}) + 2\left(\Theta\_{\rm sm}^{T}\Theta\_{\rm sm} + Q\right)\Delta \mathcal{U}\_{\rm sm} = 0\tag{49}$$

After that, one can obtain the optimal Δ*Usm* as

$$
\Delta \mathcal{U}\_{\rm sm} = \left(\Theta\_{\rm sm}^T \Theta\_{\rm sm} + \mathcal{Q}\right)^{-1} \left(\Theta\_{\rm sm}^T R\_{\rm sm}^\* - \Theta\_{\rm sm}^T F\_{\rm sm} X\_{\rm sm}\right) \tag{50}
$$

By substituting (41), (43)–(44), and (47)–(48) into Equation (50), one can finally derive the following equations:

$$
\begin{bmatrix}
\Delta u(n) \\
\Delta u(n+1)
\end{bmatrix} = \begin{bmatrix}
\frac{\mathcal{G}rm}{f\_{rm}} \\
\frac{h\_{rm}}{f\_{rm}}
\end{bmatrix} \tag{51}
$$

and

$$f\_{rm} = b\_{rm}^4 + q^2 + b\_{rm}^2 q (a\_{rm}^2 + 2a\_{rm} + 3) \tag{52}$$

$$\begin{aligned} g\_{rm} &= b\_{rm}^3 \left(\omega\_{rm}^\*(n+1) - \omega\_{rm}(n)\right) + a\_{rm} b\_{rm} q(\omega\_{rm}^\*(n+2) - \omega\_{rm}(n)) \\ &+ b\_{rm} q(\omega\_{rm}^\*(n+2) - 2\omega\_{rm}(n)) - a\_{rm} b\_{rm}^3 \Delta \omega\_{rm}(n) \\ &- b\_{rm} q \Delta \omega\_{rm}(n) \left(a\_{rm}^3 + 2a\_{rm}^2 + 2a\_{rm}\right) \end{aligned} \tag{53}$$

$$\begin{aligned} b\_{rm} &= b\_{rm}^3 \left( \omega\_{rm}^\*(n+2) - a\_{rm} \omega\_{rm}^\*(n+1) - \omega\_{rm}^\*(n+1) + a\_{rm} \omega\_{rm}(n) \right) \\ &+ b\_{rm} q \left( \omega\_{rm}^\*(n+2) - \omega\_{rm}(n) \right) - a\_{rm} b\_{rm} q \Delta \omega\_{rm}(n) \left( a\_{rm} + 1 \right) \end{aligned} \tag{54}$$

The control output in this paper can be rearranged as follows:

$$
\begin{bmatrix}
\Delta i\_q(n) \\
\Delta i\_q(n+1)
\end{bmatrix} = \begin{bmatrix}
\frac{\mathcal{G}rm}{f\_{rm}} \\
\frac{h\_{rm}}{f\_{rm}}
\end{bmatrix} \tag{55}
$$

$$\mathbf{i}\_q^\*(n) = \mathbf{i}\_q^\*(n-1) + \Delta \mathbf{i}\_q^\*(n) \tag{56}$$

$$i\_q^\*(n+1) = i\_q^\*(n) + \Delta i\_q^\*(n+1)\tag{57}$$

Combining (56) and (57), one can provide the control input q-axis current command as follows:

$$i\_{q\mathcal{C}}^\*(n) = \rho i\_q^\*(n) + (1-\rho)i\_q^\*(n+1) \tag{58}$$

where ρ is the weighting factor of the control input power. From Equations (51)–(58), one can obtain the multi-step predictive speed controller, which is displayed in Figure 7.

**Figure 7.** Multi-step predictive speed controller.

## *5.3. Predictive Current Loop Controller*

From Equation (1), one can easily derive the following d-q axis current equations:

$$\frac{d}{dt}\dot{\mathbf{i}}\_d = \frac{1}{L\_d}(v\_d - r\_s\dot{\mathbf{i}}\_d + \omega\_{re}\mathbf{L}\_q\dot{\mathbf{i}}\_q) \tag{59}$$

and

$$\frac{d}{dt}\dot{\mathbf{t}}\_q = \frac{1}{L\_q} (v\_q - r\_s \dot{\mathbf{t}}\_q - \omega\_{re} (\mathbf{L}\_d \dot{\mathbf{t}}\_d + \boldsymbol{\lambda}\_m)) \tag{60}$$

From Equations (59) and (60), and by inserting a zero-order-hold and taking the discrete form, one can obtain

$$i\_d(n+1) = e^{-\frac{T\_s}{L\_d}T\_{cu}}i\_d(n) + \frac{1 - e^{-\frac{T\_s}{L\_d}T\_{cu}}}{r\_s} \left[v\_d(n) + \omega\_{\nu\epsilon}(n)L\_q i\_q(n)\right] \tag{61}$$

and

$$i\_q(n+1) = e^{-\frac{T\_s}{L\_q}T\_{cu}}i\_q(n) + \frac{1-e^{-\frac{T\_s}{L\_q}T\_{cu}}}{r\_s} \left[v\_q(n) - \omega\_{re}(n)(L\_di\_d(k) + \lambda\_m)\right] \tag{62}$$

$$u\_d(n) = v\_d(n) + \omega\_{\text{re}}(n) L\_q i\_q(n) \tag{63}$$

and

$$u\_q(n) = v\_q(n) - \omega\_{\text{re}}(n)(L\_d i\_d(n) + \lambda\_m) \tag{64}$$

where *ud*(*n*) and *uq*(*n*) are the d- and q-axis current loop control inputs, and by finally substituting (63)–(64) into (61)–(62), one can obtain

$$i\_d(n+1) = a\_{cd}i\_d(n) + b\_{cd}
u\_d(n)\tag{65}$$

and

$$i\_q(n+1) = a\_{cq}i\_q(n) + b\_{cq}
u\_q(n)\tag{66}$$

In Equations (65) and (66), *acd*, *bcd*, *acq*, and *bcq*, which are the parameters of the IPMSM, can be expressed as follows:

$$a\_{cd} = e^{-\frac{T\_s}{L\_d}T\_c} \tag{67}$$

$$b\_{cd} = \frac{1 - e^{-\frac{T\_s}{L\_d}T\_c}}{r\_s} \tag{68}$$

$$a\_{cq} = e^{-\frac{T\_s}{L\_q}T\_c} \tag{69}$$

and

$$b\_{c\eta} = \frac{1 - e^{-\frac{T\_s}{L\_\eta}T\_c}}{r\_s} \tag{70}$$

By using similar processes [22], one can obtain

$$\boldsymbol{v}\_{d}"\left(\boldsymbol{n}\right) = \frac{\boldsymbol{z}}{z-1} \bigg( \frac{\boldsymbol{b}\_{cd}\boldsymbol{q}\left(\operatorname{id}^{\*}\left(\boldsymbol{n}+1\right)-\operatorname{i}\_{d}\left(\boldsymbol{n}\right)\right)}{\boldsymbol{b}\_{cd}^{2}\boldsymbol{q}+\boldsymbol{r}} - \frac{\boldsymbol{a}\_{cd}\boldsymbol{b}\_{cd}\boldsymbol{q}\Delta\boldsymbol{i}\_{d}\left(\boldsymbol{n}\right)}{\boldsymbol{b}\_{cd}^{2}\boldsymbol{q}+\boldsymbol{r}}\bigg) - \boldsymbol{\omega}\_{\text{rc}}\left(\boldsymbol{n}\right)\Big(\boldsymbol{L}\_{\boldsymbol{q}}\boldsymbol{i}\_{\boldsymbol{q}}\left(\boldsymbol{n}\right)\Big) \tag{71}$$

and

$$w\_q^\*(n) = \frac{z}{z-1} \bigg| \frac{b\_{cq}q\left(i\_q^\*(n+1) - i\_q(n)\right)}{b\_{q}^2q + r} - \frac{a\_{cq}b\_{cq}q\Delta i\_q(n)}{b\_{cq}^2q + r} \bigg| + a\_{\nu\varepsilon}(n)(L\_di\_d(n) + \lambda\_{\text{m}}) \tag{72}$$

Then, from Equations (71) and (72), the block diagram of the proposed predictive current control can be obtained as shown in Figure 8.

**Figure 8.** Current loop predictive controller.

## **6. Implementation**

To evaluate the correctness and feasibility of the proposed method, a field-excited flux-switching sensorless flux-switching motor drive system is implemented. The implemented block diagram of the drive system is displayed in Figure 9a, which includes a voltage-source inverter, a field-excited flux-switching motor, an H-bridge circuit that controls the excited field winding, a digital signal processor, Hall-effect sensors, and A/D converters. The switching frequency of the inverter is 10 kHz, and the switching frequency of the H-bridge circuit is 20 kHz. In addition, the sampling time of the current control is 100 μ*s* and the sampling time of the speed control is 1 *ms*. The DC voltages of the inverter and H-bridge are both 250 V. The digital signal processor (DSP) is manufactured by Texas Instruments, type TMS-320-F2808 [23]. To obtain a closed-loop high performance drive system, the DSP reads the a-phase current, b-phase current, and field winding current via A/D converters. After that, the DSP executes the rotor position estimation and all of the predictive control algorithms, and it also determines the inverter and H-bridge triggering signals.

A photograph of the implemented hardware circuit is displayed in Figure 9b, including a DSP, two drivers, two Hall-effect current sensors, an inverter, and an H-bridge circuit. The insulated gate bipolar transistor (IGBT) modules are made by the Mitsubishi Company, type CM200dy-12NF. The photo-couple gate drivers are manufactured by Avago Company, type HCPL-3120. The Hall-effect sensors are manufactured by the Swiss LEM company, type LA25-NP with a 100 kHz bandwidth. The A/D converters with fast conversion time are made by Analog Devices, type AD7655. Figure 9c displays the photograph of the flux-switching motor, which is connected to the dynamometer, in which a DC generator is used as an external load. The field-excited flux-switching motor is a seven stator salient teeth motor, with rated specifications of a power of 540 W, a speed of 600 r/min, and a current of 6 A. Its parameters include a stator resistance of 1.3 Ω, a d-axis inductance of 18.9 mH, a q-axis

inductance of 23 mH, an inertia of 0.0143 kg-m2, and a friction coefficient of 0.0047 N.m.s/rad. The motor was assembled in our laboratory because it is not available on the market.

(**c**)

**Figure 9.** Photographs of the implementation system. (**a**) Block diagram, (**b**) hardware circuits, and (**c**) motor and dynamometer.

The details of the control algorithms are shown in Figure 10a–d, which are the flowcharts of the DSP. Figure 10a is the flowchart of the main program, which waits for the zero-voltage vector to detect currents. Figure 10b is the speed-loop interrupt service routine, which uses Equation (16) to compute ωˆ*rm* and Equations (30) and (33) to determine *i* ∗ *<sup>q</sup>*. Figure 10c is the current loop interrupt service routine, which uses Equation (8) to execute the coordinate transformation, the band-pass filtering to obtain *i*α*h*, and finally computes θˆ*re*. Figure 10d is the flowchart of the H-bridge field current control. The principle is measuring the current and comparing it with the field current command. Then, proportional-integral (PI) controller is used as the current controller to obtain the field voltage. Finally, the high-frequency voltage is injected and pulse-width modulation (PWM) is executed.

**Figure 10.** *Cont*.

**Figure 10.** *Cont*.

**Figure 10.** Digital signal processor (DSP) flowcharts. (**a**) Main program, (**b**) current loop interrupt service routine, (**c**) speed-loop interrupt service routine, and (**d**) H-bridge field current control.

## **7. Experimental Results**

To verify the proposed methods, our experimental results are illustrated in this section. The input DC voltage of the inverter and the H-bridge are both 250 V. The switching frequency of the inverter is 10 kHz, and the switching frequency of the H-bridge is 20 kHz. The sampling time of the a-b-c axis current control is 100 μ*s*, and the sampling time of the speed control is 1 *ms*. The selection of the sampling intervals is based on the required computation time of the current loop and the speed loop. In addition, a synchronous relationship, in which the current loop executes 10 times and then the speed loop executes one time, is employed in this paper. The injection high-frequency voltage of the field winding is ±25*V*, which is near 10% of the DC bus voltage of the H-bridge.

Figure 11a demonstrates the measured a-phase current at 300 r/min and a 2 N.m load. The current is influenced by the different switching frequencies of the inverter and field winding. The measured a-phase current is a near sinusoidal waveform with high-frequency harmonics. The b-phase and c-phase are similar to the a-phase with a 120 degree phase shift. Figure 11b demonstrates the measured current of the field winding, which is clearly influenced by an injection of high-frequency voltage, which is 1 kHz. Figure 11c demonstrates the measured 1 kHz and 25 V high-frequency voltage, which is used to inject into the field winding.

**Figure 11.** Measured signals at 300 r/min and 2 N.m external load; (**a**) a-phase current, (**b**) field-excited winding current, and (**c**) high-frequency voltage.

Figure 12a demonstrates the measured speed responses of step input at 300 r/min. The rise time is 0.23 s, and the speed reaches its steady-state condition at 0.41 s. In addition, the estimated speed ω*rm* is very close to the speed ω*rm* . Figure 12b demonstrates the estimated and real rotor positions and both of them are very close. Figure 12c demonstrates the errors of the estimated position and they are near ±2 electrical degrees.

**Figure 12.** Measured speed and position responses using a multi-step speed predictive controller. (**a**) Speeds, (**b**) positions, and (**c**) errors of estimated position.

Figure 13a shows the different PI parameters of transient responses; Figure 13b shows the different PI parameters of 2 N.m load disturbance responses. From Figure 13a,b, the best PI parameters, considering both transient responses and load disturbance responses, can be obtained as *Kp* = 0.1 and *KI* = 4.5, which is designed by pole assignment with the following two distinct poles, *P*<sup>1</sup> = −1.9 + *j*12.4 and *P*<sup>2</sup> = −1.9 − *j*12.4.

**Figure 13.** Measured speed responses with different proportional-integral (PI) control parameters. (**a**) Transient and (**b**) load disturbance.

Figure 14a demonstrates the measured speed responses using different controllers, including a multi-step predictive speed controller (*HP* = 2), a single-step predictive speed controller (*HP* = 1), and a PI speed controller. As one can observe, the multi-step predictive speed controller has a very smooth response, but the single-step controller has a quicker response with larger overshoot. The PI controller performs the worst, with the slowest response and largest overshoot. Figure 14b demonstrates the measured load disturbance responses. The multi-step controller has a very similar response to the single-step controller. The PI controller, however, has the largest speed drop and the slowest recovery time.

**Figure 14.** Measured speed responses. (**a**) Transient and (**b**) load disturbance.

Figure 15a demonstrates the measured speed responses of a triangular command by using an encoder for comparison. The multi-step controller has a similar response to the single-step controller. Figure 15b demonstrates the measured speed responses of a triangular command by using the proposed sensorless method. The speed ripples are increased when compared with using an encoder. Figure 16a demonstrates the measured speed responses of a trapezoidal command using an encoder. The multi-step controller has similar responses to the single-step controller. Figure 16b demonstrates the measured speed responses of a trapezoidal command using the proposed sensorless method. Figure 16c shows the measured responses of a step-input command using the proposed sensorless method. The multi-step controller with a small weighting factor has a slower but smoother response than the multi-step controller, which has a weighting factor of 1.

**Figure 15.** Measured triangular speed tracking responses. (**a**) Using an encoder and (**b**) using the proposed senseless method.

Figure 17a–c demonstrate the comparison of different speed operations. Figure 17a demonstrates the highest speed by using a field-weakening control, which is 1550 r/min with an 1 N.m external load. The control method effectively reduces the d-axis flux to extend the operational speed range. In addition, Figure 17a also demonstrates the highest operational speed without using a field-weakening control, which is 1080 r/min. Figure 17b demonstrates the mid-speed operational range from 100 to 600 r/min with a 2 N.m external load. All of the different speeds have linear responses. Figure 17c demonstrates the lowest operational speed, which is near 4 r/min with a 0.5 N.m external load. The lowest speed has obvious speed ripples. Figure 18a demonstrates the measured line-to-line voltage, *vab*, which has high-frequency PWM modulation pulses. Figure 18b demonstrates the excited winding voltage, *vf* , which is controlled by an H-bridge circuit. Figure 18c demonstrates the excited winding current, *if* , during the transient time interval of the field current regulated control, which is controlled by an H-bridge circuit. The additional high-frequency flux on the core created core loss, which increases core loss and reduces the efficiency of the motor. Roughly speaking, this increased high frequency core loss is near 2% of the total losses of the motor. Comparing the torque and efficiency for analogously dimension PM and induction motors, for a 0.5 KW motor, the efficiency of a PM is near 90%, the efficiency of an induction motor is near 87%, and the efficiency of a field-excited flux-switching motor is 87%. For a 0.5 kW motor, the torque of a PM is 100% as a reference, the torque of an induction motor is 90%, and the torque of a field-excited flux-switching motor is 93%.

**Figure 16.** Measured speed tracking responses. (**a**) Trapezoidal response using an encoder, (**b**) trapezoidal response using the proposed sensorless method, and (**c**) step input using the proposed sensorless method.

**Figure 17.** Measured highest, middle, and lowest operational speeds. (**a**) Highest speed, (**b**) middle speeds, and (**c**) lowest speed.

**Figure 18.** Measured voltages and currents at 300 r/min. (**a**)-line voltage *vab*, (**b**) field-excited voltage *vf* , and (**c**) field-excited current *if* .

#### **8. Conclusions**

A rotor position estimator and three predictive controllers, which include a multi-step speed predictive controller, a single-step speed predictive controller, and a predictive current controller for a field-excited flux-switching motor drive system have been designed and implemented for this paper. A 1 kHz, high-frequency voltage is injected into the field winding to obtain the estimated rotor position. In addition, the high-frequency injection voltage does not occupy any available output voltages in the inverter. The methods in this paper provide more available voltage for the PWM modulation. The two different predictive controllers implemented in this research enhance the transient responses, decrease the speed drops of load disturbances, and improve tracking responses more effectively than previously published control methods for flux-switching drive systems. In addition, the multiple-step predictive controller provides smoother speed responses than the single-step predictive controller. These are the main contributions of this paper versus other similar research.

Experimental results show that the controllable speed range of the motor drive system discussed in this paper is from 4 r/min to 1500 r/min. Furthermore, the measured results demonstrate that the errors of the rotor position estimations are below ±2 electrical degrees. Moreover, the proposed predictive controllers provide better performance than PI controllers at different speeds. However, all the considerations in the paper are based on a circuit model with fixed parameters. In the future, methods using field methods or online measuring methods will be thoroughly investigated.

Most functions of the rotor position estimators and predictive controllers are implemented by the DSP, and only simple hardware circuits are employed. This proposed flux-switching drive system can be easily applied to meet the high dynamic requirements for applications in household appliances, such as vacuum cleaners and lawn mowers.

**Author Contributions:** Conceptualization, T.-H.L.; Methodology, Y.-H.X.; Software, Y.-H.X.; Validation, T.-H.L., M.S.M., and Y.-H.X.; Formal Analysis, T.-H.L.; Investigation, T.-H.L. and Y.-H.X.; Resources, T.-H.L.; Data Curation, M.S.M.; Writing—Original Draft Preparation, T.-H.L.; Writing—Review & Editing, T.-H.L., and M.S.M.; Visualization, M.S.M.; Supervision, T.-H.L.; Project Administration, T.-H.L.; Funding Acquisition, T.-H.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Ministry of Science and Technology under grant number MOST 108-2221-E-011-085 and -086.

**Acknowledgments:** The paper was supported by the MOST, Taiwan under Grant MOST 108-2221-E-011-085 and -086.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
