**Sensorless PMSM Drive Implementation by Introduction of Maximum E**ffi**ciency Characteristics in Reference Current Generation**

#### **Željko Planti´c 1, Tine Marˇciˇc 2, Miloš Bekovi´c <sup>3</sup> and Gorazd Štumberger 3,\***


Received: 31 July 2019; Accepted: 5 September 2019; Published: 11 September 2019

**Abstract:** This paper presents the efficiency improvement in a speed closed-loop controlled permanent magnet synchronous machine (PMSM) sensorless drive. The drive efficiency can be improved by minimizing the inverter and the PMSM losses. These can be influenced by proper selection of DC-bus voltage and switching frequency of the inverter. The direct (d-) and quadrature (q-) axis current references generation methods, discussed in this paper, further improve the efficiency of the drive. Besides zero d-axis current reference control, the maximum torque per ampere (MTPA) characteristic is normally applied to generate the d- and q-axis current references in vector controlled PMSM drives. It assures control with maximum torque per unit of current but cannot assure maximum efficiency. In order to improve efficiency of the PMSM drive, this paper proposes the generation of d- and q-axis current references based on maximum efficiency (ME) characteristic. In the case study, the MTPA and ME characteristics are theoretically evaluated and determined experimentally by measurements on discussed PMSM drive. The obtained characteristics are applied for the d- and q-axis current references generation in the speed closed-loop vector controlled PMSM drive. The measured drive efficiency clearly shows that the use of ME characteristic instead of MTPA characteristic or zero d-axis current in the current references generation improves the efficiency of PMSM drive realizations with position sensor and without it—sensorless control.

**Keywords:** magnetic loss; maximum efficiency (ME) characteristic; maximum torque per ampere (MTPA) characteristic; modeling; permanent magnet synchronous machine (PMSM); sensorless control

#### **1. Introduction**

The energy savings during the entire exploitation time of an electric drive can be achieved only through a coordinated efficiency improvement of all components of the drive, including control algorithms. Nowadays, efficiency maximization is a normal design aspect for electric drives. Thus, more and more researchers are occupied with decreasing overall power losses of electric drives and simultaneously saving component and operational costs. In some applications, minimum loss is not only one of the design criteria, but is of paramount importance, e.g., in battery-powered road-vehicles [1,2] or even more importantly in autonomous aircrafts [3]. In such applications, even a small decrease of losses can have a significant impact on overall system performance in terms of ensuring adequate operating range. Due to its high efficiency, high power density and wide operating range of speed, the permanent magnet synchronous machine (PMSM) is used very often. By employment of a speedor position-sensorless control, a more compact drive with less maintenance [4–6] can be obtained.

Besides designing a new machine, the efficiency of the drive can be improved by minimizing losses of the power inverter by proper selection of the DC-bus voltage and the switching frequency [7–9], and a proper control of the PMSM as well. Sensorless PMSM drive control algorithms vary in complexity and computational performance demanded by different target applications. The main difference between available sensorless control algorithms from the literature is that model-based methods [5,10–13] are more straightforward to implement but are not directly applicable at standstill, and therefore saliency-based methods [14–16] have been developed. Authors combine the aforementioned methods in various approaches [9,17–20], thus ensuring adequate sensorless performance from zero speed up to the field-weakening region.

In the case of vector controlled PMSM, a proper selection of the direct (d-) and quadrature (q-) axis current references can further improve drive efficiency. The simplest zero d-axis current reference control (*idr* = 0) is still employed in applications with surface mounted PMSMs [21], which exhibit similar values of d- and q-axis inductances. Although it is generally applicable, the maximum torque per ampere (MTPA) characteristic [10,15,19–27], which considers PMSM copper losses, is usually used for generation of current references in drives with interior PMSM, which exhibit different values of dand q-axis inductances and thus produce a significant reluctance torque component in addition to the permanent magnet torque component. When MTPA is applied, the PMSM drive produces maximum torque per unit of stator current but cannot reach the maximum efficiency due to the neglected iron core losses. The maximum efficiency (ME) characteristic [21,23–28] considers copper losses as well as iron core losses of PMSM. When it is applied for generation of current references, the maximum efficiency of a PMSM drive can be achieved. The concept was first evaluated by Colby et al. [29] for a surface mounted PMSM. Later it was applied by Morimoto et al. [30] in control of an interior PMSM drive as well. More recently, Pohlenz et al. [23] and Cavallaro et al. [28] reported improved efficiencies by implementing ME in the control software. While applying ME, Hassan et al. [24] analyzed PMSM and inverter losses due to changed harmonic content. Ni et al. [21] developed an iron loss formula by introducing new coefficients. However, the needed detailed PMSM data might not always be available, as discussed later.

General forms of the MTPA characteristic and the ME characteristic are shown in Figure 1a, from where it is evident that ME characteristics require higher values of the negative d-axis current in comparison to the MTPA characteristics, and MTPA characteristics in turn require higher values of the negative d-axis current than the *idr* = 0 control. Figure 1b shows the stator current vector **I***<sup>s</sup>* and its components *Id* and *Iq*, where |**I***s*| denotes the length of **I***s*.

**Figure 1.** Maximum torque per ampere (MTPA) and maximum efficiency (ME) characteristics (**a**) and stator current vector in α-β and d-q reference frame (**b**).

In this paper, the origins of losses in a PMSM drive are introduced and evaluated in Section 2. The PMSM model is used to explain the differences between the MTPA and ME characteristics, and to obtain the equations which relate the ME characteristics to the stator currents *Id* and *Iq* in Section 3. Although the ME characteristics are usually derived in relation to the magnetizing currents *Imd* and *Imq*,

the here presented Equation (20) provides the optimal current *Id* for a given load torque by considering the ME characteristic-i.e., the copper and the iron losses of the PMSM. Furthermore, this paper focuses on the drive implementation from the control designer's point of view, where detailed information of the PMSM is not available, contrarily to machine designer's point of view, where all details of the PMSM are known [21]. Unfortunately, the aforementioned equations cannot be used directly in the control implementation because the needed iron core resistance parameter value *Rc* cannot be determined directly or to be modeled in form of a function. Therefore, this paper proposes an alternative, which is suitable for control design and applicable to all types of PMSM. There the ME characteristics are determined experimentally by measurements performed on the tested vector controlled PMSM drive. The existing control algorithm is extended by a current references generator, where the determined ME characteristics are applied to generate the d- and q-axis current references. The speed closed-loop vector control of PMSM, obtained in this way, is used to measure the efficiency of the entire drive. The control is realized with the position sensor and without it—sensorless control.

Moreover, in this work the ME characteristics are implemented in sensorless control for the first time and presented in Section 4. There the rotor speed and position are determined by nonlinear phase-locked loop (PLL) observer [9,31], whose inputs are calculated by the modified active flux method [32–35]. It is shown that the drive efficiency can be improved when ME characteristics are used to generate the d- and q-axis current references, which is true for the sensorless control as well as for the control realization with position sensor. In order to facilitate a direct comparison, the space vector modulation has been utilized without any advanced modulation techniques such as discussed in [36,37]. Section 5 then concludes the paper.

#### **2. PMSM Drive Losses**

A PMSM drive normally consists of the PMSM, power inverter, load, and the control system. The load is the element, which dictates the proper selection of the other PMSM drive elements. The energy conversion process inevitably includes losses, mostly in the power inverter and PMSM. These loss components substantially influence the drive efficiency. This paper focuses on the detailed research of PMSM losses, which are minimized in the speed closed-loop controlled PMSM drive by a proper generation of d- and q-axis current references.

The employed PMSM had surface mounted magnets and was rated for 400 V and 1.6 kW at 2250 rpm.

#### *2.1. Power Losses in Power Inverter*

The power losses in a power inverter consist of conduction or Ohm's losses and switching losses. In this work they are minimized by proper selection of the DC-bus voltage, where the impact of dead times is minimized [5,8,9] and by proper selection of the switching frequency [7,24,36,37] as well. The proper values were determined experimentally by measurements performed on the tested PMSM drive.

In order to evaluate the impact of switching frequency on drive efficiency, a series of tests was performed, where the PMSM load was varied, while the switching frequency was set to 4, 6, and 9 kHz, respectively. Figure 2 presents the PMSM drive efficiency characteristics in dependence of the current vector length |**i**s|. Drive efficiencies vary because the inverter losses increase with the increase of switching frequency and simultaneously the PMSM losses decrease with the increase of switching frequency [38], and the drive efficiency represents their product. Figure 2 shows that the highest efficiency of this particular PMSM drive was obtained when the switching frequency of 6 kHz was employed. Thus, the aforementioned switching frequency value was used through all successive tests.

**Figure 2.** Measured efficiency characteristics of the permanent magnet synchronous machine (PMSM) drive at different switching frequencies.

In addition to the switching frequency, the inverter's DC-bus voltage influences the PMSM drive efficiency as well. When the desired inverter output voltage is much smaller than the DC-bus voltage, then the inverter output voltage pulses are unnecessarily narrow which increases switching losses. However, it is noted that the DC-bus voltage control is not possible if a diode rectifier is employed in the inverter. In order to evaluate the impact of DC-bus voltage on PMSM drive efficiency, a series of tests was performed, where the PMSM load was varied, while the rotor speed and DC-bus voltage were set to different values. The PMSM efficiency was not significantly affected by the change of DC-bus voltage. Thus, the significant variation of PMSM drive efficiency presented in Figure 3 was due to the inverter efficiency variation by the change of DC-bus voltage. The differences in efficiencies at low speeds were more expressed because the PMSM back-emf was smaller at low speeds. The implementation of variable DC-bus voltage could make sense in variable speed drives operating mostly at partial loads and in a broad range of speeds in applications where efficiency is crucial.

**Figure 3.** Measured efficiency characteristics of the PMSM drive at different rotor speed and different DC-bus voltage.

#### *2.2. Power Losses in PMSM*

The analysis of power losses in the PMSM is based on the two-axis steady-state PMSM model with lumped parameters, written in the d-q reference frame. The d-axis was aligned with the flux linkage vector due to the permanent magnets while the q-axis led for electrical π/2. The model is schematically shown in Figure 4.

**Figure 4.** Applied steady-state model of PMSM in d-q reference frame.

The copper losses are considered by the stator resistance *Rs* while the iron core losses are accounted for by the iron core resistance *Rc*. In the steady-state model, written in the d-q reference frame, the indices d and q refer to the d- and q-axis, *Ud* and *Uq* are the stator voltages, *Id* and *Iq* are stator currents, *Imd* and *Imq* are the magnetizing currents, *Icd* and *Icq* are the iron core currents, *Ld* and *Lq* are the stator inductances, Ψ*md* is the flux linkage due to the permanent magnets, *Te* is the electromagnetic torque, ω*<sup>e</sup>* is the electrical rotor speed and *pp* is the number of pole pairs. Considering Figure 4, the voltage Equation (1) and the torque Equation (2) can be written:

$$\begin{cases} \mathcal{U}\_d = \mathcal{R}\_s I\_d - \omega\_\varepsilon L\_q I\_{mq} \\ \mathcal{U}\_q = \mathcal{R}\_\varepsilon I\_q + \omega\_\varepsilon L\_d I\_{md} + \omega\_\varepsilon \mathcal{W}\_{md} \end{cases} \tag{1}$$

$$T\_{\varepsilon} = p\_p I\_{mq} \left( \Psi\_{md} + \left( L\_d - L\_q \right) I\_{md} \right) \tag{2}$$

where

$$\begin{cases} I\_d = I\_{md} + I\_{cd} = I\_{md} - \frac{\omega\_c I\_q I\_{mq}}{R\_c} \\ I\_q = I\_{mq} + I\_{cq} = I\_{mq} + \frac{\omega\_c}{R\_c} (L\_d I\_{md} + \Psi\_{md}) \end{cases} \tag{3}$$

The power losses in the *PMSM Ploss,PMSM* can be divided into controllable losses, which are the copper losses *Ploss,Cu* and the iron core losses *Ploss,Fe*, and uncontrollable losses, which are the mechanical and additional losses [28]. Neglecting the uncontrollable losses, (4) describes the PMSM losses.

$$P\_{loss, PMSM} = P\_{loss, Cu} + P\_{loss, Fe} \tag{4}$$

Considering Figure 4 and Equations (1) to (4), *Ploss,Cu* and *Ploss,Fe* can be expressed by (5) and (6).

$$P\_{\rm loss,Cu} = R\_5 \left( l\_d^2 + l\_q^2 \right) = R\_5 \left( I\_{md} - \frac{\alpha \nu L\_q}{R\_c} I\_{mq} \right)^2 + R\_5 \left( I\_{mq} + \frac{\alpha \nu \left( L\_d I\_{md} + \Psi\_{md} \right)}{R\_c} \right)^2 \tag{5}$$

$$P\_{loss,Fe} = R\_c \left( I\_{cd}^2 + I\_{cq}^2 \right) = \frac{\alpha\_c^2}{R\_c} \left[ \left( L\_q I\_{mq} \right)^2 + \left( L\_d I\_{md} + \Psi\_{md} \right)^2 \right] \tag{6}$$

#### **3. Improving E**ffi**ciency by Reference Current Generation**

Assuming a given inverter and a given PMSM, the PMSM drive efficiency can be improved by minimizing inverter losses and PMSM losses. In addition to advanced modulation techniques [36,37], the inverter losses can be mostly influenced by the closed-loop control or a proper choice of DC-bus voltage and by an optimal choice of switching frequency [24]. The PMSM losses can be influenced by the voltage and current total harmonic distortion that depend on the switching frequency [7,24,38], and in the case of vector controlled PMSM by the generation of the d- and q-axis current references. The generation of current references is usually based *idr* = 0 or on the MTPA characteristic [10,15,19–23], which assures the maximum torque per unit of current but not the maximum efficiency of the entire drive [21,23]. The last can be achieved when the ME characteristic is applied for generation of the dand q-axis current references.

#### *3.1. Theoretical Background of the MTPA and ME Characteristics*

Let us consider all parameters of the PMSM steady-state model, shown in Figure 4 and described by Equations (1) to (6), as constant. The starting point for derivation of the MTPA characteristic is the torque equation (Equation (2)) and the PMSM steady-state model shown in Figure 4, where the iron core resistance is neglected (*Rc* →∞). Consequently, according to Equation (3), *Icd* and *Icq* equal zero which leads to *Id* = *Imd* and *Iq* = *Imq*. Considering the angle of the stator current vector α, shown in Figure 1b, the relation between the d- and q-axis currents, that provides the maximum torque per unit of current, can be determined from Equation (2) by Equation (7), which yields Equation (8) [22].

$$\max(T\_{\varepsilon}) \Rightarrow \frac{dT\_{\varepsilon}}{d\alpha} = 0 \tag{7}$$

$$I\_{d\_{1,2}} = \frac{-\Psi\_{md}}{2\left(L\_d - L\_q\right)} \pm \sqrt{\frac{\Psi\_{md}^2}{4\left(L\_d - L\_q\right)^2} + I\_q^2} \tag{8}$$

Similarly, the ME characteristic describes the relation between the d- and q-axis current that provides the minimum of PMSM losses, which leads to the maximum efficiency of the PMSM drive. According to [24,28], the PMSM losses (Equation (4)) can be minimized for each combination of *Te*, and ω*e*. Considering Equation (2) in Equations (5) and (6), Equation (4) can be expressed as a function of *Imd*, *Te*, and ω*e*. In the case of steady-state operation, *Te*, and ω*<sup>e</sup>* are constants, which means that the PMSM power losses *Ploss,PMSM* (Equation (4)) can be described as a function of *Imd*. However, because *Te* (Equation (2)) depends on *Imd* and *Imq*, the PMSM power losses can be minimized by Equation (9), which yields Equation (10) [24] describing the ME characteristic.

$$\min\left(P\_{\text{loss},PMSM}\right) \Rightarrow \left.\frac{dP\_{\text{loss},PMSM}\left(I\_{\text{md},\epsilon}I\_{\text{mq},\epsilon}\omega\_{\epsilon}\right)}{dI\_{\text{md}}}\right|\_{T\_{\epsilon},\omega\_{\epsilon}} = 0\tag{9}$$

