*Article* **Online MTPA Trajectory Tracking of IPMSM Based on a Novel Torque Control Strategy**

#### **Jianxia Sun, Cheng Lin \*, Jilei Xing and Xiongwei Jiang**

School of Mechanical Engineering, Beijing Institute of Technology, Beijing 100081, China **\*** Correspondence: lincheng@bit.edu.cn; Tel.: +86-137-0135-3142

Received: 1 August 2019; Accepted: 22 August 2019; Published: 24 August 2019

**Abstract:** The maximum-torque-per-ampere (MTPA) scheme is widely used in the interior permanent magnet synchronous machine (IPMSM) drive system to reduce copper losses. However, MTPA trajectory is complicated to solve analytically. In order to realize online MTPA trajectory tracking, this paper proposes a novel torque control strategy. The torque control is designed to be closed form. Considering the machine reluctance torque as the torque feedback, when this is compared with the torque reference, then the excitation torque reference can be obtained. Since the excitation torque is proportional to the q-axis current, the q-axis current reference can be fed by the excitation torque reference through a proportional regulator. Once the q-axis current reference is given, the d-axis current reference can be calculated based on the per-unit model, which aims to simplify the calculation and make the control strategy independent of machine parameters. In this paper, the stability of the control system is demonstrated. Meanwhile, simulation and experiment results show this torque control strategy can realize MTPA trajectory tracking online and have success in transients.

**Keywords:** maximum-torque-per-ampere (MTPA); torque control; per unit; IPMSM

#### **1. Introduction**

Due to outstanding characteristics, such as high-power density, wide speed range, and a high torque-to-inertia ratio, interior permanent magnet synchronous machines (IPMSMs) are suitable for many industrial applications [1–3]. Differing from the surface permanent magnet synchronous machine (SPMSM), the IPMSM contains not only an excitation torque but also a reluctance torque. Moreover, the excitation torque is proportional to the q-axis current and the reluctance torque is relative to both d- and q-axis stator current. In general, the d-axis current reference is set to zero, and the reluctance torque is consequently null. Then, the torque of IPMSM is proportional to the q-axis current. As a result, this method is easy to implement, but inefficient. In order to improve the efficiency, the maximum-torque-per-ampere (MTPA) scheme, which aims to minimize copper losses of IPMSM during operating, is proposed [4,5]. Under the MTPA condition, the reluctance torque is fully employed and the amplitude of stator current is minimum.

One way for MTPA implementation is finding MTPA operating points directly. For a fixed torque, the optimal current can be solved by a quartic equation, resulting from making the differentiation of the torque equation with respect to stator current equal to zero [6]. In [7], the quartic equation, yielded by the optimal problem, was solved by utilizing Ferrari's method. However, the analytic solution is so complicated that it is hard to utilize online. Therefore, some approximate methods are proposed to avoid complex calculation. For instance, look-up-table (LUT) is a widely used method in IPMSMs. In the LUT method, the optimal current is calculated in advance and sorted in a table. During operating, the optimal current is taken from the premade table. In order to improve accuracy, the LUT method requires accurate machine parameters. In [8], a 3-D lookup table combined with a PI compensator was proposed, to eliminate the undesired deviation which results from the machine

parameter uncertainties. However, the lookup table, especially the 3-D lookup table, is huge and consumes system memory.

In order to avoid complex calculation procedures and get rid of machine parameters, some approaches have been proposed—for instance, the searching algorithm and signal injection method. Machine parameters are not required in these methods, which ensures robustness against parameter variations. In the searching algorithm [9–11], the MTPA operation is achieved by continuously disturbing the current reference, until the minimum of the stator current at the given load torque is found. The advantage of the searching algorithm is to get lower parameter dependence. However, if the convergence rate of the searching algorithm is slow, the algorithm may fail in transients. The signal injection method achieves the MTPA operation by injecting a real or a virtual signal into the reference to obtain the optimum operating points [12–15]. The signal injection method can rapidly respond to transients, whereas the injection of the perturbation signal results in disturbance and harmonics, thus, the control accuracy is reduced.