$$I\_{md} = -\frac{\omega\_\epsilon^2 L\_d \Psi\_{md} (R\_\epsilon + R\_\epsilon) + R\_\epsilon R\_\epsilon \omega\_\epsilon (L\_d - L\_q) I\_{mq}}{\omega\_\epsilon^2 L\_d^2 (R\_\epsilon + R\_\epsilon) + R\_\epsilon R\_\epsilon^2} \tag{10}$$

#### *3.2. Evaluation of the Impact of Rc on the ME Characteristic*

A sensitivity analysis was performed in order to evaluate the impact of the iron core resistance *Rc* on the ME characteristic. The parameter *Rc* was varied while the other PMSM parameters presented in Table 1 remained constant. The theoretical derivation of obtaining the necessary measurable current *Id* for every corresponding combination of *Iq* and ω*<sup>e</sup>* is presented hereinafter with the assumption of constant PMSM model parameters in the given operating points.

**Table 1.** Permanent magnet synchronous machine PMSM parameters.


The voltage equations derived from Figure 4 can be rewritten also in the form of Equation (11).

$$\begin{cases} \mathcal{U}\_d = R\_s I\_{md} - \frac{R\_s + R\_c}{R\_c} \omega\_c L\_q I\_{mq} \\ \mathcal{U}\_q = R\_s I\_{mq} + \frac{R\_s + R\_c}{R\_c} (\omega\_c L\_d I\_{md} + \omega\_c \Psi\_{md}) \end{cases} \tag{11}$$

Considering Equation (3), the currents *Imd* and *Imq* can be expressed in terms of the measurable stator currents *Id* and *Iq* (Equation (12)). Note that it is essential for control implementation that measurable quantities (such as *Id* and *Iq* instead of *Imd* and *Imq*) are included in the algorithm.

$$\begin{array}{l} I\_{md}(I\_d, I\_q) = \frac{R\_c^2 I\_d + R\_c L\_q \omega\_d I\_q - I\_q \mathcal{W}\_{md} \omega\_c^2}{R\_c^2 + L\_d I\_q \omega\_c^2} \\\ I\_{mq}(I\_d, I\_q) = \frac{R\_c^2 I\_q - R\_c L\_d \omega\_d I\_d - R\_c \mathcal{W}\_{md} \omega\_c}{R\_c^2 + L\_d I\_q \omega\_c^2} \end{array} \tag{12}$$

Similarly, the currents *Icd* and *Icq* can be expressed in terms of *Id* and *Iq* (13).

$$\begin{split} I\_{\rm cd}(I\_{d}, I\_{q}) &= \frac{L\_{d} L\_{q} \omega\_{e}^{2} I\_{d} - R\_{c} L\_{d} \omega\_{d} I\_{q} + L\_{q} \mathcal{V}\_{\rm md} \omega\_{e}^{2}}{R\_{c}^{2} + L\_{d} L\_{q} \omega\_{e}^{2}} \\ I\_{\rm cq}(I\_{d}, I\_{q}) &= \frac{R\_{c} L\_{d} \omega\_{e} I\_{d} + L\_{d} L\_{q} \omega\_{e}^{2} I\_{q} + R\_{c} \mathcal{V}\_{\rm md} \omega\_{e}}{R\_{c}^{2} + L\_{d} L\_{q} \omega\_{e}^{2}} \end{split} \tag{1.3}$$

Then the PMSM losses can be expressed in terms of *Id* and *Iq* as well (14).

$$\begin{aligned} P\_{\text{loss, PMSM}} &= M\_d^2 + B I\_d^2 + C I\_d + D I\_d + E I\_d I\_q + F\\ A &= \frac{L\_d^2 L\_d^2 (R\_s + R\_e) + R\_e^2 (R\_s R\_e^2 + 2 R\_s L\_d \omega\_d^2 + R\_e L\_d^2 \omega\_e^2)}{\left(R\_s^2 + L\_d L\_d \omega\_d^2\right)^2} \\ B &= \frac{L\_d^2 L\_q^2 \omega\_e^4 (R\_s + R\_e) + R\_e^2 \left(R\_s R\_e^2 + 2 R\_s L\_d \omega\_d \omega\_e^2 + R\_e L\_q^2 \omega\_e^2\right)}{\left(R\_s^2 + L\_d L\_d \omega\_d^2\right)^2} \\ C &= \frac{2 R\_s L\_d W\_{\text{max}} \omega\_e^2 \left(L\_d^2 \omega\_d^2 + R\_e^2\right)}{\left(R\_s^2 + L\_d L\_d \omega\_d^2\right)^2} \\ D &= \frac{2 R\_s^2 L\_q W\_{\text{max}} \omega\_e^2 \left(L\_d - L\_q\right)}{\left(R\_s^2 + L\_d L\_q \omega\_e^2\right)^2} \\ E &= \frac{2 R\_s^2 L\_d L\_q \omega\_q^2 \left(L\_d - L\_q\right)}{\left(R\_s^2 + L\_d L\_q \omega\_e^2\right)^2} \\ F &= \frac{R\_s Q\_{\text{max}}^2 \left(L\_q^2 \omega\_q^2 + R\_e^2\right)}{\left(R\_s^2 + L\_d L\_q \omega\_e^2\right)^2} \end{aligned} \tag{14}$$

If we consider Equation (15) in Equation (14), then Equation (16) is obtained.

$$\begin{aligned} I\_d &= I\_s \cos \alpha \\ I\_q &= I\_s \sin \alpha \end{aligned} \tag{15}$$

$$P\_{\rm loss, PMSM} = A l\_s^2 \cos^2 a + B l\_s^2 \sin^2 a + C l\_s \cos a + D l\_s \sin a + E l\_s^2 \sin a \cos a + F \tag{16}$$

The PMSM losses are minimized by Equation (17), yielding Equation (18).

$$\frac{dP\_{\text{loss,PMSM}}}{d\alpha} = 0\tag{17}$$

$$\begin{array}{l} -2Al\_s^2 \sin \alpha \cos \alpha + 2BI\_s^2 \sin \alpha \cos \alpha - Cl\_s \sin \alpha + \\ + DI\_s \cos \alpha + EI\_s^2 \cos^2 \alpha - EI\_s^2 \sin^2 \alpha = 0 \\ -2Al\_dI\_q + 2BI\_dI\_q - Cl\_q + DI\_d + EI\_d^2 - EI\_q^2 = 0 \end{array} \tag{18}$$

By rearranging Equation (18) we obtain Equation (19) and consequently Equation (20) which provides the necessary current *Id* for every corresponding combination of *Iq* and ω*e*.

$$0 = EI\_d^2 + \left(D + 2I\_q(B - A)\right)I\_d - \left(EI\_q^2 + CI\_q\right) \tag{19}$$

$$I\_{d1,2} = \frac{-D + 2I\_q(A - B)}{2E} \pm \frac{\sqrt{\left(D + 2I\_q(B - A)\right)^2 + 4E\left(EI\_q^2 + CI\_q\right)}}{2E} \tag{20}$$

Figure 5 presents the obtained MTPA characteristic by Equation (8), and the obtained ME characteristics by Equation (20) at different values of rotor speed and different values of the iron core resistance *Rc*. Although, this may not be evident from Figure 5 and similarly to the MTPA characteristics, the ME characteristics are nonlinear. It is however evident from Figure 5 that the iron core resistance *Rc* plays a crucial role in shaping the ME characteristic. Unfortunately, in practice its value changes depending on operating conditions and it cannot be determined directly or to be given in form of a function.

**Figure 5.** Calculated MTPA and ME characteristics at different rotor speeds and different constant values of the iron core resistance *Rc.*

Therefore, the authors have turned to an implicit but general approach for determination of the ME characteristic based on experiments, which is described hereinafter. The determined ME characteristics are then used to improve the efficiency of the PMSM drive in sensorless control as well as for the control realization with position sensor. The approach proposed here is novel, and can easily be adopted by control designers when of-the-shelf (or unknown) PMSMs are employed—contrarily to methods based, e.g., on finite element analyses which depend on detailed information about the employed geometry and materials in the PMSM.

#### *3.3. Experimental Determination of ME Characteristics*

The ME characteristics derived theoretically in the previous subsections are based on the steady-state PMSM model with constant parameters. With the aim to improve the PMSM drive efficiency and to evaluate the previously presented results, the characteristics were determined also by measurements performed on the PMSM drive, schematically shown in Figure 6. The presented approach includes all effects that are also present during real drive operation in the ME characteristics, these being effects of magnetic nonlinearities, saturation, cross-coupling, end-winding leakage, permanent magnets and slotting, converter nonlinearities and dead times in transistor switching. The PMSM was torque controlled and fed by the inverter. This means that during all tests the d- and q-axis currents were close-loop controlled, following their references set to different constant values. The vector control of PMSM was realized on the control system dSPACE PPC 1103. The rotor speed was closed-loop controlled by an active load IM. The input and output power of the inverter and PMSM, and their efficiencies, were measured by the power analyzers NORMA I (Norma D5255M) and NORMA II (Norma D6100), as shown Figure 6. The equal equipment was used throughout all successive tests, thus enabling a direct and comprehensive comparison of the PMSM drive efficiency operating in closed-loop control with or without the position sensor.

**Figure 6.** Experimental set-up—schematic presentation.

The ME characteristics were determined by changing the length and angle of the stator current vector at different values of the rotor speed. In individual steady-states, based on the measurements of voltages, currents, powers, torque and speed, the efficiencies are determined. A three-dimensional matrix of measurement results was generated, which includes length and angle of the stator current vector as well as the rotation speed. The experimentally obtained drive efficiency curves at 2250 rpm in dependence of the length and angle of the stator current vector are shown in Figure 7. The ME characteristic at 2250 rpm was then determined by connecting the points of maximum efficiencies of respective efficiency curves (black line). The determined ME characteristics at different rotor speeds are then shown in Figure 8, where they are presented in dependence of d- and q-axis currents. From Figure 8 it is evident that if we want to control the drive with maximum efficiency, the d-axis current component has to be increased by increase of the rotor speed, which is also known as flux-weakening. Unfortunately, due to hardware limitations, only currents of up to 6 A could be achieved during the experiments.

**Figure 7.** ME characteristic at 2250 rpm in dependence of the length and angle of the stator current vector.

**Figure 8.** Experimentally determined ME characteristics.

#### **4. Results and Discussion**

This section shows how the employment of the previously experimentally determined ME characteristics influences the efficiency of the discussed speed closed-loop controlled PMSM drive, when it is used to generate the d- and q-axis current references. The rotor position and speed required for the control realization are obtained in two different ways: firstly, by using sensors and secondly, by the sensorless control. In the last case, the measured PMSM currents and DC-bus voltage are used together with the PMSM dynamic model to calculate the rotor position and speed. The active flux based sensorless control method is applied. The active flux Ψ*AF* is defined as a hypothetical flux, which multiplied by the instantaneous q-axis current *iq* produces the measured instantaneous electro-magnetic torque *te* [32–35]. Consequently, the effects which are neglected in the PMSM model are thus considered by employment of Ψ*AF*.

$$t\_c = p\_p \big| \Psi\_{AF} \big| i\_q \tag{21}$$

The active flux concept uses the flux linkage vector estimates **Ψ**ˆ αβ*<sup>u</sup>* and **Ψ**ˆ αβ*<sup>i</sup>* obtained by the voltage (index *u*) and current models (index *i*) (Equation (22)), written in the α-β reference frame [32]. The difference between both estimates gives the compensation voltage *uk*, which is added in the voltage model to eliminate the DC-offset and other errors in measured stator currents, variation of stator resistance, and integrator drifts in the calculation procedure.

$$\begin{aligned} \Psi\_{\alpha\beta\iota} &= \mathop{\rm t}\limits\_{0}^{t} \Big( \mathbf{u}\_{\alpha\beta}(\tau) - R\_{\mathbf{s}} \mathbf{i}\_{\alpha\beta}(\tau) + \boldsymbol{u}\_{k}(\tau) \Big) d\tau + \Psi\_{\alpha\beta\iota}(0) \\ \Psi\_{\alpha\beta i} &= \mathop{\rm L}\limits\_{\alpha\beta} \Big[ \begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix} \Big] \mathbf{i}\_{\alpha\beta} + \Psi\_{md} \Big[ \begin{matrix} \cos\theta\_{\varepsilon} \\ \sin\theta\_{\varepsilon} \end{matrix} \Big] \end{aligned} \tag{22}$$

Considering zero initial conditions, the required estimates of the electrical rotor speed ωˆ*<sup>e</sup>* and position θˆ*<sup>e</sup>* are calculated by the nonlinear PLL-observer (Equation (23)) [31], where ε represents the error signal; *k*<sup>1</sup> and *k*<sup>2</sup> represent the observer gains.

$$\begin{aligned} \label{eq1} \phi\_{\mathfrak{c}} &= \int\_{0}^{t} k\_{1} \varepsilon(\tau) d\tau\\ \begin{array}{c} \partial\_{\mathfrak{c}} = \int\_{0}^{t} (\partial\_{\mathfrak{c}}(\tau) + k\_{2} \varepsilon(\tau)) d\tau \end{array} \end{aligned} \tag{23}$$

Thus, in the case of sensorless control, the estimates ωˆ*<sup>e</sup>* and θˆ*<sup>e</sup>* (Equation (23)) are used instead measured speed and position in the control algorithm, enabling closed loop-speed control without mechanical speed sensor.

Figure 9 shows the schematic presentation of the speed closed-loop control realizations with sensors (position and speed) and without them. These control realizations are applied in the experimental set-up, shown in Figure 6, which is used to measure drive efficiencies in cases when *idr* = 0, MTPA and ME characteristics, respectively. These characteristics (MTPA or ME) are used to determine the d- and q-axis current references *idr* and *iqr*.

Figure 10 shows the measured PMSM drive efficiencies (ηALL) at the rated speed of 2250 rpm and different loads in the case when *idr* = 0 and in the cases when *idr* and *iqr* are determined with the MTPA and ME characteristics, respectively. The results are given for the control realization with position and speed sensors and for the sensorless control. They are completed by the presentation of differences between different measured efficiency curves.

**Figure 9.** Schematic presentation of the applied speed closed-loop control.

**Figure 10.** Measured efficiencies of the speed closed-loop controlled PMSM drive and their differences for the control realization with the mechanical sensor and for the sensorless control in the cases when *idr* and *iqr* are determined with MTPA and ME characteristics and in the case when *idr* = 0.

The presented results clearly show that the lowest efficiencies are achieved when *idr* = 0. The drive efficiency is improved when the MTPA characteristic is used for the generation of *idr* and *iqr*. When it is replaced with the ME characteristic, the drive efficiency is further improved. Similar findings have been reported in [21,23,28] for an IPMSM drive with position sensor. The efficiency improvements decrease with the increase of the load torque (*Tm*). This is true for the control realization with the position and speed sensors as well as for the sensorless control realization.

#### **5. Conclusions**

This paper presents a way of improving the PMSM drive efficiency, where the more usual *idr* = 0 and MTPA characteristic for generation of reference current in vector controlled PMSM drives are compared with the proposed ME characteristic. There, both copper and iron core losses of the PMSM are considered. The equations which relate the ME characteristics to the measurable stator currents *Id* and *Iq* which are needed by control designers (instead to the magnetizing currents *Imd* and *Imq*) are derived, and consequently the equation which provides the optimal current *Id* for a given load torque by considering the ME characteristic is obtained for the first time. The problem of adequately determining the iron core resistance parameter value *Rc* is effectively mitigated by employing experimentally determined ME characteristics on the real PMSM drive system. Differences between the MTPA and ME characteristics are detailed; they are used in speed closed-loop vector control realizations with position and speed sensors and without them. The ME characteristic is implemented in sensorless control for the first time, where a nonlinear PLL observer is applied. Its inputs are calculated by the modified active flux method. The work clearly demonstrates that the use of ME characteristic instead of MTPA characteristic or *idr* = 0 in the current references generation improves the efficiency of PMSM drive realizations with position and speed sensors and without them—sensorless control.

This particular comparative case study showed that choosing an adequate switching frequency might improve the drive efficiency by 1.5–5.0%. At small load and speed values where the voltage value is small as well, an adequate DC-bus voltage control may improve the drive efficiency up to 30%. Ultimately, this particular PMSM drive efficiency was improved by 0.1–0.5% when the proposed ME was implemented in comparison to the MTPA, and by 1.0–2.0% in comparison to *idr* = 0, respectively. Therefore, it should be preferred that the PMSM drive control solutions proposed here are implemented in applications where efficiency is crucial. Moreover, the ME implementation can cut down the annual energy consumption of the studied drive by up to approximately 130 kWh, depending on the load and control case presented in Figure 11. This fact certainly cannot be neglected, since its implementation requires merely the modification of a few lines of code in the control realization.