Instead of finding the optimal current directly, some approaches propose using an outer loop controller for MTPA operation [16–18]. In these methods, the relationship between the d- and q-axis stator currents is used for MTPA operation. The d-axis current reference is generated from the q-axis current, which is provided by an outer loop controller. However, the equation derived from the relationship is complex. Besides this, the outer loop controller is usually a proportional–integral (PI) controller, which results in a high-order mathematical model of the control system, and it is difficult to design a proper PI controller which guarantees global stability.

In this paper, a novel torque control strategy is proposed to track the MTPA trajectory online. The strategy utilizes the per unit model to simplify the calculation and make the model independent of machine parameters. Combined with the MTPA condition expressed in the per unit model, the torque controller can provide the reference stator currents for the current regulators. Besides this, the torque controller is a simple proportional controller. This strategy can respond to transients rapidly. This paper is organized as follows. The per unit modeling procedure and the MTPA condition are described in Section 2. The torque control strategy, system stability analysis, and convergence rate estimation are reported in Section 3. The simulation and experiment results are discussed in Section 4. Finally, the conclusive considerations are drawn in Section 5.

#### **2. System Model and MTPA Condition**

#### *2.1. Per Unit Model of IPMSM*

The torque equation of the IPMSM in a rotor reference frame is given as follows:

$$T\_d = 1.5 P \varphi\_f i\_q + 1.5 P (L\_d - L\_q) i\_d i\_{q\_f} \tag{1}$$

where *Te* is the electromagnetic torque; *P* is the number of pole pairs; ϕ*<sup>f</sup>* is the magnetic flux linkage; *id*, *iq* are the d- and q-axis stator currents; and *Ld*, *Lq* are the d- and q-axis inductances. The electromagnetic torque *Te* consists of reluctance torque *T*<sup>1</sup> and excitation torque *T*2. Moreover, the equation of *T*<sup>1</sup> and *T*<sup>2</sup> can be expressed as:

$$\begin{cases} \quad T\_1 = 1.5P(L\_d - L\_q)i\_d i\_q\\ \quad T\_2 = 1.5Pq\_f i\_q \end{cases} \tag{2}$$

Equations (1) and (2) reveal that different machine parameters lead to varied torques. Therefore, the universality of the model is poor. In order to solve this problem, a per unit model of IPMSM is presented in this paper.

To build the per unit model, the base values must be defined. Moreover, the relationship between the base current and the base torque should correspond to (1). In this paper, the base current *Ib* is defined as <sup>ϕ</sup>*<sup>f</sup>* <sup>2</sup>(*Lq*−*Ld*) . Consequently, the base torque *Tb* is equal to:

$$T\_b = 1.5 P \varphi\_f I\_b + 1.5 P \mathbf{(L}\_d - L\_q) I\_b^2 = 0.75 P \varphi\_f I\_b. \tag{3}$$

According to (1), the torque of IPMSM in per unit model can be expressed as:

$$\begin{cases} T\_{cn} = 2i\_{qu} - i\_{qu}i\_{dn} \\ T\_{1n} = -i\_{qu}i\_{dn} \\ T\_{2n} = 2i\_{qu} \end{cases} \tag{4}$$

where, *Ten* = *Te Tb* , *<sup>T</sup>*1*<sup>n</sup>* = *<sup>T</sup>*<sup>1</sup> *Tb* , *<sup>T</sup>*2*<sup>n</sup>* = *<sup>T</sup>*<sup>2</sup> *Tb* , *idn* <sup>=</sup> *id Ib* , and *iqn* <sup>=</sup> *iq Ib* . *Ten* is the per unit value of *Te*, *T*1*<sup>n</sup>* is the per unit value of *T*1, *T*2*<sup>n</sup>* is the per unit value of *T*2, *idn* is the per unit value of *id*, and *iqn* is the per unit value of *iq*.

The per unit model can be applied to all IPMSMs, because it is independent of machine parameters. The machine parameters can be measured offline or estimated online [19].

#### *2.2. MTPA Condition*

The magnitude of stator current *Is* is equal to *i* 2 *<sup>d</sup>* + *<sup>i</sup>* 2 *<sup>q</sup>*, which can be represented in per unit values as:

$$
\dot{l}\_{sn}^2 = \dot{l}\_{dn}^2 + \dot{l}\_{qn\prime}^2 \tag{5}
$$

where *Isn* = *Is Ib* is the per unit value of *Is*. According to (4), the relationship between *iqn* and *idn* can be expressed as:

$$i\_{qn} = \frac{T\_{cn}}{2 - i\_{dn}}.\tag{6}$$

Equation (6) reveals that the relationship between *idn* and *iqn* is fixed for a given torque *Ten*. Inserting (6) into (5), the relationship between *Isn* and *idn* can be obtained as:

$$
\dot{l}\_{sn}^2 = \dot{i}\_{dn}^2 + \frac{T\_{cn}^2}{\left(2 - i\_{dn}\right)^2}.\tag{7}
$$

The minimum value of *Isn* can be obtained by differentiating (7) with respect to *idn*, and setting the resulting expression to zero yields,

$$\frac{dI\_{sn}^2}{di\_{dn}} = 2i\_{dn} + \frac{2T\_{cn}^2}{\left(2 - i\_{dn}\right)^3} = 0.\tag{8}$$

Equation (8) is a quartic equation and hard to solve, i.e., it is difficult to get the optimal solution of *idn* and *iqn* for a given torque *Ten* in MTPA condition. In order to solve this problem, a novel torque control strategy is proposed in this paper, which will be described in the following sections.

#### **3. Torque Control and Stability Analysis**

#### *3.1. Torque Control with MTPA Trajectory Tracking*

For the sake of MTPA trajectory tracking, a novel torque control strategy is introduced (Figure 1). The reference torque *Te*\_*re f* is given by an outer controller, and only positive reference torque will be taken into account in this paper. Considering the machine reluctance torque *T*<sup>1</sup> as the torque feedback, comparing with the torque reference *Te*\_*re f* , then the excitation torque reference *T*<sup>2</sup> can be

obtained. Since the excitation torque is proportional to the q-axis current, the q-axis current reference *i* ∗ *<sup>q</sup>* can be fed by the excitation torque reference *T*<sup>2</sup> through a proportional regulator.

**Figure 1.** The torque control strategy structure with maximum-torque-per-ampere (MTPA) conditions.

For the sake of simplicity, the per unit model is used in this strategy. Combing (2) and (4), *iqn* is equal to *T*2/2*Tb*. Once *iqn* is given, the per unit value *idn* in MTPA condition can be solved as follows.

By combining (6), Equation (8) can be simplified to:

$$(i\_{qn}^2 + (2 - i\_{dn})i\_{du} = 0.\tag{9}$$

Solving *idn* gives:

$$i\_{du} = 1 - \sqrt{1 + i\_{qu}^2}.\tag{10}$$

Since *idn* and *iqn* have been solved, the currents reference vector (*i* ∗ *d*, *i* ∗ *<sup>q</sup>*) can be derived easily. Within the effects of torque control, (*i* ∗ *d* , *i* ∗ *<sup>q</sup>*) will converge asymptotically to the MTPA operating point (*i* ∗ *<sup>d</sup>*,*mtpa*, *i* ∗ *<sup>q</sup>*,*mtpa*).

The dynamic analysis findings are as follows: if the *q*-axis current reference *i* ∗ *<sup>q</sup>* > *i* ∗ *<sup>q</sup>*,*mtpa*, then the *d*-axis current reference *i* ∗ *<sup>d</sup>* < *i* ∗ *<sup>d</sup>*,*mtpa*, the machine reluctance torque *T*<sup>1</sup> increases, then excitation torque reference *T*<sup>2</sup> will decrease, and *i* ∗ *<sup>q</sup>* will decrease consequently; in the case of *i* ∗ *<sup>q</sup>* < *i* ∗ *<sup>q</sup>*,*mtpa*, it is exactly the same way. This torque controller keeps adjusting *i* ∗ *<sup>q</sup>* until it is equal to *i* ∗ *<sup>q</sup>*,*mtpa* and then *i* ∗ *<sup>d</sup>* is equal to *i* ∗ *<sup>d</sup>*,*mtpa*, i.e., (*i* ∗ *<sup>d</sup>*,*mtpa*, *i* ∗ *<sup>q</sup>*,*mtpa*) is the equilibrium point of (*i* ∗ *d*, *i* ∗ *<sup>q</sup>*). The relationship between the given torque and MTPA operating point is shown in Figure 2.

**Figure 2.** The relationship between the given torque and MTPA operating point.

The reference currents can be used properly in the field-oriented control (FOC) scheme to achieve MTPA trajectory tracking.

#### *3.2. Stability Analysis and Convergence Rate Estimation*

The dynamic analysis describes roughly the process of tracking MTPA operating points. In this section, the stable performance of torque control strategy will be proven as follows.

Within the proper design of vector control, the d- and q-axis current closed loop systems can be approximated by a low pass filter [20], i.e.,:

$$\begin{cases} \dot{i}\_d = \frac{1}{\pi} (\dot{i}\_d^\* - \dot{i}\_d) \\ \dot{i}\_q = \frac{1}{\pi} (\dot{i}\_q^\* - \dot{i}\_q) \end{cases} \tag{11}$$

where, τ is the time constant of current closed loop. Denoting the torque *T*∗ *<sup>e</sup>*, *T*<sup>∗</sup> <sup>1</sup> and *T*<sup>∗</sup> <sup>2</sup> are given by *i* ∗ *<sup>d</sup>* and *i* ∗ *<sup>q</sup>*. Meanwhile, their relationships should correspond to (1) and (2), i.e., *T*<sup>∗</sup> <sup>1</sup> = 1.5*P* - *Ld* − *Lq i* ∗ *d i* ∗ *q*, *T*∗ <sup>2</sup> = 1.5*P*ϕ*<sup>f</sup> i* ∗ *<sup>q</sup>*, and *T*<sup>∗</sup> *<sup>e</sup>* = *T*<sup>∗</sup> <sup>1</sup> + *T*<sup>∗</sup> 2.

Making the time derivative of *T*<sup>1</sup> and combining the result with (11), the following equation can be obtained:

$$
\tau \dot{T}\_1 = -2T\_1 + 1.5P(L\_d - L\_q)i\_d i\_q^\* + 1.5P(L\_d - L\_q)i\_q i\_d^\*.\tag{12}
$$

Differentiating (12) with respect to time again and then inserting the result into (12), it follows that:

$$
\tau^2 \ddot{T}\_1 + 3\tau \dot{T}\_1 + 2T\_1 = 2T\_1^\*. \tag{13}
$$

As expressed in Figure 1, *T*∗ <sup>2</sup> equals to *<sup>T</sup>*2, and *<sup>T</sup>*<sup>2</sup> equals to *Te*\_*re f* <sup>−</sup> *<sup>T</sup>*1. Thus, . *T* ∗ <sup>2</sup> <sup>=</sup> <sup>−</sup> . *T*1, .. *T* ∗ <sup>2</sup> = − .. *T*1, and *T*<sup>1</sup> = *Te*\_*re f* − *T*<sup>∗</sup> <sup>2</sup>. Inserting the results into (13), it holds that:

$$
\tau^2 \ddot{T}\_2^\* + 3\tau \dot{T}\_2^\* = 2(T\_{e\\_ref} - T\_e^\*).\tag{14}
$$

According to the singular perturbation approximation principle [9], the second order term <sup>τ</sup><sup>2</sup> .. *T* ∗ 2 can be neglected because τ<sup>2</sup> is far less than 3τ. Combing with *T*<sup>∗</sup> <sup>2</sup> = 1.5*P*ϕ*<sup>f</sup> i* ∗ *<sup>q</sup>*, the derivative of *i* ∗ *<sup>q</sup>* is given by:

$$\dot{i}\_q^\* = \frac{4(T\_{\varepsilon\,ref} - T\_{\varepsilon}^\*)}{9\pi P q\_f}.\tag{15}$$

Based on (10), the derivative of *i* ∗ *<sup>d</sup>* can be expressed as:

$$\dot{i}\_d^\* = -\frac{8\dot{i}\_q^\*(T\_{c\\_ref} - T\_c^\*)}{9\pi P \wp\_f \sqrt{1 + \left(\dot{i}\_q^\* / I\_b\right)^2}}.\tag{16}$$

Apparently, *Te*\_*re f* = *T*<sup>∗</sup> *<sup>e</sup>* when . *i* ∗ *<sup>q</sup>* <sup>=</sup> 0 and . *i* ∗ *<sup>d</sup>* = 0, i.e., *Te*\_*re f* = *T*<sup>∗</sup> *<sup>e</sup>* when (*i* ∗ *d*, *i* ∗ *<sup>q</sup>*) converge to the MTPA operating point. Considering the Lyapunov function *V* = <sup>1</sup> 2 - *Te*\_*re f* − *T*<sup>∗</sup> *e* 2 , the derivative of *V* is given by (17): .

$$\begin{split} \dot{V} &= \dot{T}\_{\varepsilon}^{\*} \Big( T\_{\varepsilon}^{\*} - T\_{\varepsilon \text{-} ref} \Big) = \begin{pmatrix} \frac{\partial T^{\*}}{\partial \dot{t}\_{d}^{\*}} \frac{d \dot{t}\_{d}^{\*}}{dt} + \frac{\partial T^{\*}}{\partial \dot{t}\_{q}^{\*}} \frac{d \dot{t}\_{q}^{\*}}{dt} \Big) \Big( T\_{\varepsilon}^{\*} - T\_{\varepsilon \text{-} ref} \Big) \\ &= -\left[ \frac{3}{2} P \big( L\_{d} - L\_{q} \big) \dot{t}\_{d}^{\*} + \frac{3}{2} P q\_{f} - \frac{3 P \big( L\_{d} - L\_{q} \big) \dot{t}\_{q}^{\*2}}{\sqrt{1 + \left( \dot{t}\_{q}^{\*} / I\_{b} \right)^{2}}} \right] \frac{4 \left( T\_{\varepsilon, ref} - T\_{\varepsilon}^{\*} \right)^{2}}{9 \pi P q\_{f}} \end{split} \tag{17}$$

Since only positive reference torque *Te*\_*re f* is taken into account, *i* ∗ *<sup>q</sup>* > 0 and *i* ∗ *<sup>d</sup>* < 0. Besides, *Ld* < *Lq*. Therefore, *<sup>V</sup>* <sup>≥</sup> 0 and . *V* ≤ 0, i.e., this system is asymptotically stable.

According to (17), . *<sup>V</sup>* ≤ −4*<sup>V</sup>* <sup>3</sup><sup>τ</sup> , solving *V* gives:

$$
\Delta V \le V\_0 e^{-\frac{4(t-t\_0)}{3\tau}},
\tag{18}
$$

where, *t*<sup>0</sup> is the initial time, *V*<sup>0</sup> is the initial value of *V*. This shows that the convergence rate of *V* is higher than <sup>4</sup> 3τ .

According to (11), solving (*id*, *iq*) gives:

$$\begin{cases} \ i\_d = i\_d^\* - \left( i\_d^\* - i\_{d0} \right) e^{-\frac{1}{\tau} \left( t - t\_0 \right)} \\\ i\_q = i\_q^\* - \left( i\_q^\* - i\_{q0} \right) e^{-\frac{1}{\tau} \left( t - t\_0 \right)} \end{cases} \tag{19}$$

where (*id*0, *iq*0) are the initial values of (*id*, *iq*). Equation (19) shows the convergence rate of (*id*, *iq*) equals <sup>1</sup> <sup>τ</sup> . Obviously, the convergence rate of *V* is higher than (*id*, *iq*). That means the convergence rate of the torque loop is higher than the current loop, i.e., *T*∗ *<sup>e</sup>* coverages to *Te*\_*re f* and (*i* ∗ *d* , *i* ∗ *<sup>q</sup>*) converge to (*i* ∗ *<sup>d</sup>*,*mtpa*, *i* ∗ *<sup>q</sup>*,*mtpa*) before (*id*, *iq*) converge to (*i* ∗ *d* , *i* ∗ *<sup>q</sup>*). Therefore, this torque control strategy is suitable for dynamic MTPA trajectory tracking.

#### **4. Results**

#### *4.1. Simulation Results*

To validate the proposed torque control strategy, the whole system was simulated in Matlab/Simulink. The parameters of interior permanent magnet synchronous machine (IPMSM) are tabulated in Table 1.


**Table 1.** Parameters of IPMSM.