**Figure 11.** Annual energy savings achieved by different methods of reference current generation in the speed closed-loop controlled PMSM drive for the realization with the mechanical sensor and for the sensorless control.

**Author Contributions:** Ž.P. performed derivations, experimental work and analyses. T.M. and M.B. prepared the final form of the paper, including editing of the text, equations and figures. G.Š. prepared the experiments and coordinated the research work and paper preparation.

**Funding:** This work was supported by the Slovenian Research Agency [grant numbers J2-1742, P2-0115].

**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/).

## *Article* **Induction Machine Control for a Wide Range of Drive Requirements**

#### **Bojan Grˇcar 1,\*, Anton Hofer <sup>2</sup> and Gorazd Štumberger <sup>1</sup>**


Received: 24 October 2019; Accepted: 23 December 2019; Published: 31 December 2019

**Abstract:** In this paper, a method for induction machine (IM) torque/speed tracking control derived from the 3-D non-holonomic integrator including drift terms is proposed. The proposition builds on a previous result derived in the form of a single loop non-linear state controller providing implicit rotor flux linkage vector tracking. This concept was appropriate only for piecewise constant references and assured minimal norm of the stator current vector during steady-states. The extended proposition introduces a second control loop for the rotor flux linkage vector magnitude that can be either constant, programmed, or optimized to achieve either maximum torque per amp ratio or high dynamic response. It should be emphasized that the same structure of the controller can be used either for torque control or for speed control. Additionally, it turns out that the proposed controller can be easily adapted to meet different objectives posed on the drive system. The introduced control concept assures stability of the closed loop system and significantly improves tracking performance for bounded but arbitrary torque/speed references. Moreover, the singularity problem near zero rotor flux linkage vector length is easily avoided. The presented analyses include nonlinear effects due to magnetic saturation. The overall IM control scheme includes cascaded high-gain current controllers based on measured electrical and mechanical quantities together with a rotor flux linkage vector estimator. Simulation and experimental results illustrate the main characteristics of the proposed control.

**Keywords:** induction machines; nonlinear control; energy efficiency; torque/speed control

#### **1. Introduction**

Squirrel cage induction machines (IM) are widely used in various industrial applications, machining, transportation, and electrical vehicles. They are relatively simply constructed and therefore easy to produce but quite difficult to control, especially if high performance and energy efficiency are required. Field oriented control (FOC) and direct torque control (DTC) schemes with many modifications introduced to fulfill different control objectives are widely accepted by the industry. Following the progress in non-linear control theory, several propositions based on feedback linearization, passivity, and flatness [1–4] were introduced in the past. These concepts facilitate deeper understanding of nonlinear IM dynamics and improve feedback performance and robustness but are often quite difficult to design and to implement. As a consequence, a great number of modifications and extensions have been proposed for basic FOC and DTC to improve overall performance, efficiency, and robustness along with new estimation techniques [5–10].

It is a long term goal to develop IM control strategies which are able to meet increasing efficiency requirements, i.e., controllers which can produce prescribed torque trajectories while minimizing the power losses. For example, in [11] an optimal torque controller for an IM with a linear magnetic characteristic was proposed based on the calculus of variations. A quite different approach was presented in [12]. Introducing a general *steady state* model suitable for induction motors, permanent-magnet synchronous motors, reluctance synchronous motors, and DC motors, the problem of power loss minimization including core saturation is considered. In [13] a control strategy is proposed based on ideas similar to [11]. The problem of missing information about the future reference torque trajectory is removed assuming simple first order transients between constant reference torque phases. The online implementation requires the solution of a small static optimization problem. The problem of power loss minimization for step changes of the torque reference trajectory is also considered in [14]. In addition to [13] saturation effects of the main inductance are taken into account and a simplified suboptimal solution avoiding an online optimization is presented. Recently a quite different control scheme for a permanent magnet synchronous machine was given in [15]. Based on an linearized parameter dependent mathematical model for the motor which also includes magnetic saturation, a model predictive control concept is developed. The online optimization problem is solved by the projected fast gradient method.

In this paper, a new control concept based on 3-D non-holonomic (NH) integrator including drift terms is proposed. Using standard modeling assumptions and assuming "perfect" current control, a reduced IM (current) model can be written in the general NH integrator form. From the control standpoint NH systems, which originate in mechanics, are quite challenging since they fail to fulfill the necessary conditions for stabilization with a smooth state feedback [16]. For this class of systems, control schemes based on discontinuous feedback (hybrid, sliding mode, and time optimal), time-varying state feedback (trajectory planning, back-stepping, pattern generation) were proposed in [17–20]. In [21] the authors introduced a non-linear state controller based on NH system analysis applied to IM torque control. Considering the non-holonomic constraint requiring periodic orbits of the rotor flux linkage and stator current vectors, the proposed feedback assured global asymptotic stability and maximal torque/amp ratio in steady states for the nominal parameter case. An essential characteristic of the proposed control scheme is implicit rotor flux linkage tracking provided by the adjusted amplitude and frequency modulation of the stator current vector. Two main limitations were observed for the proposed control: only piecewise *constant* references were allowed and the problem of singularity (zero rotor flux) required an additional time optimal flux controller and a switching logic to provide soft transitions between the two regimes [22].

This contribution is based on the results given in [23] for the control of the general 3-D-nonholonomic integrator with drift. The main focus of the present paper is the adaption of the control concept presented in [23] for the usage in IM torque and speed control applications. For this purpose a nonlinear flux model is included and the problem of flux optimization in order to achieve the maximum torque per ampere feature is considered from the practical point of view. It is shown that a simple suboptimal strategy for the adjustent of the flux magnitude can be applied to a wide range of torque reference scenarios as well as to speed control tasks. The control input is consequently composed of two orthogonal components, the first providing the required rotor flux linkage vector and the second producing the desired torque. In this way the overall control structure becomes partially similar to FOC but without relying on the uncertain rotor flux reference frame. The proposed structure enables the rotor flux linkage vector to be either constant, programmed, or optimized for energy efficiency with respect to drive characteristics and requirements. Due to similarities with widespread FOC, the proposed control scheme can be easily implemented on a standard industrial hardware.

The presentation of the results is structured as follows. In Section 2, the basic IM model is introduced in the rotor reference frame. The resulting reduced current model along with a dynamic torque estimation constitute a general 3-D non-holonomic integrator that is used to introduce the basic nonlinear state controller. Section 3 shows how the basic control structure can be extended in order to provide improved tracking performance. A proof of the stability of the resulting closed loop system is given. In Section 4 the mathematical description of the IM is extended by a simple nonlinear flux model covering the phenomenon of magnetic saturation. Next, the problem of rotor flux magnitude optimization is addressed, such that torque tracking of the IM is performed efficiently. Different classes of torque reference signals are considered and appropriate solutions are given. Section 5 is dedicated to extensions of the controller structure and the task of speed control. In Section 6, implementation aspects and experimental results confirming the improved tracking performance for torque and speed control are discussed. Concluding remarks are given in Section 7.

#### **2. Im Model and Basic Controller**

Using a standard modeling approach the two-axis dynamic model of a 2-pole IM model written in an arbitrary *d*-*q* reference frame is given by

$$\begin{aligned} \dot{\mathbf{i}}\_{d,q} &= \left(-\tau\_{\sigma}^{-1}\mathbf{I} - \mathbf{J}\omega\_{\mathbf{i}}\right)\mathbf{i}\_{d,q} + \left(\frac{\mathcal{M}}{L\_{\tau}L\_{\sigma}\tau\_{\mathbf{r}}}\mathbf{I} - \frac{\omega\_{m}\mathcal{M}}{L\_{\tau}L\_{\sigma}}\mathbf{J}\right)\Psi\_{d,q} \\ &+ L\_{\sigma}^{-1}\mathbf{u}\_{d,q} \\ \dot{\Psi}\_{d,q} &= \left(-\tau\_{\tau}^{-1}\mathbf{I} - \mathbf{J}(\omega\_{s} - \omega\_{m})\right)\Psi\_{d,q} + \frac{\mathcal{M}}{\tau\_{\tau}}\mathbf{i}\_{d,q} \\ T\_{el} &= \frac{\mathcal{M}}{L\_{\tau}}\left(\psi\_{d}i\_{q} - \psi\_{q}i\_{d}\right), \end{aligned} \tag{1}$$

where **i***d*,*q*, **u***d*,*q*, *ψd*,*<sup>q</sup>* are the stator current vector, the stator voltage vector, and the rotor flux linkage vector written in the *d*-*q* reference frame. The mutual inductance is denoted by *M*, *Lr* is the rotor inductance, *Rr* is the rotor resistance, *ω<sup>m</sup>* is the rotor speed, *ω<sup>s</sup>* is the speed of the *d*-*q* reference frame with respect to the stator, *Tel* is the electric torque while **I** is the 2 × 2 identity matrix and **J** is the 2 × 2 skew symmetric rotational matrix. The leakage inductance *Lσ*, the resistance *Rσ*, the rotor time constant *<sup>τ</sup><sup>r</sup>* and the leakage time constant *τσ* are defined as *<sup>L</sup><sup>σ</sup>* = *Ls* − *<sup>M</sup>*2/*Lr*, *<sup>R</sup><sup>σ</sup>* = *Rs* − (*RrM*)2/*L*<sup>2</sup> *r*, *τ<sup>r</sup>* = *Lr*/*Rr* and *τσ* = *Lσ*/*Rσ*. Throughout the text references, estimated variables and errors values are denoted as (.)∗,( \$.),( %.).

Under the assumption of a perfect current control, i.e., **i***d*,*<sup>q</sup>* = **i** ∗ *<sup>d</sup>*,*q*, the simplified IM (current) model written in the rotor reference frame *γ*-*δ* is obtained from the model written in (1) by considering *ω<sup>m</sup>* = *ω<sup>s</sup>* and changing *d* to *γ* and *q* to *δ*:

$$
\dot{\Psi}\_{\gamma,\delta} = -\mathbf{r}\_r^{-1} \Psi\_{\gamma,\delta} + \frac{M}{\mathbf{r}\_r} \mathbf{i}\_{\gamma,\delta} \tag{2}
$$

where **i***γ*,*<sup>δ</sup>* are control inputs, while the machine torque represents the controlled variable. For torque control the output error is defined as *T*˜ *el* := *T*<sup>∗</sup> *el* <sup>−</sup> *<sup>T</sup>*<sup>ˆ</sup> *el*, where the actual machine torque *Tel* in (1) is replaced by a filtered estimate

$$
\tau\_f \dot{\Upsilon}\_{el} = -\Upsilon\_{el} + \frac{M}{L\_I} \left( \psi\_\gamma i\_\delta - \psi\_\delta i\_\gamma \right). \tag{3}
$$

In this model *τ<sup>f</sup> τ<sup>r</sup>* is the filter time constant, that is introduced to capture the estimation dynamics and to remove higher frequency components generated by estimation process and measurement of stator currents. A standard singular perturbation argument involving the estimation time constant *τ<sup>f</sup>* recovers the algebraic expression for the actual machine torque. Equations (2) and (3) constitute the general 3-D non-holonomic integrator that includes drift terms.

Using the output error *T*˜ *el* and the abbreviation *s*<sup>3</sup> := sign(*T*<sup>∗</sup> *el*) with the convention *T*<sup>∗</sup> *el* = 0 → *s*<sup>3</sup> = 0, the basic nonlinear state controller was proposed in [21] as:

$$
\begin{bmatrix} i\_{\gamma} \\ i\_{\delta} \end{bmatrix} = \frac{\tau\_r}{\mathcal{M}} \begin{bmatrix} \frac{1}{\tau\_r} + s\_3 k\_p \dot{T}\_{el} & -\left(\frac{s\_3}{\tau\_r} + k\_p \dot{T}\_{el}\right) \\\ \frac{s\_3}{\tau\_r} + k\_p \dot{T}\_{el} & \frac{1}{\tau\_r} + s\_3 k\_p \dot{T}\_{el} \end{bmatrix} \begin{bmatrix} \psi\_{\gamma} \\ \Psi\_{\delta} \end{bmatrix} \tag{4}
$$

where *kp* > 0 denotes a controller gain. The terms that are output error independent compensate the drift terms and provide purely oscillatory dynamics of the rotor flux linkage vector *ψγ*,*<sup>δ</sup>* as the output error *T*˜ *el* vanishes. On the other hand, the terms that are output error dependent provide the adjusted amplitude and frequency modulation of the input current vector. The feedback system is obtained from (2) to (4) in the following form:

$$\begin{aligned} \dot{\psi}\_{\gamma} &= s\_3 k\_p \mathcal{T}\_{el} \psi\_{\gamma} - \left( s\_3 \boldsymbol{\tau}\_r^{-1} + k\_p \boldsymbol{\mathcal{T}}\_{cl} \right) \psi\_{\delta} \\ \dot{\psi}\_{\delta} &= \left( s\_3 \boldsymbol{\tau}\_r^{-1} + k\_p \boldsymbol{\vec{\mathcal{T}}}\_{cl} \right) \psi\_{\gamma} + s\_3 k\_p \boldsymbol{\vec{\mathcal{T}}}\_{cl} \psi\_{\delta} \\ \dot{\boldsymbol{\mathcal{T}}}\_{cl} &= -\boldsymbol{\tau}\_f^{-1} \boldsymbol{\vec{\mathcal{T}}}\_{cl} + \boldsymbol{\tau}\_f^{-1} \frac{\boldsymbol{\tau}\_r}{L\_r} \left( s\_3 \boldsymbol{\tau}\_r^{-1} + k\_p \boldsymbol{\vec{\mathcal{T}}}\_{cl} \right) \left( \boldsymbol{\Psi}\_{\gamma}^2 + \boldsymbol{\Psi}\_{\delta}^2 \right) \end{aligned} \tag{5}$$

In [21] it was proven that (5) has the following properties: For any constant *T*∗ *el* = 0 , any *T*ˆ *el*(0) and *ψγ*,*δ*(0) = **0** the system converges to a state where *T*ˆ *el* = *T*<sup>∗</sup> *el* and the components of *ψγ*,*δ*(*t*) are harmonic functions with the amplitude

$$\|\|\psi\_{\gamma,\delta}\|\|\_{2} = \sqrt{L\_{r} \left|T\_{cl}^{\*}\right|}.$$

and the angular frequency *s*3*τ*−<sup>1</sup> *<sup>r</sup>* . Furthermore, a distinctive characteristic of (5) is that, for a given piecewise constant output reference *T*∗ *el* the minimal norm of the control input vector **i***γ*,*<sup>δ</sup>* = [ *i<sup>γ</sup> i δ* ] *T* is assured in steady states

$$\|\mathbf{i}\_{\gamma\mathcal{A}}\|\_2 = \sqrt{2\frac{L\_r}{M^2}\left|T\_{cl}^\*\right|}.$$

The angle between the rotor flux vector and the stator current is *s*3*π*/4 *permanently* as was already presented in [21]. The rotation speed of the stator current vector is obtained as *ω<sup>I</sup>* = *ω<sup>r</sup>* + *ωm*, where the introduced "relative" speed *ω<sup>r</sup>* corresponds to the inverse time rotor constant *τ*−<sup>1</sup> *<sup>r</sup>* . The relative speed exactly matches the slip speed in steady states. Further remarks and characteristics of the controlled system (5) regarding stability, robustness, singularity avoidance, and overall performance can be found in [21,22].

#### **3. Improved Torque Tracking Controller**

Since from now on all vectors are expressed as vectors in the rotor reference frame the index combination *γ*, *δ* is omitted. The main characteristic of the presented basic controller, namely the implicit rotor flux tracking without separate control loop works perfectly for a piecewise constant torque references assuring also maximal torque-per-amp ratio. The main obstacle for a good tracking performance in general is observed especially at zero crossings of *T*∗ *el* with finite slope since the magnitude of the rotor flux linkage becomes too small or even zero. This results in the output *T*ˆ *el* getting out of control for a short period as can easily be seen from the third equation in (5).

#### *3.1. Controller Modification*

To cope with this problem the basic controller scheme (4) has to be conceptually modified. Based on the ideas outlined in [23] the control input is composed of two orthogonal parts

$$\mathbf{i} = \mathbf{i}\_{\Psi} + \mathbf{i}\_{\tau} \tag{6}$$

where the vector **i***ψ* influences only the rotor flux linkage magnitude *ψ*, whereas the vector **i***τ* facilitates machine torque. Considering the basic control structure it is intended that the effect of the diagonal elements in (4) is replaced by an appropriately selected vector **i***ψ*. The presentation of the control design can be simplified if polar coordinates are introduced for the rotor flux vector

$$
\psi\_{\gamma} = \psi \cos \varphi \qquad \psi\_{\delta} = \psi \sin \varphi \tag{7}
$$

with

$$\psi = \left\| \Psi\_{\gamma,\delta} \right\|\_2 = \sqrt{\Psi\_{\gamma}^2 + \Psi\_{\delta}^2} \ \ \ \ \varphi = \arctan \left( \frac{\psi\_{\delta}}{\psi\_{\gamma}} \right) . \tag{8}$$

The plant Equations (2) and (3) are then obtained as:

$$\psi = -\tau\_r^{-1}\psi + \tau\_r^{-1}M \left( i\_\gamma \cos\varphi + i\_\delta \sin\varphi \right) \tag{9}$$

$$
\psi \phi = \pi\_r^{-1} M \left( -i\_\gamma \sin \varphi + i\_\delta \cos \varphi \right) \tag{10}
$$

$$\dot{\hat{T}}\_{cl} = -\tau\_f^{-1}\hat{T}\_{cl} + \tau\_f^{-1}\frac{M}{L\_r}\left(-i\_\gamma\sin\varphi + i\_\delta\cos\varphi\right)\psi\tag{11}$$

$$T\_{cl} = \frac{M}{L\_r} \left( -i\_\gamma \sin \varphi + i\_\delta \cos \varphi \right) \psi \tag{12}$$

Now it is obvious that the control vector **i***ψ* must be of the form

$$\mathbf{i}\_{\emptyset} = \begin{bmatrix} i\_{\emptyset \gamma} \\ i\_{\emptyset \delta} \end{bmatrix} = \begin{bmatrix} I \cos \varphi \\\ I \sin \varphi \end{bmatrix}$$

with a suitable selected magnitude *I* in order to influence the flux magnitude *ψ* only. Introducing the reference signal *ψ*∗ for the rotor flux linkage magnitude *ψ* a simple proportional control is proposed for the "magnetizing" vector **i***ψ*

$$\mathbf{i}\_{\psi} = \frac{1}{M} \begin{bmatrix} \left(\psi^\* + k\_{\psi}(\psi^\* - \psi)\right) \cos \varrho\\ \left(\psi^\* + k\_{\psi}(\psi^\* - \psi)\right) \sin \varrho \end{bmatrix} \tag{13}$$

where *k<sup>ψ</sup>* > 0 is a tuning parameter and *ψ*<sup>∗</sup> > 0 is assumed. By inserting (13) in (9) it follows that the rotor flux magnitude is *ψ* driven by the linear first order dynamics with unit gain:

$$
\dot{\Psi} = -\tau\_r^{-1} (1 + k\_{\Psi}) \psi + \tau\_r^{-1} (1 + k\_{\Psi}) \psi^\* \tag{14}
$$

The time constant of (14) can be simply influenced by an appropriate selection of the parameter *kψ*. It is important to note that (14) is completely decoupled from the other two states *ϕ*, *T*ˆ *el*.

Since **i***<sup>ψ</sup>* should replace the effect of the diagonal elements in (4) the second control vector **i***<sup>τ</sup>* influencing only the output error results from the remaining part of (4):

$$\begin{aligned} \mathbf{i}\_{\tau} &= \frac{\mathbf{r}\_{r}}{\mathcal{M}} \begin{bmatrix} - \left( \mathbf{s}\_{3} \mathbf{r}\_{r}^{-1} + k\_{p} \mathbf{\mathcal{T}}\_{cl} \right) \boldsymbol{\Psi}\_{\delta} \\ \left( \mathbf{s}\_{3} \mathbf{\mathcal{T}}\_{r}^{-1} + k\_{p} \mathbf{\mathcal{T}}\_{cl} \right) \boldsymbol{\Psi}\_{\boldsymbol{\mathcal{V}}} \end{bmatrix} \\ &= \frac{\mathbf{r}\_{r}}{\mathcal{M}} \begin{bmatrix} - \left( \mathbf{s}\_{3} \mathbf{\mathcal{T}}\_{r}^{-1} + k\_{p} \mathbf{\mathcal{T}}\_{cl} \right) \boldsymbol{\Psi} \sin \boldsymbol{\Psi} \\ \left( \mathbf{s}\_{3} \mathbf{\mathcal{T}}\_{r}^{-1} + k\_{p} \mathbf{\mathcal{T}}\_{cl} \right) \boldsymbol{\Psi} \cos \boldsymbol{\Psi} \end{bmatrix} \end{aligned} \tag{15}$$

Obviously **i***τ* has no effect on the flux magnitude *ψ*. By inserting (15) into (10), it follows

$$
\phi = s\_3 \pi\_r^{-1} + k\_p \mathcal{T}\_{el} \dots
$$

When the torque error *T*˜ *el* vanishes the rotor flux linkage dynamics results in a harmonic oscillator with an angular speed *ϕ*˙ = *s*3*τ*−<sup>1</sup> *<sup>r</sup>* . This property which is implicitly contained in the basic controller (4) is reasonable only for piecewise *constant* references *T*∗ *el* but it is too restrictive in general. Therefore, it is proposed to replace the term *s*3*τ*−<sup>1</sup> *<sup>r</sup>* in (15) with the expression

$$R\_{\mathcal{I}} \frac{T\_{cl}^{\*}}{\left(\psi^{\*}\right)^{2}} \quad \text{(with the restriction } \psi^{\*} > 0\text{)}\text{.}$$

*Energies* **2020**, *13*, 175

One of the reasons for this replacement is the fact that in the case of a constant reference *T*∗ *el* and the choice *ψ*∗ = " *Lr T*<sup>∗</sup> *el* the same steady state is obtained as by using the basic controller (4). This means that the important property of (4), i.e., producing a constant torque with minimum norm of the control input is recovered by the modified controller. Finally, by combining (13) and (15) (with the proposed modification), a new control law is obtained:

$$\begin{split} \dot{i}\_{\mathcal{V}} &= \frac{\mathsf{T}\_{\mathcal{r}}}{\mathcal{M}} \Big[ - \left( R\_{\mathcal{r}} \frac{\mathsf{T}\_{\mathcal{cl}}^{\*}}{(\mathsf{\psi}^{\*})^{2}} + k\_{\mathcal{P}} \tilde{\mathsf{T}}\_{\mathcal{cl}} \right) \boldsymbol{\psi} \sin \boldsymbol{\varrho}. \\ & \quad + \mathsf{r}\_{\mathcal{r}}^{-1} \left( \boldsymbol{\psi}^{\*} + k\_{\boldsymbol{\Psi}} (\boldsymbol{\psi}^{\*} - \boldsymbol{\psi}) \right) \cos \boldsymbol{\varrho} \Big] \\ \dot{i}\_{\mathcal{\delta}} &= \frac{\mathsf{T}\_{\mathcal{r}}}{\mathcal{M}} \Big[ \left( R\_{\mathcal{r}} \frac{\mathsf{T}\_{\mathcal{cl}}^{\*}}{(\mathsf{\psi}^{\*})^{2}} + k\_{\mathcal{P}} \tilde{\mathsf{T}}\_{\mathcal{cl}} \right) \boldsymbol{\psi} \cos \boldsymbol{\varrho}. \\ & \quad + \mathsf{r}\_{\mathcal{r}}^{-1} \left( \boldsymbol{\psi}^{\*} + k\_{\boldsymbol{\Psi}} (\boldsymbol{\psi}^{\*} - \boldsymbol{\psi}) \right) \sin \boldsymbol{\varrho} \Big] \end{split} \tag{16}$$

The resulting closed loop system built up by (9)–(11), (16) is given as:

$$\begin{aligned} \dot{\psi} &= -\tau\_r^{-1} (1 + k\_{\Psi}) \psi + \tau\_r^{-1} (1 + k\_{\Psi}) \psi^\* \\ \dot{\varphi} &= R\_r \frac{T\_{el}^\*}{\left(\psi^\*\right)^2} + k\_p T\_{el} \\ \dot{T}\_{el} &= -\tau\_f^{-1} (1 + \frac{k\_p}{R\_r} \psi^2) \hat{T}\_{el} + \tau\_f^{-1} (\frac{\psi^2}{\left(\psi^\*\right)^2} + \frac{k\_p}{R\_r} \psi^2) T\_{el}^\* \end{aligned} \tag{17}$$

Of course the control law (16) can be transformed back into the original state variable form by the inverse relations (7) and (8). Which, in combination with the plant model (2) and (3) leads to the following description of the closed loop system:

$$\begin{aligned} \psi\_{\gamma} &= -\tau\_{r}^{-1} \Psi\_{\gamma} - \left( R\_{r} \frac{T\_{cl}^{\*}}{\left(\Psi^{\*}\right)^{2}} + k\_{p} \mathcal{T}\_{cl} \right) \Psi\_{\delta} \\ &+ \tau\_{r}^{-1} \left( \psi^{\*} + k \varphi \left( \psi^{\*} - \sqrt{\Psi\_{\gamma}^{2} + \Psi\_{\delta}^{2}} \right) \right) \frac{\Psi\_{\gamma}}{\sqrt{\Psi\_{\gamma}^{2} + \Psi\_{\delta}^{2}}} \\ \psi\_{\delta} &= -\tau\_{r}^{-1} \Psi\_{\delta} - \left( R\_{r} \frac{T\_{cl}^{\*}}{\left(\Psi^{\*}\right)^{2}} + k\_{p} T\_{cl} \right) \Psi\_{\gamma} \\ &+ \tau\_{r}^{-1} \left( \psi^{\*} + k \varphi \left( \psi^{\*} - \sqrt{\Psi\_{\gamma}^{2} + \Psi\_{\delta}^{2}} \right) \right) \frac{\Psi\_{\delta}}{\sqrt{\Psi\_{\gamma}^{2} + \Psi\_{\delta}^{2}}} \\ \dot{\mathcal{T}}\_{cl} &= -\tau\_{f}^{-1} \mathcal{T}\_{cl} + \frac{\tau\_{f}^{-1}}{R\_{r}} \left( R\_{r} \frac{T\_{cl}^{\*}}{\left(\Psi^{\*}\right)^{2}} + k\_{p} T\_{cl} \right) \left( \Psi\_{\gamma}^{2} + \Psi\_{\delta}^{2} \right) \end{aligned} (18)$$

#### *3.2. Stability Analysis*

Since the closed loop system has two external inputs it must be proven that the norm of the state vector **x** = ! *ψγ ψδ T*ˆ *el <sup>T</sup>* remains bounded whenever the norm of the input vector defined as **r** := ! *ψ*∗ *T*∗ *el <sup>T</sup>* and **x**(**0**)<sup>2</sup> are bounded. This task becomes much easier, if the model (17) is used. Obviously we have