The rated speed was set to 1000 rpm. The machine speed was controlled by an outer speed controller, which provided the reference torque *Te*\_*re f* . In order to maintain a constant speed, *Te*\_*re f* was varied with the load torque *TL*. Based on the different state of *TL*, the simulation was divided into three stages. Stage 1: 0 s~0.2 s, *TL* increased from 0 Nm to 50 Nm equably. Stage 2: 0.2 s~0.4 s, *TL* maintained at 50 Nm. Stage 3: 0.4 s~1 s, *TL* dropped to 40 Nm at 0.4 s, then remained constant until 1 s. Figure 3 shows the torque and current response of MTPA trajectory tracking under the torque control strategy. Figure 3a shows the drop to 40 Nm at 0.4 s, then remaining constant until 1 s. Figure 3 shows the torque and current response of MTPA trajectory tracking under the torque control strategy. Figure 3a shows the response of the machine speed. Figure 3b shows the response of reference torque *Te*\_*re f* and machine torque *Te*. Figure 3c shows the response of the d-axis reference current *i* ∗ *<sup>d</sup>* and stator current *id*. Figure 3d shows the response of the q-axis reference current *i* ∗ *<sup>q</sup>* and stator current *iq*. It can be seen that *Te* followed *Te*\_*re f* closely in these load conditions, even if *TL* changed continuously. Besides, *id* and *iq* varied with *Te*\_*re f* .

**Figure 3.** The simulation results: (**a**) machine speed; (**b**) reference and machine torque; (**c**) d-axis reference current and stator current; (**d**) q-axis reference current and stator current.

Figure 4 reports the transient current locus and the desired MTPA trajectory. It can be seen that (*id*, *iq*) converged to the MTPA operating point asymptotically. The simulation results reveal that this online MTPA trajectory tracking strategy is excellent in steady and transient conditions.

**Figure 4.** The transient current locus and the desired MTPA trajectory.

#### *4.2. Experimental Results*

The proposed control approach was validated through experimental investigations, which were conducted on a laboratory IPMSM drive system, shown in Figure 5.

The experiment platform contained a 3 kW IPMSM, which was used for applying the proposed strategy and an asynchronous motor operated as a programmable load, which was controlled by the PC. The tested motor was driven by the three-phase voltage source PWM (Pulse Width Modulation) inverter and the controller was implemented in the TMS320F28335 DSP (Digital Signal Processor). The position of the tested motor was estimated by a sliding-mode observer. The machine reluctance

torque was calculated based on the machine parameters. The parameters of the IPMSM were tabulated in Table 1.

**Figure 5.** Experimental setup.

The rated speed was set to 1000 rpm. The machine speed response is shown in Figure 6a. The response of the reference torque *Te*\_*re f* is shown in Figure 6b. The response of the d-axis reference current *i* ∗ *<sup>d</sup>* and stator current *id* are shown in Figure 6c. The response of the q-axis reference current *i* ∗ *q* and stator current *iq* are shown in Figure 6d. Figure 7 shows the transient response of the reference torque and q-axis reference current. It can be seen that the stator current can follow the current reference closely. The torque reference was adjusted with speed and the current reference was adjusted with torque reference. The experiment results reveal that this strategy is active.

**Figure 6.** The experiment result: (**a**) machine speed; (**b**) reference torque; (**c**) d-axis reference current and d-axis stator current; (**d**) q-axis reference current and q-axis stator current.

**Figure 7.** The experiment result: (**a**) reference torque; (**b**) q-axis reference current.

#### **5. Conclusions**

A novel torque controller was presented in this paper. Combined with the MTPA condition, this torque control strategy can track the MTPA trajectory online. The reference torque consisted of two parts, excitation torque and reluctance torque. The excitation torque was fed to the proportional regulator, which produced the q-axis current reference. Once the q-axis current reference was given, the d-axis current reference could be calculated online based on the MTPA condition in the per unit model. Moreover, the implementation of the per unit model can simplify the MTPA condition, reduce the calculation burden, and make the machine parameters independent. In addition, parameter estimation algorithms can be conveniently combined with the torque controller to improve accuracy. This strategy combined the MTPA condition with the torque controller to realize MTPA trajectory tracking online. Moreover, the stability and dynamic performance were proven analytically. The simulation and experiment results show the excellent performance of the strategy in steady and transient conditions.

**Author Contributions:** Data curation, X.J.; Formal analysis, J.S.; Investigation, J.S.; Methodology, J.S.; Project administration, C.L.; Software, J.X.; Supervision, C.L.

**Funding:** This work was supported by Ministry of Science and Technology of the People's Republic of China under Grant 2017YFB0103801.

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

#### **References**


© 2019 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/).