$$\|\mathbf{x}\|\_2^2 = \boldsymbol{\psi}^2 + \boldsymbol{\Upsilon}\_{el'}^2 \tag{19}$$

therefore only the first and the third equation of (17) have to be considered. Next, it is assumed that the following bounds hold for the external inputs

$$0 < \psi\_{\text{min}}^{\*} \le \psi^{\*}(t) \le \psi\_{\text{max}}^{\*} \text{ for all } t \ge 0 \tag{20}$$

$$|T\_{cl}^\*(t)| \le T\_{cl, \text{max}}^\* \text{ for all } t \ge 0,\tag{21}$$

where *T*∗ *el*,max, *ψ*<sup>∗</sup> min, *ψ*<sup>∗</sup> max are some finite positive constants. With the restriction *k<sup>ψ</sup>* > 0 it follows immediately from the first equation in (17) that

$$
\psi(t) \le e^{-\tau\_r^{-1}(1+k\_\psi)t} \psi(0) + \psi^\*\_{\text{max}} \le \psi(0) + \psi^\*\_{\text{max}}.\tag{22}
$$

Thus *ψ*(*t*) remains bounded for any finite initial value *ψ*(0) ≥ 0 and any external input *ψ*∗(*t*) satisfying (20).

Now the third equation in (17) can be regarded as a linear *timevariant* first order system with the external input *T*∗ *el* and with the output *<sup>T</sup>*<sup>ˆ</sup> *el*. For the stability analysis of this system the homogeneous case (i.e., *T*∗ *el* = 0) is considered first which leads to:

$$\begin{aligned} \mathcal{T}\_{cl}(t) &= \mathcal{e}^{\int\_{t\_0}^t -\tau\_f^{-1}(1 + \frac{k\_p}{\mathcal{R}\_r}\psi^2)d\tau} \mathcal{T}\_{cl}(t\_0) \\ &= \mathcal{e}^{-\tau\_f^{-1}(t - t\_0)} \mathcal{e}^{\int\_{t\_0}^t -\tau\_f^{-1}\frac{k\_p}{\mathcal{R}\_r}\psi^2 d\tau} \mathcal{T}\_{cl}(t\_0) \end{aligned}$$

Taking into account *kp* > 0, it follows

$$\left| \left| \mathcal{T}\_{cl}(t) \right| \right| \leq e^{-\tau\_{f}^{-1}(t - t\_{0})} \left| \mathcal{T}\_{cl}(t\_{0}) \right| \cdot $$

which shows that the homogeneous system is uniformly exponentially stable. Furthermore from (20) to (22) we can derive the result

$$\pi\_f^{-1} \left| \frac{\psi^2}{\left(\psi^\*\right)^2} + \frac{k\_p}{R\_r} \psi^2 \right| \le c \text{ for all } t \ge 0$$

where *c* denotes some finite positive constant. This inequality combined with the fact that the homogeneous part is uniformly asymptotically stable assures that <sup>|</sup>*T*<sup>ˆ</sup> *el*(*t*)| remains bounded, whenever *T*<sup>∗</sup> *el*(*t*) and *<sup>T</sup>*<sup>ˆ</sup> *el*(0) are bounded [24]. Finally from (19) the desired result follows that **x**<sup>2</sup> remains bounded. Furthermore, it is important to note that also the angular rotation speed *ϕ*˙ of the rotor flux vector remains bounded (see second equation in (17)).

#### **4. Magnetic Saturation and Flux Reference Optimization**

The new control law (16) provides the flux reference signal *ψ*∗ as an additional external input. In this section the problem how to choose *ψ*∗ in order to achieve energy efficient tracking of a desired torque *T*∗ *el* will be tackled. Since we would like to take full advantage of the dynamic capabilities of the induction machine, magnetic saturation must be considered. A correct mathematical description of this phenomenon would lead to a very complicated model, useless for control design. Therefore, a simple approach proposed in [25] for field oriented control is adopted which is also applied in [26]. It is based on the magnetization curve of the induction machine given in the form

$$
\psi = f(i\_{\psi}),
\tag{23}
$$

where the magnetizing current *iψ* is defined as

$$\dot{a}\_{\psi} := \dot{\iota}\_{\gamma} \cos \varrho + \dot{\iota}\_{\delta} \sin \varrho. \tag{24}$$

It should be mentioned that *iψ* <sup>=</sup> & &**i***<sup>ψ</sup>* & & <sup>2</sup> is valid. Simple geometric considerations based on the phasors **i***γ*,*<sup>δ</sup>* and *ψγ*,*<sup>δ</sup>* reveal that this relation holds, where **i***<sup>ψ</sup>* is defined in (6) and (13) respectively. The magnetization curve can be obtained by measurements (e.g., [26]). Throughout the rest of the paper *M* and *Lr* are the nominal values of the mutual inductance and rotor inductance below saturation.

The following nonlinear flux model was introduced in [25]

$$
\psi = \pi\_r^{-1} M \left( -f^{-1}(\psi) + i\_{\psi} \right),
\tag{25}
$$

where *f* <sup>−</sup>1(*ψ*) is the inverse magnetizing function. For the next considerations it is useful also to introduce the torque producing current component

$$\dot{a}\_{\tau} := -\dot{\iota}\_{\gamma} \sin \varphi + \dot{\iota}\_{\delta} \cos \varphi. \tag{26}$$

(with the property |*iτ*| = **i***<sup>τ</sup>*<sup>2</sup>, given by the same geometric considerations as mentioned above), which allows us to write (12) in the form:

$$T\_{cl} = \frac{M}{L\_r} i\_\pi \psi \tag{27}$$

Note that the real electrical torque *Tel* is considered in the present case. The Equations (25) and (26) are the basis for the flux optimization problem. In the following subsection a piecewise constant torque reference *T*∗ *el* is considered and the corresponding optimization problem is solved. It is assumed that the desired torque can be produced by the IM within the given current and voltage limitations. The resulting optimal flux magnitude *ψopt* can be used as flux reference signal *ψ*∗ in the proposed control structure.

#### *4.1. Constant Reference T*∗ *El*

In the first case, the following question is considered: How should the *constant* flux *ψ* be selected, so that a desired *constant* torque *Tel* = *T*<sup>∗</sup> *el* = 0 is achieved with minimum norm of the current vector in steady state? In other words, we want to get the *maximum torque per amp* feature. If magnetic saturation is neglected (i.e., for the linear flux model) the solution is simple

$$
\psi\_{opt} = \sqrt{L\_r \left| T\_{el}^\* \right|} \tag{28}
$$

was already mentioned in Section 3. In the saturation case we get from (23)

$$i\_{\emptyset} = f^{-1}(\psi)\_{\prime}$$

which is the necessary value of *iψ* in order to achieve the constant *ψ* in steady state. From (27) it follows

$$i\_{\tau} = \frac{L\_{\tau}}{M} \frac{T\_{dl}^{\*}}{\psi}$$

and using the fact **<sup>i</sup>**<sup>2</sup> <sup>2</sup> = *i* 2 *<sup>ψ</sup>* + *i* 2 *<sup>τ</sup>* the problem can be formulated as:

$$\psi\_{opt} = \underset{\psi}{\text{argmin}} \left( \left( f^{-1}(\psi) \right)^2 + \left( \frac{L\_r}{M} \frac{T\_{cl}^\*}{\psi} \right)^2 \right) \tag{29}$$

By differentiating with respect to *ψ* the necessary condition is obtained

$$2f^{-1}(\psi)\frac{df^{-1}(\psi)}{d\psi} - \frac{2L\_r^2}{M^2}\frac{T\_{cl}^{\*2}}{\psi^3} \stackrel{!}{=} 0..$$

*Energies* **2020**, *13*, 175

Unfortunately, this equation cannot be solved for *ψ* but this causes no severe problems since it can be written as:

$$\mathcal{L}\_{\tau} \left| T\_{cl}^{\*} \right| = M \sqrt{\psi^3 f^{-1}(\psi) \frac{df^{-1}(\psi)}{d\psi}} =: \mathcal{g}(\psi) \tag{30}$$

The function *g*(*ψ*) can be evaluated and stored in a lookup table so that the desired result can be obtained numerically in the form

$$\psi\_{opt} = \text{g}^{-1}(L\_{\mathbb{F}} \left| T\_{el}^{\*} \right|). \tag{31}$$

A comparison between the optimal flux magnitude based on (28) and (31) is shown in Figure 1b. A similar result was presented in [14] where power losses for torque steps were considered and a function, which determines the optimal constant flux magnitude for constant torque levels, is given.

**Figure 1.** Magnetization curve (**a**); Optimal flux magnitude for constant torque (**b**).

#### *4.2. General Case*

In general, an arbitrary but bounded reference signal *T*∗ *el*(*t*) defined over a finite time interval 0 ≤ *t* ≤ *T*<sup>0</sup> is considered. In contrast to the forementioned case, a time function *ψ*(*t*) should be determined now, so that the criterion

$$E = \int\_0^{T\_0} \left( i\_{\psi}^2 + i\_{\tau}^2 \right) dt \tag{32}$$

is minimized for *Tel*(*t*) = *T*<sup>∗</sup> *el*(*t*). This means that the optimal flux trajectory *ψ*(*t*) for an arbitrary torque reference signal should be determined. Obviously the value of *E* will also depend on the initial condition *ψ*(0). In order to eliminate this influence, a periodic continuation of *T*∗ *el*(*t*) is assumed and the steady state solution of the problem is considered. In this case *ψ*(*t*) is a periodic time function too. From (27), it follows

$$i\_{\pi}(t) = \frac{L\_r}{M} \frac{T\_{el}^\*(t)}{\psi(t)}\tag{33}$$

and therefore the optimization problem can be formulated as:

$$\begin{aligned} \text{minimize} & \int\_{0}^{T\_{\text{\textquotedblleft}}} \left( i\_{\psi}(t)^{2} + \left( \frac{L\_{r}}{M} \frac{T\_{\text{cl}}^{\*}(t)}{\psi(t)} \right)^{2} \right) dt \\ \text{subject to} & \\ & \psi = \tau\_{r}^{-1} \mathcal{M} \left( -f^{-1}(\psi) + i\_{\Psi} \right) \\ \psi(0) = \psi(T\_{0}) \qquad i\_{\Psi}(0) = i\_{\Psi}(T\_{0}). \end{aligned} \tag{34}$$

The boundary conditions are given by the fact that *ψ*(*t*) and *iψ*(*t*) must also be continuous functions. This flux optimization problem can be treated by the classical methods of optimal control (see e.g., [27]). First the Hamilton function is introduced (the argument time is omitted for brevity)

$$H = i\_{\Psi}^{2} + \left(\frac{L\_{r}}{M} \frac{T\_{el}^{\*}}{\Psi}\right)^{2} + \lambda \tau\_{r}^{-1} M \left(-f^{-1}(\Psi) + i\_{\Psi}\right),\tag{35}$$

where *λ* denotes a Lagrange multiplier. Now the solution of the problem must satisfy the following necessary conditions:

$$0 = \frac{\partial H}{\partial \dot{i}\_{\Psi}} = 2\dot{i}\_{\Psi} + \lambda \pi\_r^{-1} M \quad \rightarrow \quad \dot{i}\_{\Psi} = -\frac{\lambda \pi\_r^{-1} M}{2} \tag{36}$$

$$\lambda = -\frac{\partial H}{\partial \psi} = 2\left(\frac{L\_r}{M}T\_{el}^\*\right)^2 \frac{1}{\Psi^3} + \lambda \,\tau\_r^{-1}M\frac{df^{-1}(\psi)}{d\psi} \tag{37}$$

$$\dot{\psi} = \frac{\partial H}{\partial \lambda} = \pi\_r^{-1} M \left( -f^{-1}(\psi) + i\_{\psi} \right) \tag{38}$$

Substituting *iψ* in Equation (38) with the expression from Equation (36), the resulting boundary value problem is obtained:

$$\begin{aligned} \lambda &= 2 \left( \frac{L\_r}{M} T\_{cl}^\* \right)^2 \frac{1}{\psi^3} + \lambda \pi\_r^{-1} M \frac{df^{-1}(\psi)}{d\psi} \\ \dot{\psi} &= \pi\_r^{-1} M \left( -f^{-1}(\psi) - \frac{\lambda \pi\_r^{-1} M}{2} \right) \\ \psi(0) &= \psi(T\_0) \qquad \lambda(0) = \lambda(T\_0) \end{aligned} \tag{39}$$

This problem can be solved numerically using appropriate software e.g., [28]. Finally, the optimal solution *ψopt*(*t*) can be applied as flux reference signal *ψ*∗(*t*). This solution could be implemented only if the reference trajectory is known in advance. A less restrictve assumption is discussed next.

In order to give more information about the solution of the boundary value problem (39), an example is presented. The parameters of the IM are taken from Table 1 and a reference torque *T*∗ *el*(*t*) is chosen as shown in Figure 2a. The MATLAB function bvp4c was used to solve the boundary value problem. The resulting optimal flux magnitude *ψopt*(*t*) is depicted in Figure 2b,c; it contains the corresponding current *iψ*(*t*) and also the component *iτ*(*t*) which is determined from the optimal solution by (33) with *ψ*(*t*) replaced by *ψopt*(*t*). The resulting value of the performance criterion (32) is given as *E* = 336.04.

**Figure 2.** (**a**) Reference torque; (**b**) optimal flux magnitude; (**c**) optimal currents.


**Table 1.** Nominal motor data and controller parameters.

#### **5. Final Controller Modifications**

In the previous section two different scenarios concerning the torque reference signal *T*∗ *el* were investigated. The general case—being the most interesting one—needs the solution of a nonlinear boundary value problem. In this case the torque reference signal must be known in advance which drastically limits the applicability. Therefore, a suboptimal but easily implementable method for the selection of the flux reference signal *ψ*∗ is preferred. A very simple but nevertheless reasonable possibility can be derived by a comparison of the optimization problems (29) and (34). In addition to the boundary conditions, the two problems differ mainly in the fact that the *dynamic* flux model is taken into account in (34). If replaced by the *static* flux model (23) a similar functional as in (29) is obtained. Based on the result (31) for problem (29) the following choice

$$\psi^\*(t) = \mathcal{g}^{-1}\left(L\_r \left| T\_{el}^\*(t) \right| \right)$$

is introduced for the general case. In [14] the special case of step changes in the reference torque (or step changes filtered by a stable linear second order system) are considered. The optimization problem for the flux magnitude in order to minimize the power losses is formulated. Finally, the optimization problem is simplified in the way that only the constant power losses in steady state are taken into account. For proposed solution it should be emphasized that *ψ*∗(*t*) is applied for general torque trajectories. Obviously this selection of the flux reference is optimal in the case of piecewise constant reference torques but it requires certain limitations in order to meet the requirements in the general case. Therefore, it is proposed to use

$$\psi^\*(t) = \text{sat}\left(\mathcal{g}^{-1}\left(L\_r\left|T\_{cl}^\*(t)\right|\right), \psi\_{\text{min}}^\* \psi\_{\text{max}}^\*\right),\tag{40}$$

where the saturation function is defined as

$$y = \text{sat}\left(\mathbf{x}, \mathbf{x}\_{\text{min}}, \mathbf{x}\_{\text{max}}\right) := \begin{cases} \mathbf{x}\_{\text{min}} & \text{if } \mathbf{x} < \mathbf{x}\_{\text{min}} \\ \mathbf{x} & \text{if } \mathbf{x}\_{\text{min}} \le \mathbf{x} \le \mathbf{x}\_{\text{max}} \\ \mathbf{x}\_{\text{max}} & \text{if } \mathbf{x} > \mathbf{x}\_{\text{max}} \end{cases} \dots$$

The lower bound *ψ*∗ min > 0 has to be introduced in order to avoid the problems at zero crossings of *T*∗ *el*(*t*) as mentioned in Section 3. Of course it is not possible to give an exact bound on the "suboptimality" of the proposed flux reference strategy for an arbitrary torque reference signal which is not given in advance. However, the choice of the lower bound *ψ*∗ min offers a simple way to adjust the concept to the expected reference torque trajectories. Especially in the case of slowly varying reference torques the efficiency improvement can be significantly compared to the constant flux reference case as can be seen in Figure 2. Upper bound *ψ*∗ max should prohibit unduly high magnetization or can be used in the field weakening range being defined as a function of the rotor speed *ωm*. The introduction of the lower bound *ψ*∗ *min* for the flux reference signal is helpful to meet different drive requirements. In cases where the dynamics of the torque response is of main interest an operation with constant flux magnitude may be preferable which can be obtained by setting *ψ*∗ *min* = *ψ*<sup>∗</sup> *max*. On the other hand, if efficiency of the drive is emphasized a small value of *ψ*∗ *min* will be a suitable choice. Of course, in cases where the required torque *T*∗ *el*(*t*) is known in advance (e.g., in robotic applications) the optimal flux reference *ψ*∗(*t*) resulting from the boundary value problem (39) can be applied directly, i.e., the first equation of the controller (41) can be omitted.

The results obtained by this simple approach are shown in Figure 3. The same torque reference signal *T*∗ *el*(*t*) as used in the previous example leads to a flux reference signal *ψ*∗(*t*) nearly coinciding with the optimal flux *ψopt*(*t*) most of the time. Here the bounds were chosen somewhat arbitrarily as *ψ*∗ min = 0.35 and *ψ*<sup>∗</sup> max = 1.4. Interestingly, the performance criterion only marginally deteriorates compared to the optimal solution. In this case the value is given as *E* = 337.6. Considering the torque reference signal in Figure 2a and the deviation from the optimal flux in Figure 3c it follows that high dynamic changes in the reference torque require an increasing value of the lower bound *ψ*∗ min in order to get a good approximation of the optimal flux magnitude.

**Figure 3.** (**a**) Optimal flux magnitude; (**b**) flux reference resulting from (40); (**c**) deviation from optimal flux.

The inclusion of the nonlinear flux model (25) into the control concept requires a minor modification of (13). The feedforward part *ψ*∗/*M* in this scheme must be replaced by the inverse magnetization function in order to achieve the correct flux magnitude.

Thus, the torque controller is proposed in its final form as follows: given the reference torque *T*∗ *el*(*t*) as external input, the estimated torque *T*ˆ *el*(*t*) and the flux coordinates *ψ*(*t*), *ϕ*(*t*) (or their estimated values *ψ*ˆ(*t*), *ϕ*ˆ(*t*) from an observer/estimator) the controller outputs *iγ*(*t*), *iδ*(*t*) (i.e., the reference signals for the current controller) are determined by the relations

$$\begin{aligned} \psi^\*(t) &= \text{sat}\left(g^{-1}\left(L\_{\tau}\left|T^\*\_{cl}(t)\right|\right), \psi^\*\_{\text{min}}, \psi^\*\_{\text{max}}\right) \\ i\_{\psi}(t) &= f^{-1}\left(\psi^\*(t)\right) + \frac{k\_{\psi}}{M}\left(\psi^\*(t) - \hat{\psi}(t)\right) \\ i\_{\tau}(t) &= \left(\frac{L\_{\tau}}{M}\frac{T^\*\_{cl}(t)}{\psi^\*(t)^2} + \frac{k\_{p}L\_{\tau}}{R\_{\tau}M}\left(T^\*\_{cl}(t) - \hat{T}\_{cl}(t)\right)\right)\psi(t) \\ i\_{\gamma}(t) &= i\_{\psi}(t)\cos\phi(t) - i\_{\tau}(t)\sin\phi(t) \\ i\_{\delta}(t) &= i\_{\psi}(t)\sin\phi(t) + i\_{\tau}(t)\cos\phi(t). \end{aligned} \tag{41}$$

**Remark 1.** *If the nonlinear flux model (25) is used in combination with the modified equation for iψ from (41) the first equation of the closed loop model (17) changes to*

$$\dot{\Psi} = -\tau\_r^{-1}(Mf^{-1}(\psi) + k\_\Psi \psi) + \tau\_r^{-1}(Mf^{-1}(\psi^\*) + k\_\Psi \psi^\*). \tag{42}$$

*Since the inequality M f* <sup>−</sup>1(*ψ*) ≥ *<sup>ψ</sup> holds for all <sup>ψ</sup>* ≥ <sup>0</sup> *it follows immediately that <sup>ψ</sup>* = <sup>0</sup> *is an asymptotically stable equilibrium of the unforced (i.e., for ψ*∗(*t*) = 0*) system (42). Additionally, using the Lyapunov functon V* = <sup>1</sup> <sup>2</sup>*ψ*<sup>2</sup> *it can be shown that (42) has the so-called input to state stability property which means that <sup>ψ</sup>*(*t*) *remains bounded whenever ψ*(0) *and ψ*∗(*t*) *are bounded [29]. Therefore, all results from the stability analysis of Section 3.2 remain valid.*

#### *5.1. Guidelines for Selection of Controller Parameters*

The final control law (41) contains two parameters *k<sup>ψ</sup>* and *kp* which must have positive values in order to achieve stability of the closed loop system. Considering the dynamic capabilities of the controlled induction machine it would be desirable to select large values for *k<sup>ψ</sup>* and *kp* but of course this choice is limited by the constraint posed on the stator current vector. Therefore, some guidelines should be given next. Starting from the given maximum possible flux magnitude *ψ*max we can determine the associated magnetization current amplitude *iψ*,*<sup>m</sup>* via the measured magnetization curve, i.e., *iψ*,*<sup>m</sup>* = *f* <sup>−</sup>1(*ψ*max). In order to achieve a desired dynamic flux control loop we may put the following constraint on the magnetization current

$$i\_{\Psi} \le i\_{\Psi, \text{max}} := 2i\_{\Psi \mathcal{M}}.$$

From (41) we have

$$i\_{\Psi} = f^{-1}(\psi^\*) + \frac{k\_{\Psi}}{M}(\psi^\* - \hat{\psi})^2$$

which immediately leads to the conservative constraint

$$f^{-1}(\psi\_{\text{max}}) + \frac{k\_{\psi}}{M} \psi\_{\text{max}} \le i\_{\psi, \text{max}}.$$

From this inequality we can derive a bound for the parameter *kψ* as

$$k\_{\psi} \le \frac{M}{\psi\_{\text{max}}} (i\_{\psi,\text{max}} - f^{-1}(\psi\_{\text{max}})) = \frac{M}{\psi\_{\text{max}}} i\_{\psi,\text{tr}}.\tag{43}$$

Given the maximum possible stator current amplitude *i*max we can determine the maximum amplitude *iτ*,max of the torque producing current by

$$i\_{\tau,\text{max}} = \sqrt{i\_{\text{max}}^2 - i\_{\psi,\text{max}}^2}.$$

Now considering (15) it follows

$$\left|\mathbf{i}\_{\mathbf{r}}\right| = \left\|\mathbf{i}\_{\mathbf{r}}\right\|\_{2} = \frac{\mathbf{r}\_{\mathbf{r}}}{M} \left(\mathbf{s}\_{3}\mathbf{r}\_{r}^{-1} + k\_{p}\vec{T}\_{el}\right)\boldsymbol{\Psi}^{\dagger}$$

from which the inequality

$$(1 + k\_p \tau\_r \mathcal{T}\_{cl}) \le i\_{\tau, \text{max}} \frac{M}{\psi\_{\text{max}}}$$

can be derived. Under the assumption

$$|\vec{T}\_{cl}| \le 2T\_{cl, \text{max}}\rho$$

where *Tel*,max denotes the maximum possible torque, we get the desired bound

$$k\_P \le \frac{1}{2\tau\_r T\_{cl,\text{max}}} \left( i\_{\text{\textpi,max}} \frac{M}{\psi\_{\text{max}}} - 1 \right). \tag{44}$$

It should be emphasized that (43) and (44) are conservative bounds; therefore, the applied parameters may be increased to some extent.

#### *5.2. Speed Control*

The proposed torque control concept can be easily extended to speed control of an IM. In this case the mathematical model of the plant has to be augmented with the differential equation for the mechanical part of the system

$$J\dot{\omega}\_{\text{fl}} = -c\omega\_{\text{fl}} + T\_{\text{el}} - T\_{\text{l}}.\tag{45}$$

The inertia of the drive is denoted by *J*, *c* is a constant for viscous friction, and *Tl* is an unknown load torque. It is assumed that the mechanical speed *ω<sup>m</sup>* can be measured. Now the proposed torque controller (41) can be simply extended by a speed controller whose output serves as reference torque input *T*∗ *el* for the torque controller. A standard PI-controller may be used for this purpose. The resulting structure of the controller including the flux and torque observer/estimator is shown in Figure 4. The reference signal for the mechanical speed is denoted by *ω*∗ *m*.

**Figure 4.** Controller structure for speed control.

For high performance servo applications it is preferable to choose *ψ*∗ below the saturation point to ensure the required dynamic response. On the other hand, for slowly changing or constant speed references *ω*∗ the flux command *ψ*∗ could be optimized as proposed.

If we dispense with the maximum torque per ampere feature a simplified version of the speed controller can be derived from (41) in the following way: First (3) of the plant model is replaced by (45), next a suitable *constant* flux reference value *ψ*∗ is chosen and finally the current component *iτ* in (41) has to be modified such that it depends on the speed error. In this way the speed controller is given as:

$$\begin{aligned} i\_{\Psi}(t) &= f^{-1}(\psi^\*) + \frac{k\_{\Psi}}{M} \left( \psi^\* - \hat{\Psi}(t) \right) \\ i\_{\tau}(t) &= \left( \frac{cL\_r}{M} \frac{\omega\_m^\*(t)}{\psi^{\*2}} + \frac{k\_p L\_r}{R\_r M} \left( \omega\_m^\*(t) - \omega\_m(t) \right) \right) \hat{\psi}(t) \\ i\_{\gamma}(t) &= i\_{\Psi}(t) \cos \hat{\phi}(t) - i\_{\tau}(t) \sin \hat{\phi}(t) \\ i\_{\delta}(t) &= i\_{\Psi}(t) \sin \phi(t) + i\_{\tau}(t) \cos \phi(t) .\end{aligned} \tag{46}$$

The closed loop system composed of (2), (45), and (46) is described by the relation (42) and

$$\begin{split} \dot{\boldsymbol{\psi}} &= \boldsymbol{c} \boldsymbol{R}\_{r} \frac{\boldsymbol{\omega}\_{\text{m}}^{\*}}{\boldsymbol{\Psi}^{\*2}} + \boldsymbol{k}\_{p} \left( \boldsymbol{\omega}\_{\text{m}}^{\*} - \boldsymbol{\omega}\_{\text{m}} \right) \\ \dot{\boldsymbol{\omega}}\_{\text{m}} &= -\frac{1}{\boldsymbol{I}} \left( \boldsymbol{c} + \frac{k\_{p}}{R\_{r}} \boldsymbol{\Psi}^{2} \right) \boldsymbol{\omega}\_{\text{m}} \\ &+ \frac{1}{\boldsymbol{I}} \left( \boldsymbol{c} \frac{\boldsymbol{\Psi}^{2}}{\boldsymbol{\Psi}^{\*2}} + \frac{k\_{p}}{R\_{r}} \boldsymbol{\Psi}^{2} \right) \boldsymbol{\omega}\_{\text{m}}^{\*} - \frac{1}{\boldsymbol{I}} T\_{\text{I}} \end{split}$$

#### **6. Control Implementation and Experimental Results**

In the following, the practical implementation of the controller combined with the appropriate inner current control along with a rotor flux linkage estimator is considered. The current control that provides reference tracking, disturbance suppression, and voltage drop compensation of the Voltage Source Inverter (VSI) was realized in the rotor reference frame. The current error

$$\vec{\mathbf{i}}\_{\gamma,\delta} = \mathbf{i}\_{\gamma,\delta}^\* - \mathbf{i}\_{\gamma,\delta}$$

was calculated from the reference currents obtained from the torque controller and the machine stator currents, that were directly measured. Based on an internal model principle a PI controller, that provides reference tracking and compensation of the back-EMF, was used in the following form:

$$\mathbf{u}\_{\gamma,\delta} = k\_{p i} \mathbf{\tilde{i}}\_{\gamma,\delta} + k\_{p i} k\_{i i} \int\_{0}^{t} \left( \frac{\mathbf{J} \omega\_{m}}{k\_{i i}} + \mathbf{I} \right) \mathbf{\tilde{i}}\_{\gamma,\delta} d\tau \tag{47}$$

The constants *kpi*, *kii* > 0 are free design parameters. In this application *kii* = *τ*−<sup>1</sup> *<sup>σ</sup>* was set. The implementation of the proposed control law clearly requires a rotor flux linkage estimator or observer. In experiments, a simplified magnetically nonlinear estimator was used based on the rotor model (25) in the following form

$$\begin{aligned} \dot{i}\_{\Psi} &= i\_{\gamma} \cos \hat{\phi} + i\_{\delta} \sin \hat{\phi} \\ \dot{i}\_{\tau} &= -i\_{\gamma} \sin \phi + i\_{\delta} \cos \phi \\ \dot{\Psi} &= \tau\_{r}^{-1} M \left( -f^{-1}(\dot{\Psi}) + i\_{\Psi} \right) \\ \dot{\hat{\rho}} &= \tau\_{r}^{-1} M \frac{\dot{i}\_{\tau}}{\dot{\Psi}} \\ \dot{\Upsilon}\_{cl} &= -\tau\_{f}^{-1} \Upsilon\_{cl} + \tau\_{f}^{-1} \frac{M}{L\_{r}} \Psi i\_{\tau} . \end{aligned} \tag{48}$$

The magnetization curve *ψ*ˆ = *f*(*iψ*) was implemented as lookup table.

The experimental setup consisted of the tested IM coupled with a separately controlled DC motor serving as variable load, a dSPACE PPC Controller board, a host PC with installed development environment, an incremental encoder with 5000 pulses per revolution, current sensors, voltage inverter, and a torque transducer. During tests, data acquisition, transformations, the flux, and torque estimator and the control algorithm were executed on the PowerPC with a sampling time of 250 μs, while a slave DSP was used for the vector modulation executed at 4 kHz. The program codes for the PowerPC and for the slave DSP were developed with the Real Time Interface and Simulink. Characteristic motor data and controller parameters for torque control are given in Table 1.

The first experiment presented in Figure 5 shows the machine response for a torque control with a piecewise constant torque reference and with adjusted rotor flux linkage magnitude as proposed in (31). The diagram (a) shows the estimated torque, the reference torque, and the filtered measured torque. For comparison on diagram (b), the same experiment is repeated using just the basic FOC (constant flux command, no compensation of magnetic nonlinearity, constant parameter flux estimator, current controllers in d-q reference frame, PI speed controller). Since no compensation of nonlinearity was used, some mismatch between observed, measured, and reference torque can be observed at higher torque values. Diagrams (c) and (d) show measured and reference currents while on the fifth diagram (e) the rotor flux linkage vector trajectory is shown. The last diagram (f) shows the reference and adjusted magnitude of the rotor flux linkage vector considering the nonlinear magnetizing characteristic.

Considering our previous proposition of the basic controller in (4), the advantages of the proposed extended control become evident if the reference changes sign with finite slope. This is illustrated for the reference profile that includes torque reversal in Figure 6. The proposed controller enables precise torque tracking also at zero crossings (a). On diagram (b) the comparison with FOC shows some deviation in the measured torque at higher reference values. Slightly better performance can be observed for the proposed control mostly due to the fact that just basic FOC was implemented. As in the previous experiment reference and measured currents (c,d), rotor flux linkage vector trajectory (e) and adjusted magnitude of the rotor flux linkage vector (f) are shown for the proposed control. Diagrams (g,h) show the stator current vector norms and corresponding integrals for a constant and optimized rotor flux. The current norm and the quantity E were obtained numerically from the measurements. For a constant rotor flux obtained with the proposed control and FOC the difference in E is a consequence of slightly bigger magnetizing current that was applied for a proposed control (curves 2 and 3). As expected, E is minimized for optimized rotor flux (curve 1). Similar performance could be achieved also with FOC for optimized rotor flux.

**Figure 5.** Torque control with adjusted rotor flux magnitude command and piecewise constant torque reference: torque response for proposed method (**a**), torque response for FOC (**b**), stator current components gamma and delta (**c**,**d**), rotor flux vector trajectory (**e**), rotor flux vector magnitude (**f**).

**Figure 6.** Torque control with optimized rotor flux magnitude command for the torque reversal including zero crossing with a finite slope (**a**–**f**) and comparison of stator current norm and corresponding integrals (**g**,**h**) for the proposed controller with optimized rotor flux (1) and with constant rotor flux (2) as well as for field oriented control (FOC) with constant rotor flux (3).

In Figure 7 the performance comparison for low speed reference is shown. No significant differences can be observed between proposed control and FOC. In the experiment shown in Figure 8 tracking performance for a harmonic speed reference with simultaneous constant load torque *Tl* = −8 Nm starting at *t* = 1.9 s is given. Actual and reference speed are shown in diagram (a), measured torque is given in diagram (b), while in diagrams (c) and (d) reference and measured currents *i<sup>γ</sup>* and *i<sup>δ</sup>* are shown. In Figure 9, the same variables as in the previous experiment are shown for a speed reversal command.

**Figure 7.** Performance comparison at low speed.

**Figure 8.** Speed response for harmonic speed reference and a constant load torque starting at *t* = 1.9 *s*: reference and measured speed (**a**), measured torque (**b**), reference and measured stator current vector components (**c**,**d**).

**Figure 9.** Control performance for the speed reversal command: speed response (**a**), measured torques (**b**), and reference and measured stator current vector components (**c**,**d**).

#### **7. Conclusions**

In this paper we proposed an extended controller for an induction machine based on the 3-D non-holonomic integrator that assures tracking of a general torque reference. This was achieved by extending a basic non-linear state controller with a separate control of the rotor flux linkage magnitude. The modified control results in a structure that is similar to FOC, although a completely different approach was used. It should be pointed out that brief additional results obtained with classical FOC are given exclusively to show that similar overall performance could be obtained with the proposed approach and are not meant for detail performance evaluation or comparison. The additional control loop in proposed control scheme provides extra flexibility to meet performance requirements posed on the drive. As occasion demands the rotor flux linkage magnitude can be optimized to assure energy efficiency or can be kept constant to provide high dynamic tracking performance in servo applications. It follows from our experience that both requirements are hard to satisfy simultaneously without prior information about the class of reference signals. Simultaneous online optimization of the flux magnitude for completely arbitrary torque reference is computationally still too demanding for existing hardware. With the proposed controller, singularity problems of the original controller (4) are completely removed. Since the nonlinear state controller requires a rotor flux linkage estimate, further improvement could be achieved by implementing well known estimation techniques in the proposed control.

**Author Contributions:** Conceptualization, B.G. and A.H.; writing—original draft preparation, B.G. and A.H.; experimental validation, G.Š.; writing—review and editing, G.Š. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the Slovenian Research Agency, project no. P2-0115 and J2-1742. **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 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/).

## *Article* **Water and Energy Efficiency Improvement of Steel Wire Manufacturing by Circuit Modelling and Optimisation**

#### **Muriel Iten 1,\*, Miguel Oliveira 2, Diogo Costa <sup>3</sup> and Jochen Michels <sup>4</sup>**


Received: 29 November 2018; Accepted: 3 January 2019; Published: 11 January 2019

**Abstract:** Industrial water circuits (IWC) are frequently neglected as they are auxiliary circuits of industrial processes, leading to a missing awareness of their energy- and water-saving potential. Industrial sectors such as steel, chemicals, paper and food processing are notable in their water-related energy requirements. Improvement of energy efficiency in industrial processes saves resources and reduces manufacturing costs. The paper presents a cooling IWC of a steel wire processing plant in which steel billets are transformed into wire. The circuit was built in object-oriented language in OpenModelica and validated with real plant data. Several improvement measures have been identified and an optimisation methodology has been proposed. A techno-economic analysis has been carried out to estimate the energy savings and payback time for the proposed improvement measures. The suggested measures allow energy savings up to 29% in less than 3 years' payback time and water consumption savings of approximately 7.5%.

**Keywords:** energy efficiency; industry; water circuits; OpenModelica; optimisation

#### **1. Introduction**

Industrial Water Circuits (IWCs) encompass all the water streams of a factory unit, including the most common uses of water industrially such as water cooling and water treatment. Such processes are responsible for significant consumption of water as well as energy. The industrial sector is responsible for 20% of the total water consumption in the world [1], while in terms of energy it represents 30% of the global final end-use of energy [2]. The European industry in particular is responsible for 40% of the total water abstractions [3] and 25% of the final end-use of energy in the European Union is attributed to this sector [4]. Electric motors, in turn, represent 65% of the total electric energy consumption in European industry [5] and pumping systems in particular account for 21% of the aforementioned share of electricity consumption [6]. These figures show that it is crucial to improve the energy efficiency of IWCs.

The 2012 Energy Efficiency Directive [7] and implementing directive for the Ecodesign requirements for water pumps [8] focus on reducing energy consumption and consequently, water consumption. These directives emerged in the scope of the climate change and energy targets of the Europe 2020 Strategy [9]. Hence, increasing attention have been given to alternative strategies for the reduction of energy and water consumptions.

The European project WaterWatt [10] aims to develop new tools and guidelines to reduce energy consumption in IWCs. IWCs are generally treated as secondary streams in the production process

of a plant and therefore neglected. Furthermore, the associated energy consumption is frequently despised, as the circuits only represent a minor part of the total electricity consumption. In this alignment, the present paper presents an object-oriented model, developed in the open source OpenModelica Software in order to analyse and improve the energy consumption of an IWC case study part of the WaterWatt project [11].

The optimisation problems approached in this work are based on the application of improvement measures to reduce electric energy consumption and water consumption in IWC. The measures for energy efficiency improvement in IWC, and similar systems such water networks and water distribution systems, have been proposed by several authors. Cabrera et al. [11] proposed two sets of strategies for energy efficiency improvement, an operational-type set and structural-type set. In operational group, it was proposed: the operation of the pumping systems at the best efficient point (BEP), the improvement of the system regulation to avoid surplus energy, the minimization of leaks and the minimization of friction losses. For the structural group, it were proposed: the use of more efficient pumps, the decoupling of energy sectors and the improvement of designs and layouts.

As part of the WaterWatt project and this work specifically, the optimization based on energy efficiency improvement of the IWC will be performed taking into account the knowledge about the proposed measures for this type of systems and the knowledge about the achievement of energy savings by the adjustment of the operation of electric motors. The application of variable speed drives (VSD) in electric motors is proposed to attain considerable reductions in both water and energy consumptions, which in the case of the pumping systems it has been revealing as an excellent alternative to currently applied improvement measures.

In this paper, different improvement measures are proposed, and energy savings and payback associated to these are estimated. The case study corresponds to water cooling circuit of a steel wire processing plant in Germany in which steel billets are transformed into wire. In general, iron and steel industry present an electrical energy consumption of 2327 ktoe/year, which represents about 12% of the total energy consumption of the industrial sector namely in Germany [12]. Moreover, the production of basic metals in the manufacturing industry in Germany has a total water consumption of about 581 dam3, in which 85% corresponds to water cooling [13].

#### **2. Model Assembling**

#### *2.1. Circuit Description*

The cooling circuit under analysis includes a hot roiling mil, four cooling towers, three pump groups, four sand filters, one oil separator and a cyclone separator (Figure 1). The pump group 1 is constituted by two pumps that transport cooled water from a tank (Cold Water Tank) to the rolling mill. Then, the water stream flows down the channel to an oil separator and cyclone (Cyclone Tank) for the pre-treatment of water. The pre-treated water is distributed to two sand filters by the action of pump group 2, constituted by two pumps in operation. From the sand filters, the water stream flows back to the cooling tower tank. The designated pump group 3 is part of a secondary circuit, constituted by two pumps. Pump group 3 corresponds to the same type as pump groups 1 and 2 to reduce maintenance and logistic effort. This separate water stream is transported from the cooling tower tank to the four cooling towers and the outlet water streams from the cooling tower are then, transported back to the tank. This secondary circuit is part of a recent modernisation drive in the factory. Such modernization aims to improve energy efficiency. It was installed to allow the operation of the cooling towers only when required (i.e., when the temperature at the cold water tank increases above the operational target-value by 0.5 ◦C) and as soon as temperature drops the circuit stops. For instance, this may constitute a constraint for the process optimization based on energy efficiency improvement, in which the temperature of the cold water tank must not be exceeded by 0.5 ◦C relatively to its target-value. Otherwise, this operational requirement is considered to be surpassed. The separation of the cooling towers from the main circuit also allows a closer control on the cooling process as the water demand in the circuit is variable. This is

because the billets are produced in batches and due to the existence of variabilities in the material and the rolling mill operational conditions, such as rolling speed and cooling speed.

The inlet water stream in the rolling mill is used for the main purpose of direct and indirect cooling. Nonetheless, it used in the remaining processes of hot forming, such as descaling and scale transport. For that reason, a water treatment section is installed after the hot forming process. This section is constituted by the filtration units, constituted each one, by two sand filters. Moreover, it is preceded by an oil separator and a cyclone, constituting a pre-treatment section.

**Figure 1.** Schematics of the case study cooling circuit.

Table 1 summarizes the specifications for the design and assembling of the model circuit (described Section 2.2).


**Table 1.** Specifications for the assembling of the cooling circuit.


**Table 1.** *Cont.*

The heat input at the hot rolling mill (*qinput*,*RM*) is determined by:

$$
\rho\_{RMinput} = Q\_1 \times \rho\_w \times c\_{pw} \times \left( T\_{RMinlet} - T\_{RMoutlet} \right) \tag{1}
$$

where *Q*<sup>1</sup> and *Q*<sup>2</sup> corresponds to the water flowrate of the main and secondary circuit respectively (Table 2). The density of the water (*ρw*) corresponds to 999 kg/m3 and the specific heat (*cpw*) to 4.186 kJ/kg◦C. The *Tinlet\_RM* is determined by:

$$Q\_1 \times \rho\_W \times C\_{P,w} \times (T\_{in,1} - T\_{out,1}) = Q\_2 \times \rho\_W \times C\_{P,w} \times (T\_{in,2} - T\_{out,2}) \tag{2}$$

whereas:

$$T\_{out,1} = T\_{RMintet} \tag{3}$$

These two additional inputs of the circuit are listed in Table 2.

**Table 2.** Additional inputs in the IWC.


The measurements at the plant included the water temperature, water flowrate and pressure drop across pumps and sand filters (Table 3). In order to measure the flowrate, ultrasonic signals have been used, employing the transit time difference principle. The ultrasonic signals are emitted by a transducer installed on the pipe and received by a second transducer. The transit time difference is measured and allows the flowmeter to determine the average flow velocity along the propagation path of the ultrasonic signals. A flow profile correction is then performed to obtain the area averaged flow velocity, which is proportional to the volumetric flow rate. Two integrated microprocessors control the entire measuring process, removing disturbance signals, and checking each received ultrasonic wave for its validity which reduces noise. Moreover, the temperature of the water circulating within the pipes was measured by coupling two resistance thermometer (RTD) temperature probe (applicable to the supply and return pipes simultaneously) whose outputs are integrated with the flow measurements at the meter. The pressure drops across pumps and filters have been measured by the coupling of manometers before and after each component. Also, the power of pumps have been recorded during

representative period of time. This has been performed by a power analyser connected to a data logger, taking three phase power measurements by coupling clamp meters into the phases. Table 1 summarizes the specifications for the design and assembling of the model circuit (described Section 2.2).


**Table 3.** Specifications for each type of measurement.

#### *2.2. Design and Assembling of IWC Model*

The IWC was built in the object-oriented language in OpenModelica, used for the modelling of such physical systems and supports most of the equipment (Figure 2).

**Figure 2.** Diagram of the circuit of a steel wire processing plant in OpenModelica (Legend: (1) Pump Group 1, (2) Hot Rolling Mills, (3) Cyclone Tank, (4) Pump Groups 2, (5) Pipes from pumps to filters, (6) Sand Filters, (7) Pipes from filters to tank, (8) Pump Group 3, (9) Pipes from pumps to cooling towers, (10) Cooling Towers, (11) Cold Water Tank, (12) Reference System).

Currently, a number of free and commercial component libraries in different domains are available, including electrical, mechanical, thermo-fluid and physical-chemical. The current case study was built using the ThermoPower library [14] that has been developed with the goal of providing a solid foundation for dynamic thermal and power modelling. It includes a list of components, such as pumps, tanks, cooling towers, pipes and filters, that can be "drag and drop" for the modelling assembling. Each component corresponds to singular model to be setup accordingly with the required parameters and therefore, replicability is assured.

The selected components and respective models correspond to the main components of the circuit namely: tank model, pump model, pipe model, filter model, cooling tower model and rolling mill as further described. Each model was accordingly setup using the measured data collected in the case study circuit and presented in the beginning of this section.

#### 2.2.1. Tank Model

The model represents a free-surface water tank corresponding to the vertical axis in the region where the fluid level is supposed to lie within. The tank geometry is described the input parameters: *V*<sup>0</sup> (volume of the fluid when the level is at reference zero height), *A* (cross section of the free surface) and *y*<sup>0</sup> (height of "zero-level"). The model also requires the input parameters *ymin* and *ymax*, which should be appropriate to avoid underflow or overflow of the tank. The pressure at the inlet and outlet ports accounts for the ambient pressure and the static head due to the level, while the pressure at the inlet port corresponds to the ambient pressure.

Concerning the mass and energy balances the model assumes there is no heat transfer except through the inlet and outlet flows. It is to note that, due to the inexistence of the performance of heat transfer directly in the modelling of the cold water tank in the present IWC, the sand filters outlet water streams are cooled through the contact with the water stream from the secondary circuit at a lower temperature.

#### 2.2.2. Pump Model

The model represents a centrifugal pump, or rather a group of centrifugal pumps in parallel. The input parameters are related with the pump curves therefore defined in the model for each single pump or pump group. The component characteristics are given for nominal conditions of flowrate and rotational speed. These characteristics may be then adapted according to similarity equations, namely the affinity laws for the pumps to different flow conditions (Equation (7)).

There are two input functions: *flowCharacteristic* that corresponds to the relationship of the volumetric flowrate and the head, while the *efficiencyCharacteristic*, corresponds to the relationship between the volumetric flowrate and the hydraulic efficiency. One of the output parameter is the electric power of the pump (Pump). This also accounts for the mechanical efficiency of the pump motor (*motorEfficiency*) also set as input.

In addition to the inlet and outlet ports of the water stream, the pump models also include inlet ports for the motor speed (*in\_n*), the electricity tariff (*tax*) and the number of pumps in parallel (*in\_NP*).

#### 2.2.3. Pipe Model

The model represents the piping system of the circuit, either as a single tube or by N parallel tubes, representing the water flow streams of the factory unit. It is built with certain assumptions, such as: one-phase fluid state, uniform velocity across the section, constant turbulent friction, longitudinal heat diffusion is neglected and uniform pressure distribution in the energy balance equation.

The finite volume method is used to discretize the mass, momentum and energy balance equation, taking as the state variables one pressure, one flowrate and N-1 specific enthalpies. The dynamic momentum is not accounted by default.

The input parameters for the pipe model are: the pressure losses on the pipes due to the friction phenomena (*dpnom*) and the static head (*h*) as well as the mass flowrate (*wnom*). In addition to the inlet and outlet ports of the water stream, the pipe model also contains a wall-type inlet (*wall*), which can be specified for a certain heat input. Taking into account that the circuit pipes are all installed at the same level, this is, it does not exist pressure losses due to elevation, the total pressure drop of the sections of the piping system are determined considering the head losses due the friction of the pipes, according to:

$$
\Delta P\_f = \rho\_w g f \frac{L}{D} \frac{\mu^2}{2g} \tag{4}
$$

The friction factor (*f*) was determined using the Moody pipe friction chart, according to a transitional flow regime (*Re* = 24,300). Moreover, for the total head losses the minor losses are also considered. The minor losses consist of the local pressure drop (Equation (5)) due to sudden or gradual expansion or contraction, the presence of bends, elbows, tees and valves (open or partially closed).

$$
\Delta P\_L = \rho\_w g f K\_L \frac{\mu^2}{2g} \tag{5}
$$

The *KL* is determined according to tabulated values for each type of local pressure drop [15].

#### 2.2.4. Filter Model

The model is set to induce and determine the pressure loss in the circuit. The pressure drop across the inlet and outlet connectors of the components is proportional to the square velocity of the fluid, defined according to a turbulent friction model. The input parameter corresponds to the pressure drop across de filter (*dpnom*).

#### 2.2.5. Cooling Tower Model

The model represents an open cooling tower with fans and it is filled with metallic packaging. The humid ambient air is used to cool the incoming water stream. The fluid is modelled using the *IF-97 water-steam model* as an ideal mixture of dry air and steam.

The hold-up of water in the packaging is modelled assuming a linear relationship between the hold-up in each volume and the corresponding outgoing flow. The finite volume method is used to discretize the corresponding 1-D counter-current equations. The Merkel's equation [16] is used for the mass and heat transfer phenomena from the hot water to the humid air (Equation (6)):

$$\frac{k\text{ S}}{Q\_{Mw}} = \text{C}\_{pw} \int\_{T\_{w,air}}^{T\_{w,in}} \frac{1}{(H\_s - H\_{air})} dT \tag{6}$$

Note that *k* designates the overall heat and mass transfer coefficient and *S* the surface of packing. The required input parameters correspond to: the overall mass and heat transfer (*k\_wa\_tot*), nominal water mass flowrate (*wlnom*), nominal air volume flow rate (*qanom*), normal air density (*rhoanom*), surface of packing (*S*), nominal fan rotational speed (*rpm\_nom*) and the nominal power consumption (Wnom). The fan behaviour is modelled by kinematic similarity: the air flow is proportional to the fan speed (*rpm\_nom*) and the electric power consumption (*Wnom*) is proportional to the cube of the speed, according to the affinity laws for the fans (Equations (7) and (8)).

The cooling tower model also contain an inlet port for the fan speed (*fanRpm*). In which, similarly to the pump motor rotational speed, the fan speed can be changed relatively to the nominal value (*in\_n*). The outcome parameters include the inlet (*TCTinlet*) and outlet (*TCToutlet*) ports of the water stream, there is also an outlet port for the electric power consumption of the fan (*powerConsumption*).

#### 2.2.6. Rolling Mill

The rolling mill connects two models from the ThermoPower library, namely the heat source model and the pipe model already described above. Regarding the heat source model it corresponds to an ideal tubular heat flow source, with uniform heat flux. The actual heating power (*P*) is provided as input by the power signal connector (*in\_p*). The finite volume method is used for discretisation. The output parameter corresponds to the water outlet temperature (*TRMoutlet*).

While the comprehension of the thermal phenomena intrinsic in the analysis of this type of circuits have been performed for the design and calibration assuming the more accurate values possible, it is to note that the cooling system of a rolling mill component has been simplified. This represents a thermal power input in the circuit, allowing the increase of temperature proportional to that heat input and water flow rate. Therefore, it still lacks a more complex understanding of the occurrence of phenomena at a more dynamic basis.

#### 2.2.7. System and Time Tables

The System component must be added in any circuit in order to define the system condition, namely water temperature and ambient pressure. In this study, the temperature was considered approximately 15 ◦C according to the average water temperature provided by the utility company in Germany. The Time Table component allows the user to specify a certain input value over time for a specific component. In the current case study, they are connected to the pump motor and to the fan of the cooling tower to specify the rotational speed of these equipment(rpm). Furthermore, it can also be connected to the rolling mill to specify the thermal input change with time. Additionally, a time table can be connected to specify the electricity tariff price with time.

#### **3. Model Validation**

The model results have been validated with real casa data from the case study. For the energy and water consumption the following parameters have been identified for validation: pumps power (*Pump1*, *Pump4*, *Pump8*), the inlet and outlet water temperature at the cooling towers (*TCTinlet* and *TCToutlet*) and the water outlet temperature of the rolling mills (*TRMoutlet*). The outputs of the numerical results for these parameters are presented from Figures 3–5.

**Figure 3.** *Cont*.

**Figure 3.** Numerical and measured pump power (**a**) pump group 1, (**b**) pump group 2 and (**c**) pump group 3.

Furthermore, it is also relevant to analyse the values obtained for the temperatures in different points of the circuit, namely the inlet and outlet water stream temperature at the cooling tower (Figure 4) and outlet water temperature at the rolling mill (Figure 5).

The deviation that has been observed for each pump group to comparison to the nominal value was listed in Table 3 above. The measured electric power used for the operation of the pumps was analysed and validated against plant data and summarized in Table 1. It is observed that the values obtained in the simulation are overall consistent with plant data. Moreover, the value outlet water stream temperature of the hot rolling mill and at the inlet of the cooling towers have minor deviations compared to the real values. The slightly major deviation is related with the outlet temperature at the cooling towers (Table 4). This may be related with the fact that the water circuit in the cooling towers corresponds to a separated circuit from the main circuit. Overall, the thermal phenomena are consistent with plant data and the required cooling load is given by the model.

**Figure 4.** *Cont*.

**Figure 4.** Water temperature at cooling tower (**a**) inlet: *TCTinlet* and (**b**) outlet: *TCToutlet*.

**Figure 5.** Water outlet temperature of rolling mill (*TRMoutlet*).

**Table 4.** Deviations obtained for the pumps measured power.


#### **4. Sensitivity Analysis**

A sensitivity analysis of the model was performed. It corresponds to a crucial analysis for the solving of the optimisation problems. The methodology to perform this analysis is distinct from the ones elaborated in dynamic modelling and simulation. As there is a high number of required inputs to be given to the model in order to run, a pre-selection of variables to analysis is required. The selected variables for the analysis correspond to the speeds of pump motor and the cooling tower fan which have a direct influence on power consumption and operational requirements of the IWC. Different scenarios have been defined for these variables (Tables 5 and 6) to understand their influence on the power consumption (Figure 6). The scenarios have been assembled following a methodology similar to the Monte Carlo method. For this study, is based on random sampling, in which the results for power consumption are obtained for random rotational speeds values. Thus, different values for power consumption are obtained for different scenarios. To note that the best scenarios have been selected based on the lower power consumption of all the pump groups and also without surpassing the operational requirements.


**Table 5.** Scenarios for the variation of Pump Rotational Speed.

Moreover, it is also necessary to evaluate the influence of those input parameters on simulations results related to temperature and hydraulic parameters. In this study, the influence of the change of the motors speed of pumps and fans (cooling tower) on the water temperature at the cold water tank will be analysed (Figures 7 and 8).

For the sensitivity analysis related with the change of the pump motor rotational speeds, the power consumptions of both pump groups 1 and 2 were analysed simultaneously, since both are part of the main circuit, and having the same technical characteristics. Figure 6 presents the results ten scenarios for the power consumption of pump groups 1, 2 and 3.

**Figure 6.** Influence of the rotational speed of the pump motors on the power consumption.

From the ten created scenarios, it is necessary to select the ones associated to higher reduction in power consumption of both pump group 1 and 2 and pump group 3. This analysis focuses on the favorable scenarios on the viewpoint of energy efficiency and the possibility of its application. Therefore, observing Figure 6, the three scenarios correspond to 5, 8 and 10 and to be selected to

proceed with the sensitivity analysis (Figure 7). The remaining scenarios were neglected, as they do not present significant reduction in power consumption.

**Figure 7.** Influence of rotational speed of the pump motors on the water temperature of the cold water.

**Table 6.** Scenarios for the variation of Cooling Tower Fans Rotational Speed.


The scenarios for the cooling towers (Table 6) were created considering the share of reduction in power consumption by the decrease of the rotational speed of the cooling tower fans, following the affinity laws of the fans (Equation (8)). Thus, three scenarios designated 1, 2 and 3 were created bearing in mind shares of reduction in power consumption of, respectively, 20%, 27% and 50%.

**Figure 8.** Influence of rotational speed of the cooling tower fans (nCT) on the water temperature of the cold water tank (T).

The most efficient scenarios for each case: promoting the major energy reductions in power consumption respecting the operational requirements in temperature, will be analysed furtherly.

#### **5. Identification of Improvement Measures**

The efficiency of water circuits can be reached by the dynamic adjustment of the circuit operation to the production process, through the reduction of pressure demand and pressure losses in the circuit and by improvement of energy efficiency of pumps and cooling tower. Such improvements include matching the scale of the motor service to the work demand (i.e., speed control devices such as adjustable speed drives); reducing demand for energy services (e.g., substituting a blower for compressed air or turning off steam supplied to inactive equipment) and the replacement to high efficiency motors. The improvement measures should account for the technical specifications of the process: avoiding interruptions, improvement of production reliability (maintaining required flow, pressure and temperature) and improved maintenance practices. The measures do also account for economic and environmental aspects namely cost reduction, energy savings and reduction of CO2 impact which are the primary aspect that industrials consider regarding energy efficiency.

The present model allows a prior analysis of efficiency measures to be implemented in the IWC avoiding extra costs with technical work. To note that the structural related input variables such as: length and diameter of the pipes, heat inputs in the case of non-looped circuits, despite being highly influential and part of the formulation of certain optimisation problems, they are set as design variables, meaning constant parameters of the assembling of the circuits. Therefore, the analysis of such measures, as related to energy efficiency measures at the design stage and therefore not in the scope of this work.

Furthermore, aim of this work does not focuses on the analysis of process control strategies. Noteworthy such strategies are intrinsic in the operation of all the processes in industry, including the water circuits. Thus, a more accurate analysis of these models would require the comprehension of the models related to types of valves, flowmeters, electric meters and temperature sensors. Such extensions, would possibly assist the operational improvement measures, as well as the application of control related measures. Hence, this paper analyses the application of five improvement measures to the case study circuit in order to reduce the energy and water consumption as well as guaranteeing the required operation conditions.

#### *5.1. Variable Speed Drives (VSD) in Pumps*

The coupling of VSD's allow energy savings of 20% to 25% [17]. The installation of VSD in pump motors allows the automatic flow adjustment to the needs of the process. The change of the pump speed adjusts the pump rotation frequency to the optimal efficiency point. The performance of this measure relies on the affinity laws for the pumps. The automatic flow adjustment is described by Equation (7) [18]:

$$\frac{Q\_{w,nom}}{Q\_{w,opt}} = \frac{Q\_{pump,nom}}{Q\_{pump,opt}}\tag{7}$$

For the cooling circuit under analysis, the operation of the pump motors was initially set at rotational speed of 1450 rpm for the nominal flowrate. Applying the Equation (7) for the real running flowrate, the optimised rotational speed corresponds to 1350 rpm. New simulations have been performed for the optimal point speed and savings have been estimated and presented in Table 7.

#### *5.2. Refurbishment of Pumps*

The refurbishment of pumps corresponds to the mechanical cleaning and overall of a pump to restore its initial functioning, namely its energy efficiency. Pumps without maintenance over the years generate a lower flow. The improvement measure for the pumping system of the case study circuit considered the refurbishment of all pumps. For the calculations, the authors have assumed an improvement of 10% in the hydraulic efficiency of the pump, set in IMechE [19].

#### *5.3. Change of Filters*

In water treatment process, sand filters present considerable impact on energy efficiency due to lower pressure drop induced into the water. The improvement measures for the treatment section was to. In replace the two sand filters with a pressure drop of 1.3 bar to new ones with a pressure drop of 0.5 bar. The results are presented in Table 7.

#### *5.4. Change of Electric Motors of the Pumps*

The change of electric motors in pumping systems, specifically to electric motors categorized as IE3 Premium Efficiency allow high energy savings. A higher efficiency motor is more expensive than a conventional motor however its lifespan is much longer as it operates at a lower temperature and hence heat losses are lower [20]. In the current circuit, the pump motors with an initial mechanical efficiency of 90% were exchanged to IE3 electric motors with a mechanical efficiency of 95.4%. The savings are presented in Table 7.

#### *5.5. VSD in the Cooling Tower Fans*

The installation of VSDs in the cooling tower fans allow a dynamic adjustment of the airflow respecting the required temperature values in the circuit. The power consumption associated to the operation of the fan is directly linked to the fan speed, therefore high energy savings can be achieved with adjustments decreasing the fan speed [21]. The performance of this measure relies on the affinity laws for the fans. The automatic air flow adjustment is described by Equation (7) [22]. Moreover, the change on the power of the fans is described by Equation (8) [18]:

$$\frac{Q\_{air,nom}}{Q\_{air,opt}} = \frac{\eta\_{fan,nom}}{\eta\_{fan,opt}}\tag{8}$$

$$\frac{P\_{fan,nom}}{P\_{fan,opt}} = \left(\frac{\eta\_{fan,nom}}{\eta\_{fan,nom}}\right)^3 \tag{9}$$

In the circuit under analysis the operation of the four cooling tower fans operates at an initial rotational speed of 1450 rpm. This has been changed to an optimal speed of 1350 rpm. The results are presented in the Table 7.

#### **6. Definition of Optimization Methodology**

An optimisation methodology must be developed with the aim to improve energy efficiency in this IWC. This methodology is primarily based on the assembling of several settings for each case, in which a certain input parameter is modified with the aim to achieve lower energy consumption. These are primarily related with the operational conditions of the pumps and cooling towers. The definition of this methodology is necessary for the techno-economic evaluation of the analysed improvement measures. Based on the scenarios assembled in the sensitivity analysis, the objective functions are attained by the changes on the decision variables, although that change must respect the operational requirement of the circuit which are translated by the constraints.

The main objective of the presented optimisation problems is the reduction on electric energy consumption. Although thermal energy consumption is also a relevant concern, and heat is a direct input in certain components of these models, the thermal phenomena are mostly treated as imposed variables of the model, namely as input parameters. Hence, the optimisation methodology includes the definition of objective functions, decision variables and constraints.

The secondary circuit is only switched on when the temperature of the water inside the cold water tank (*TCWT*) reaches more than 0.5 ◦C relatively to the target value (22.1 ◦C). From this, the constraint is defined as a requirement of not exceeding 0.5 ◦C of the temperature at the cold water tank. This allows to maintain the cooled water temperature close to the required in the cooling process.

The objective functions correspond to the reduction of energy consumption. The decision variables are the same as the ones selected for the sensitivity analysis: rotational speeds of pump motors (*nPG*1, *nPG*<sup>2</sup> and *nPG*3), cooling tower fans (*nCT*) and the hydraulic and mechanical efficiencies of the pumps (designated by *ηhydr* and *ηmech*, respectively) as well as the pressure drop in the filters. Table 7 summarizes the objective functions, decision variables and constraints of the optimisation problems.


**Table 7.** Deviations obtained for the pumps measured power.

The scenarios presented in the sensitivity analysis enables high energy savings by respecting the circuit operational conditions. For the solving of optimisation problems, the scenarios associated to the lower energy consumption respecting the process limit conditions, this is, its constraints, will be selected and analysed for energy efficiency improvements. The application of VDS's in pump motors is related to scenario 8 for the variation of pump motor speed. The application of VSD's in cooling tower fans are related with the scenario 2 for the variation of fan speed.

#### **7. Techno-Economic Evaluation**

A techno-economic analysis was undertaken for the implementation of specific improvement measures. The economic savings of an improvement measure, determined in a per year basis, are, in general, calculated according to Equation (9), in which *Savings* designates power savings, *top* the annual operational time, *C* the electricity cost per energy unit and *Pnom* and *Popt* the power at initial conditions and optimized conditions, respectively:

$$\text{Savings} = \text{t}\_{op} \times \text{C}\_{\text{elec}} \times (P\_{\text{man}} - P\_{\text{opt}}) \tag{10}$$

In the case of the implementation of electric motors with higher efficiency, the method to calculate economic savings by Equation (10) takes the form of Equation (11):

$$Savings = t\_{op} \times \mathbb{C} \times \left(\frac{P\_{IE1}}{\eta\_{IE1}} - \frac{P\_{IE3}}{\eta\_{IE3}}\right) \tag{11}$$

The *PIE*<sup>1</sup> and *ηIE*<sup>1</sup> correspond to the pump power and its overall efficiency at initial conditions of a standard efficiency motor *E*1. While *PIE*<sup>3</sup> and *ηIE*<sup>3</sup> correspond to the pump power and its overall efficiency of a premium efficiency motor *IE*3.

In the case of the application of VSD's in electric motors, the savings are calculated according to Equation (12), in which *Pnom* and *ηnom* designate, the power of a pump and its overall efficiency at initial conditions respectively; *Popt* and *ηopt* correspond to the power of a pump and its mechanical efficiency, respectively, by the installation of VSD's and *i* to the numeral designation of a load regime:

$$Savings = \sum\_{i} t\_{op} \times \mathbb{C} \times \left(\frac{P\_{nom}}{\eta\_{nom}} - \frac{P\_{VSD}}{\eta\_{VSD}}\right) \tag{12}$$

Since the electric motors in question, namely pump and fan motors, are operated in a continuous load regime, Equation (12) is simplified to:

$$Savings = t\_{op} \times \mathbb{C} \times \left(\frac{P\_{nom}}{\eta\_{nom}} - \frac{P\_{VSD}}{\eta\_{VSD}}\right) \tag{13}$$

The investment payback period of a given measure is calculated according to Equation (14), in which *PB* designates the payback period and *Inv* the investment cost:

$$PB = \frac{Inv}{Savings} \tag{14}$$

The implementation of an improvement measure is considered economically viable under the patterns of European industry if the payback period is below 2–3 years. The economic evaluation was proceeded by considering an average electricity price in Portugal of 0.0923 €/kWh [23] and a simple payback, not considering an inflation rate of the investment costs. The cost of the electric motors has been obtained from the manufacturer catalogue [24] as well as the cost of the variable speed drives [25]. The cost of the refurbishment of each pump was estimated to be 2 full days of working hours of a maintenance technician with a monthly salary of 800 Euros. It is assumed that the material required for repairs and exchanges of equipment are intrinsic to the company hence no cost was considered. The energy savings, a share of energy savings, investment cost and the payback values are presented in Table 8.

**Table 8.** Energy savings, investment cost and payback period of proposed measures.


Furthermore, the application of VSD's in pump motors also allows the achievement of savings in water consumption by the decrease of water flow rate (Table 9).

**Table 9.** Savings in Water Consumption.


#### **8. Discussion**

Overall all the applied improvement measures allowed energy savings. The payback period was generally lower than 3 years, with exception to measure iv) and largely overcoming that period and therefore not considered as an efficient economic measure.

In general, the application of VSD in the pumps (measure i.1, i.2, i.3) allows energy consumption reductions from 19% to 20%. This converges with the typical values of energy savings with such application [19]. Furthermore, it reveals itself as favorable measure due to its low payback time of approximately from 11 mouths to 12 months for all pump groups. The coupling of VSD in the cooling tower fans (measure v) is equally an attractive measure, with energy savings of 27% and payback period of nearly a year.

The refurbishment of pumps (measure ii) is peculiarly advantageous because despite saving less than the applying of VSD it presents a low payback time. This is due to the low investment cost of the implementation of such measure. Adding the application of VSD in pumps and its refurbishment is shown as the most propitious measures since it produces the larger reduction in energy consumption associated to an acceptable payback period.

The replacement of conventional pump motors to the IE3 premium efficiency pumps has not the same propitiousness as the aforementioned measures. This is also applied to the change of filters. However, to note that motors and filters are changed at the end of their lifecycle and savings can be expected at that time. If the replacement is joined with the coupling of VSD however it becomes a very attractive measure presenting a payback of 2.5 years for pump group 1 and 2 and 2.2 years for pump group 3. Such approach is also beneficial on a technical point of view as the application of VSD reduces of the motors lifetime and replacing to a high efficiency motor enables to overcome such limitation [17].

From an environmental perspective, the application of VSD in pumps allows to reduce fresh water consumption. The current nominal water flow rate corresponds to 600 m3/h for the pump group 3. From the simulation, it was observed that required the cooling load was achieved for a flowrate of 555 m3/h. Taking into account that average annual water consumption per inhabitant in Germany in 2016 was 45 m<sup>3</sup> [26] the water savings could meet the needs of 6601 inhabitants.

#### **9. Conclusions**

The present work presents the modelling and optimisation of an IWC, with the following steps: model assembling, modelling, simulation results, model validation, sensitivity analysis, formulation of optimisation problems, implementation of energy efficiency improvement measures and a techno-economic evaluation.

The simulation results for the presented IWC have been successfully calibrated against real measured data. The deviations corresponded to: pump group 1 (0.74%), pump group 2 (0.77%), pump group 3 (1.56%), the cooling tower inlet (2.62%) and outlet (4.75%) and the rolling mill water outlet temperature (3.29%).

Overall, the hydraulic phenomena, namely the pump power demonstrate lower deviations than the thermal phenomena. This is due the hydraulic phenomena being directly related to energy use of the pumping systems that correspond to the main power consumers of the IWC. Thus, they are also associated to the main improvement measures. While the thermal phenomena are related with temperatures requirements that are associated to the limit conditions. Furthermore, several improvement measures have been identified and their implementation allowed up to 29% share of energy savings. For this case study, indeed there is a large potential for energy savings and improvement of the efficiency in related industrial water circuits (IWC). This has been envisaged in general by the WaterWatt project [7] and demonstrated through the case studies part of the project.

The presented methodology has been successfully replicated for the other steel sector case-studies of the WaterWatt project [7] including: closed cooling circuit of an inductive furnace and closed cooling circuit of a blast furnace. Similarly, for each case-study, simulation results have been validated against real measured data, optimisation problems have been formulated and sets of scenarios have been defined for the improvement measures as well as for the techno-economic analysis.

**Author Contributions:** M.I., M.O. and D.C. conceived and designed the experiments; M.O. and D.C. performed the experiments; M.I. and J.M. analyzed the data; J.M. contributed with materials and analysis tools; M.I. and M.O. wrote the paper.

**Funding:** This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement "No 695820".

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

#### **Nomenclature**



#### **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/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Energies* Editorial Office E-mail: energies@mdpi.com www.mdpi.com/journal/energies

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18