Optimal control of electrical machines using Pontryagin’s maximum principle

Abstract

This paper presents an optimization method that deal with control techniques for electrical machines, in order to impose a displacement with a minimum of energy consumption. It’s a generic method; However, its application to electrical machines is complex for implementation. The present work provides the shape of the optimal solution. Moreover, a resolution technique related to the nature of the model is proposed, which facilitate the numerical resolution of costate problem. The method is applied to hybrid stepper motor and the obtained results are then analyzed.

Keywords

Efficiency optimal control dynamic programming Pontryagin’s maximum principle hybrid stepper motor

1. Introduction

The improvement of energetic efficiency of electrical machines is a field of study of high interest. In fact, the electrical consumption of electrical machines is about 46% of the total electrical consumption in the world. A 1% increase in energetic efficiency of electrical machines would save almost 80 million of kWh per year, nearly 14 million barrels of oil. Therefore, the energetic efficiency of machines is a crucial issue for the preservation of energetic resources, the environment and sustainable development [1,2].

Energy consumption problem of electrical machines could be formulated in three different ways. The first approach called, optimal design problem, consists in acting over the design of the machine and/or on the nature of materials in order to find the best design to reach the desired performances with the minimization of a certain criterion. In this case, dimensions and shapes evolve [3]. A second approach is to find optimal controls that ensure good performances while minimizing energy consumption subject to certain constraints. This is an optimal control problem. In this case, the main purpose is to find optimal control variables while dimensions are fixed [4]. Finally, it’s possible to simultaneously search for sizing and optimal control. This is the global optimization. In this sense, both dimensions and control are sought simultaneously. This approach leads to a more complex problem but more suitable for application [5].

Whatever the considered approach, there are usually two major difficulties. First, the problem should be correctly formulated by defining considered variables and the objective function. Then, the chosen optimization method should find the global minimum within an acceptable calculation time. Given that, the purpose is to conceive an electrical machine with an optimal control and design, the global optimization problem turns out to be the most interesting. However, multi-level approaches, which are widely used for this kind of problem could be inefficient because of their deterrent computation time [6]. That is why, to achieve our goal, it is necessary to adopt new methods that allows, both the integration of control parameters and the design optimization with better computation time. In this paper, we herein propose a realistic model of an electrical motor, and then we formulated an optimal control problem with the willing to propose a global optimum.

The optimal control theory is a branch of variation calculation. It appeared after the Second World War to solve practical problems especially in the field of aeronautic. The key of this theory is the Pontryagin’s Maximum Principle, PMP, introduced by L. S. Pontryagin in 1956. The PMP provides a necessary condition for optimality and allows the calculation of optimal trajectories [7,8]. This theory is applied in various fields. e.g. in [9], the PMP is used to solve problems of orbital transfer. In [10] and [11], the control theory is applied to find strategies of energetic optimization of hybrid electrical vehicles. Investigations in [12] used the theory in order to find an optimal control strategy that minimizes the flux of power in a small electrical network. However, applying optimal control theories on electromagnetic devices is not frequent, because of devices models that are often complex and non-linear [13]. Nevertheless, some investigations get interested on the application of the PMP for the study of electrical machines, with using first order machine models. In [4], e.g. a position control algorithm is developed while minimizing time and energy of an induction motor.

The present work applied the theory of optimal control on a magnetodynamic model based on a reluctance network model of a hybrid stepper motor, HSM. The main purpose is to move the load to a certain position while minimizing energy. HSM motors are characterized by a tridimensional distribution of flux. The rotor corporates a permanent magnet axially magnetized and magnetic teeth. The stator carries the winding and has small teeth at the periphery of the bore. Varieties of models of HSM were developed in many researches. In [14], Matsui and al elaborated one of the first examples of models. These models are based on an equivalent magnetic circuit composed of permeances for the air gap, magneto mortice forces for the coils and of the magnet. The model is an entirely analytical and the relative permeability of iron is assumed infinite. Kang and Lieu [15] presented another approach. In their research, they calculated the air gap permeances using the 2D finite elements method. A deeper study was presented by Kuet and al in [16], were the elaborated magnetic model take in account the saturation phenomena and the leakage flux. From 2002, advances on the numerical calculations have made the study of 3D finite of complex structures possible. In [17], Koasaka and al, presented a 3D finite elements model of an HSM. Investigation in [18] by Stuebig and Ponick proposes a comparison between different modeling methods.

The model of electrical machine used in the present study is based on Matsui and al [14] work, with an improvement of the accuracy of air gap permeances calculation. The model choice is built over its compromise between the accuracy and calculation time.

The paper is organized as follows; Section 2 presents the chosen model structure under the form of state variable. In Section 3, the PMP is applied to find the optimal control profile that minimizes the necessary energy for displacement. A classical control obtained with an optimization on a step is illustrated in Section 4. Finally, the last section is dedicated to compare and analyze obtained results of two controls.

2. Modeling

Figure 1 illustrates a section of the hybrid stepper motor structure. It presents 8 powered stator’s plot and has two phases 𝛼 and 𝛽. Each stud has 6 small teeth. The rotor is composed of two rings; each one has fifty teeth and are connected between them by a permanent magnet. The flux path is then tridimensional in the rotor.

Fig. 1.

Structure of the hybrid stepper motor.

A 3D magnetic model of reluctance network was developed in [17]. This model computes inductances L_𝛼(θ) and L_𝛽(θ), flux φ_𝛼(θ) and φ_𝛽(θ) provided by the magnet. These values are calculated for each mechanical position of the rotor θ depending of dimensions and machine’s parameters. These values are injected in the electrical equations, in order to determine phases currents I_𝛼(t) and I_𝛽(t). The electromagnetic torque C_em(t) is deduced and is applied further in the magnetic model. Figure 2 shows the magnetodynamic global model and relations between models. This model is parametrized and contains the coupling of 3 physics: electrical, magnetic and mechanical. With l_m: length of the magnet (m), l_r: length of a rotor ring (m), g: gap length (m), t_w: tooth width (m), N_s: number of turns, N_t: number of tooth, and H_C: coercive field of the magnet (A/m). θ(rd): electric angle, k_fv: viscous friction coefficient, J: inertia of the system (kg⋅m²).

Fig. 2.

Electromechanical model of the stepper motor.

Thus, the dynamic model of the device could be expressed as in (1): $\begin{eqnarray}\displaystyle \left(\begin{array}{@{}c@{}}\displaystyle {\displaystyle \frac{dI_{{\alpha}}}{dt}}\\[11.0pt] \displaystyle \frac{dI_{{\beta}}}{dt}\\[11.0pt] \displaystyle \frac{d{\Omega}}{dt}\\[11.0pt] \displaystyle \frac{d{\theta}}{dt}\end{array}\right)=\left(\begin{array}{@{}c@{}}\displaystyle -{\displaystyle \frac{1}{L_{{\alpha}}}}I_{{\alpha}}\left(R+{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}{\Omega}\right)-{\displaystyle \frac{1}{L_{{\alpha}}}}{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}{\Omega}\\[13.0pt] \displaystyle -{\displaystyle \frac{1}{L_{{\beta}}}}I_{{\beta}}\left(R+{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}{\Omega}\right)-{\displaystyle \frac{1}{L_{{\beta}}}}{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}{\Omega}\\[13.0pt] \displaystyle \frac{1}{J}\left(\frac{1}{2}I_{{\alpha}}^{2}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}+\frac{1}{2}I_{{\beta}}^{2}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}+I_{{\alpha}}{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}+I_{{\beta}}{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}-k_{fv}{\Omega}\right)\\[5.0pt] {\Omega}\end{array}\right)+\left(\begin{array}{@{}cc@{}}\displaystyle \frac{1}{L_{{\alpha}}} & 0\\[10.0pt] 0 & \displaystyle \frac{1}{L_{{\beta}}}\\[12.0pt] 0 & 0\\ 0 & 0\end{array}\right)\left(\begin{array}{@{}c@{}}V_{{\alpha}}\\[4.0pt] V_{{\beta}}\end{array}\right). & & \displaystyle\end{eqnarray}$ (1) From the supply V_𝛼 and V_𝛽, the 4 state variables of the system differential equation could be obtained.

3. Optimal energy

3.1. Optimization problem

Starting from an initial position θ₀, a final position θ_f should be reached with a final speed 𝛺_f in an instant T. This displacement has to be operated with the minimum of energy. The problem is therefore given by ((2)). $$\begin{eqnarray} \displaystyle \begin{array}{@{} l@{}}\displaystyle \min _{V_{{\alpha}},V_{{\beta}}}\mathrm{Obj}(t)=\int _{0}^{T}(|V_{{\alpha}}|+|V_{{\beta}}|)∼\mathrm{d}t\\ [6.0pt] \text{With:}\\ [4.0pt] \qquad V_{{\alpha}\min }\leq V_{{\alpha}}\leq V_{{\alpha}\max }\\ [4.0pt] \qquad V_{{\beta}\min }\leq V_{{\beta}}\leq V_{{\beta}\max }\\ [6.0pt] {\theta}(t_{0})={\theta}_{0},\quad {\theta}(T)={\theta}_{f}\\ [4.0pt] {\Omega}(t_{0})={\Omega}_{0},\quad {\Omega}(T)={\Omega}_{f}\\ [4.0pt] I_{{\alpha}}(t_{0})=I_{{\alpha}0},\quad I_{{\alpha}}(T)=I_{{\alpha}f}\\ [4.0pt] I_{{\beta}}(t_{0})=I_{{\beta}0},\quad I_{{\beta}}(T)=I_{{\beta}f}. \end{array}& & \displaystyle \end{eqnarray}$$(2) The result of Eqs ((2)) presents voltages trajectories. Voltage values for each instant for V_𝛼(t) and V_𝛽(t). With this formulation, the hypothesis is that the decrease of the control variable, which means the voltage, produces a decrease of currents and then energy. This is valid at low speeds [4]. This formulation allows avoiding stiff phenomenon while resolving the Hamiltonian system.

3.2. Hamiltonian model

The Hamiltonian of system (1) is given by: $\begin{eqnarray}\displaystyle H\text{} & =\text{} & \displaystyle \displaystyle (|V_{{\alpha}}|+|V_{{\beta}}|)-{\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}\left(RI_{{\alpha}}+I_{{\alpha}}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}{\Omega}+{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}{\Omega}-V_{{\alpha}}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \displaystyle -∼{\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}\left(RI_{{\beta}}+I_{{\beta}}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}{\Omega}+{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}{\Omega}-V_{{\beta}}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \displaystyle +∼{\psi}_{3}{\displaystyle \frac{1}{J}}\left({\displaystyle \frac{1}{2}}I_{{\alpha}}^{2}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}+{\displaystyle \frac{1}{2}}I_{{\beta}}^{2}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}+I_{{\alpha}}{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}+I_{{\beta}}{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}-k_{fv}{\Omega}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼{\psi}_{4}{\Omega}\end{eqnarray}$ (3) 𝜓₁, 𝜓₂, 𝜓₃ and 𝜓₄ are costates variables that provide the necessary conditions for the optimality. Thus, the optimal control problem can be stated as follows: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}\displaystyle \frac{d{\psi}_{1}}{dt}=-{\displaystyle \frac{\partial H}{\partial I_{{\alpha}}}}\\[10.0pt] \displaystyle \frac{d{\psi}_{2}}{dt}=-{\displaystyle \frac{\partial H}{\partial I_{{\beta}}}}\\[10.0pt] \displaystyle \frac{d{\psi}_{3}}{dt}=-{\displaystyle \frac{\partial H}{\partial {\Omega}}}\\[10.0pt] \displaystyle \frac{d{\psi}_{4}}{dt}=-{\displaystyle \frac{\partial H}{\partial {\theta}}}.\end{array} & & \displaystyle\end{eqnarray}$ (4)

According to the problem formulation, the optimal control must minimize the Hamiltonian of the system. Thus, we should study the variation of the Hamiltonian with respect to each control, and then we pass through the derivative of H according to each control, as H_{V 𝛼} and H_{V 𝛽}. $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle H_{V{\alpha}}={\displaystyle \frac{\partial H}{\partial V_{{\alpha}}}}={\displaystyle \frac{V_{{\alpha}}}{|V_{{\alpha}}|}}+{\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}\end{eqnarray}$ (5) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle H_{V{\beta}}={\displaystyle \frac{\partial H}{\partial V_{{\beta}}}}={\displaystyle \frac{V_{{\beta}}}{|V_{{\beta}}|}}+{\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}.\end{eqnarray}$ (6) The study of the derivative sign enables finding the minimum of the Hamiltonian. E.g. for 𝜓₁(1∕L_𝛼) > 1, H_{V 𝛼} will be positive, in this case, H (V_𝛼) is strictly increasing and the minimum of the Hamiltonian is reached at V_𝛼 = V_𝛼min. The application of Pontryagin’s maximum principle enables us to determine the control in terms of costate variables 𝜓₁(t) and 𝜓₂(t). After analysis, (7) and (8) are subtracted: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle V_{{\alpha}}^{\ast }=\left\{\begin{array}{@{}ll@{}}\displaystyle -V_{{\alpha}\max }\text{sign}\left({\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}\right) & \displaystyle \text{if}∼\left|{\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}\right|>1\\[12.0pt] 0 & \displaystyle \text{if}∼\left|{\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}\right|\leq 1\end{array}\right.\end{eqnarray}$ (7) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle V_{{\beta}}^{\ast }=\left\{\begin{array}{@{}ll@{}}\displaystyle -V_{{\beta}\max }\text{sign}\left({\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}\right) & \displaystyle \text{if}∼\left|{\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}\right|>1\\[12.0pt] 0 & \displaystyle \text{if}∼\left|{\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}\right|\leq 1.\end{array}\right.\end{eqnarray}$ (8)

Equations (7) and (8) are called switching functions that switch according to costate variables. These types of controls are called bang-off-bang controls, while they assume values only in the boundary points of the control interval. This result is interesting for the study of electrical machines. In fact, we firstly proved that for minimizing the energetic consumption, we must deal with bang-off-bang controls. In a second row, this result is obtained after explicitly expressing the control, considering parameters in Fig. 2.

Equations (1), (4), (7) and (8) gather the Hamiltonian model composed of the costate model (9) and of the state model (10). The link between the two models is obtained with the control, expressed in (7) and (8). We must emphasize here that one of the drawbacks of this method is the heaviness of the analytical calculations. $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\displaystyle \frac{d{\psi}_{1}}{dt}}={\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}\left(R+{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}{\Omega}\right)-{\psi}_{3}{\displaystyle \frac{1}{J}}\left(I_{{\alpha}}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}+{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle {\displaystyle \frac{d{\psi}_{2}}{dt}}={\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}\left(R+{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}{\Omega}\right)-{\psi}_{3}{\displaystyle \frac{1}{J}}\left(I_{{\beta}}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}+{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle {\displaystyle \frac{d{\psi}_{3}}{dt}}={\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}\left(I_{{\alpha}}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}+{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}\right)+{\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}\left(I_{{\beta}}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}+{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}\right)+{\psi}_{3}{\displaystyle \frac{k_{fv}}{J}}-{\psi}_{4}\nonumber\end{eqnarray}$ $\begin{eqnarray}\displaystyle {\displaystyle \frac{d{\psi}_{4}}{dt}}\text{} & =\text{} & \displaystyle {\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}{\Omega}\left(-I_{{\alpha}}{\displaystyle \frac{1}{L_{{\alpha}}}}\left({\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}\right)^{2}+I_{{\alpha}}{\displaystyle \frac{\partial ^{2}L_{{\alpha}}}{\partial {\theta}^{2}}}-{\displaystyle \frac{1}{L_{{\alpha}}}}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}+{\displaystyle \frac{\partial ^{2}{\varphi}_{{\alpha}}}{\partial {\theta}^{2}}}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼{\psi}_{1}{\displaystyle \frac{1}{L_{{\alpha}}}}\left(-RI_{{\alpha}}{\displaystyle \frac{1}{L_{{\alpha}}}}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}+V_{{\alpha}}^{\ast }{\displaystyle \frac{1}{L_{{\alpha}}}}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼{\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}{\Omega}\left(-I_{{\beta}}{\displaystyle \frac{1}{L_{{\beta}}}}\left({\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}\right)^{2}+I_{{\beta}}{\displaystyle \frac{\partial ^{2}L_{{\beta}}}{\partial {\theta}^{2}}}-{\displaystyle \frac{1}{L_{{\alpha}}}}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}+{\displaystyle \frac{\partial ^{2}{\varphi}_{{\beta}}}{\partial {\theta}^{2}}}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼{\psi}_{2}{\displaystyle \frac{1}{L_{{\beta}}}}\left(-RI_{{\beta}}{\displaystyle \frac{1}{L_{{\beta}}}}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}+V_{{\beta}}^{\ast }{\displaystyle \frac{1}{L_{{\beta}}}}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}\right)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle -∼{\psi}_{3}{\displaystyle \frac{1}{J}}\left({\displaystyle \frac{1}{2}}I_{{\alpha}}^{2}{\displaystyle \frac{\partial ^{2}L_{{\alpha}}}{\partial {\theta}^{2}}}+{\displaystyle \frac{1}{2}}I_{{\beta}}^{2}{\displaystyle \frac{\partial ^{2}L_{{\beta}}}{\partial {\theta}^{2}}}+I_{{\alpha}}{\displaystyle \frac{\partial ^{2}{\varphi}_{{\alpha}}}{\partial {\theta}^{2}}}+I_{{\beta}}{\displaystyle \frac{\partial ^{2}{\varphi}_{{\beta}}}{\partial {\theta}^{2}}}\right)\end{eqnarray}$ (9) $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}\displaystyle {\displaystyle \frac{dI_{{\alpha}}}{dt}}=-{\displaystyle \frac{1}{L_{{\alpha}}}}I_{{\alpha}}\left(R+{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}{\Omega}\right)-{\displaystyle \frac{1}{L_{{\alpha}}}}{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}{\Omega}+{\displaystyle \frac{1}{L_{{\alpha}}}}V_{{\alpha}}^{\ast }\\[14.0pt] \displaystyle {\displaystyle \frac{dI_{{\beta}}}{dt}}=-{\displaystyle \frac{1}{L_{{\beta}}}}I_{{\beta}}\left(R+{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}{\Omega}\right)-{\displaystyle \frac{1}{L_{{\beta}}}}{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}{\Omega}+{\displaystyle \frac{1}{L_{{\beta}}}}V_{{\beta}}^{\ast }\\[14.0pt] \displaystyle {\displaystyle \frac{d{\Omega}}{dt}}={\displaystyle \frac{1}{J}}\left({\displaystyle \frac{1}{2}}I_{{\alpha}}^{2}{\displaystyle \frac{\partial L_{{\alpha}}}{\partial {\theta}}}+{\displaystyle \frac{1}{2}}I_{{\beta}}^{2}{\displaystyle \frac{\partial L_{{\beta}}}{\partial {\theta}}}+I_{{\alpha}}{\displaystyle \frac{\partial {\varphi}_{{\alpha}}}{\partial {\theta}}}+I_{{\beta}}{\displaystyle \frac{\partial {\varphi}_{{\beta}}}{\partial {\theta}}}\right)\\[12.0pt] \displaystyle {\displaystyle \frac{d{\theta}}{dt}}={\Omega}.\end{array} & & \displaystyle\end{eqnarray}$ (10) After that, we numerically solved the Hamiltonian system (9) and (10) using the single shooting method, as the problem is a two-point boundary value problem for which both the initial states of the state variables (10) and the final states of the co-state variables (9) are known. A function S, called shooting function, is defined in (11): $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle S:\left(\begin{array}{@{}c@{}}{\psi}_{1}(t_{0})\\[3.0pt] {\psi}_{2}(t_{0})\\[3.0pt] {\psi}_{3}(t_{0})\\[3.0pt] {\psi}_{4}(t_{0})\end{array}\right)\quad \stackrel{(\text{ode15s})}{\overrightarrow{\phantom{(ode15s)}}}\quad \left(\begin{array}{@{}c@{}}I_{{\alpha}}(T)\\[3.0pt] I_{{\beta}}(T)\\[3.0pt] {\Omega}(T)\\[3.0pt] {\theta}(T)\end{array}\right)\end{eqnarray}$ (11) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \text{Find}∼\left(\begin{array}{@{}c@{}}{\psi}_{1}(t_{0})\\[3.0pt] {\psi}_{2}(t_{0})\\[3.0pt] {\psi}_{3}(t_{0})\\[3.0pt] {\psi}_{4}(t_{0})\end{array}\right)\quad \text{as:}∼S\left(\begin{array}{@{}c@{}}{\psi}_{1}(t_{0})\\[3.0pt] {\psi}_{2}(t_{0})\\[3.0pt] {\psi}_{3}(t_{0})\\[3.0pt] {\psi}_{4}(t_{0})\end{array}\right)=\left(\begin{array}{@{}c@{}}I_{{\alpha}}(T)\\[3.0pt] I_{{\beta}}(T)\\[3.0pt] {\Omega}(T)\\[3.0pt] {\theta}(T)\end{array}\right).\end{eqnarray}$ (12) This method is accurate but its major inconvenient is its high sensitivity when it comes to the choice of initial conditions [7,8]. In the case of an electrical machine, all chosen initial conditions are naturally the solution of a positioning problem, defined to reach progressively the desired solution. Thus, the solution of each sub-problem will be the initial condition of the next problem, which often allow the convergence of shooting method. In the case of an electrical machine, we propose an easier and faster solution; a dichotomy technique is applied on the initial condition of the most influential costate variable. In our case, it’s the variable related to the position 𝜓₄(t₀). This resolution method is well adapted for the study of electrical machines and allows a rapid achievement of a solution. Figure 3 shows the integration of the Hamiltonian model and the methodology of resolution.

Fig. 3.

Resolution process with the Hamiltonian model.

3.3. Results

The method is applied on a small hybrid stepper motor of 200 steps per spin (1 step = 0.0314 rd). The machine is at standstill and must drive a load, reach a position with an accuracy of 10⁻⁴ rd, in this example, the angle 0.6776 rd and must be at standstill in 1.3 s. The problem resolution takes about 10 min with an ordinary computer. Graphs in Figs 4 and 5 illustrate respectively, the obtained optimal control $V_{{\alpha}}^{\ast }$ and $V_{{\beta}}^{\ast }$ and the corresponding switching functions.

Fig. 4.

Optimal control $V_{{\alpha}}^{\ast }$ .

Fig. 5.

Optimal control $V_{{\beta}}^{\ast }$ .

The optimal control generated is a bang-off-bang control with different pulse widths that decrease progressively until reaching the desired final position. In Table 1 different widths values of positive and negative pulse of the control $V_{{\alpha}}^{\ast }$ are given:

Table 1

Pulse width

Bandwidth	Value (s)
t1	0.1082
t2	0.1060
t3	0.1052
t4	0.1038
t5	0.1013
t6	0.0979
t7	0.0943
t8	0.0903
t9	0.0878
t10	0.0842
t11	0.0815

The motor makes approximately 21 to 22 steps. This optimal number of steps and optimal phase switching are generated automatically.

Figure 6 shows the rotation speed and the evolution of the mechanical angle.

Fig. 6.

Evolution of rotation speed and mechanical angle.

After the start-up, the speed varies between 0.67 rd/s and 0.33 rd/s, the switches of phases and the evolution of steps cause the drops of speed. The mean speed over a step increases progressively during the first 13 steps then drops to be equal to zero on the final position. This behavior is caused by the inertia effect that exert a high impact on the energy consumption. Thanks to the inertia, the motor needs less and less energy to reach the final position, which elucidate the decrease of pulse widths, given in Table 1 and the shape of the evolution of the position that illustrates a curvature.

Figure 7 outlines the evolution of costate variables 𝜓₁(t) and 𝜓₂(t), respectively related to currents I_𝛼(t) and I_𝛽(t). We naturally find the same shapes and the same phase shift as the switching functions with a difference in amplitudes due to à 1∕L_𝛼 and 1∕L_𝛽 beta (5), (6).

Fig. 7.

Costate variables 𝜓₁(t) and 𝜓₂(t).

Figure 8 illustrates the evolution of the costate variables 𝜓₃(t) and 𝜓₄(t) related to the speed 𝛺(t) and the position θ(t), respectively.

Fig. 8.

Costate variables 𝜓₃(t) and 𝜓₄(t).

Figure 9 shows the evolution of the currents I_𝛼(t) and I_𝛽(t).

Fig. 9.

Evolution of the excitation current.

3.4. Dichotomy

A dichotomy method applied to the costate variable according to the position 𝜓₄(t₀) enable reaching rapidly the desired position with respect to the accuracy. This method is fast, easy to apply and avoids problem of convergence of the shooting method. Figure 10 depicts the functions of switch and the evolution of the mechanical angle for different values of 𝜓₄(t₀).

Fig. 10.

Switch function and mechanical angle for different value of 𝜓₄(t₀).

The method of dichotomy on the initial value of 𝜓₄(t₀) enable the settle of the accuracy on the final value of the position. In our considered example, for a value of 𝜓₄(t₀) = −2005.8 allows reaching the final position of 0.6776 rd with accuracy and less iteration.

4. Classic control

The purpose is to compare the obtained results with the PMP to a classical control. In the case of a classical control, pulse widths are unchanged. To find a classical control that enable reaching the desired position with a minimum of energy, we proceed by using time-scaling method [19]:

We consider x₁ the width of each pulse, for the two controls, and x₂ (resp. x₃) the width of the last pulse of the control V_𝛼(t) (resp. V_𝛽(t)). We aim to find x₁, x₂, and x₃ which allow us to achieve an angle 0.6776 rd and to standstill in 1.3 s with the minimum of energy. The problem could be formulated as in ((13)): $$\begin{eqnarray} \displaystyle \begin{array}{@{} l@{}}\displaystyle \min _{x1,x2,x3}\text{Obj}(t)=\int _{0}^{T}(|V_{{\alpha}}|+|V_{{\beta}}|)∼\mathrm{d}t\\ [10.0pt] \text{With:}\\ [4.0pt] \qquad x_{1\min }\leq x_{1}\leq x_{1\max }\\ [4.0pt] \qquad x_{2\min }\leq x_{2}\leq x_{2\max }\\ [4.0pt] \qquad x_{3\min }\leq x_{3}\leq x_{3\max }\\ [6.0pt] {\theta}(t_{0})={\theta}_{0},\quad {\theta}(T)={\theta}_{f}\\ [4.0pt] {\Omega}(t_{0})={\Omega}_{0},\quad {\Omega}(T)={\Omega}_{f}\\ [4.0pt] I_{{\alpha}}(t_{0})=I_{{\alpha}0},\quad I_{{\alpha}}(T)=I_{{\alpha}f}\\ [4.0pt] I_{{\beta}}(t_{0})=I_{{\beta}0},\quad I_{{\beta}}(T)=I_{{\beta}f}. \end{array}& & \displaystyle \end{eqnarray}$$(13)

The used solver is the fmincon function of Matlab^® with the algorithm SQP. The pulse widths obtained by optimization are given in Table 2.

Table 2

Pulse widths

Bandwidth	Value (s)
x1	0.1080
x2	0.0825
x3	0.1192

Figure 11 shows the classical controls V_𝛼(t) and V_𝛽(t) obtained with the optimization by ((13)).

Fig. 11.

Classic optimal control.

5. Comparison of results

In Fig. 12 are illustrated the rotation speed and the evolution of the mechanical angle for the two types of control.

Fig. 12.

Evolution of rotation speed and mechanical angle for both methods.

The classic control enables reaching the position faster with a more regular speed evolution. Figure 13 presents the evolution of the voltages sum at each time instant for the two cases of control. The value evolve between 2V_max and V_max for the control obtained by the PMP, while for the classic control, the value stand constant at 2V_max at each step from the start, thus each step consumes the same amount of energy.

Fig. 13.

Evolution of objective function.

As shown in Fig. 14 of the evolution of the power developed by the two controls, we notice the effect of the inertia and the gradually decrease of the power consumption for the optimal control. The machine develops then the minimum of torque necessary for driving the load to the final position. Figure 15 represents the electromagnetic torque obtained for the two types of controls. For the classic control the same torque is developed for each step, whereas, for the optimal control the torque decreases progressively.

Fig. 14.

Evolution of input power.

Fig. 15.

Evolution of electromagnetic torque.

Table 2 introduces a comparison between the energetic performances of the HSM for the types of control.

Table 3

Energetic performance comparison for the classic control and optimal control

	Classic control	Optimal control	Gain %
Objective function	812 535	765 860	5.75%
Input power	21.094	20.025	5.06%

The optimal control obtained with the PMP offered a gain of 5.7% compared to the classical control for the studied problem.

6. Conclusion

In this paper, we proposed the application of the Pontryagin’s maximum principle of to a magneto-dynamic model based on a reluctance network of a hybrid stepper motor. The main purpose was to find an optimal control that minimizes the energy. The applying of PMP enable to express explicitly the control according to the parameters of the model and to have a Hamiltonian model that generates automatically an optimal control. The study has showed that the problem of sensitivity to the initial conditions encountered at the resolution of the Hamiltonian model could be solved by a simple dichotomy when it comes to the control problem of an electrical machine. The presented results showed a gain of 5.7% compared to a classic method thanks to the effect of inertia that has an important impact on the energy consumption and that is present in all the electrical systems.

References

Lei

Zhu

and Guo

, Multydisciplinary Design Optimisation Methods for Electrical Machines and Drive Systems, Springer-Verlag, Berlin, Heidelberg, 2016.

Dorrell

, A review of the methods for improving the efficiency of drive motors to meet IE4 efficiency standards, Journal of Power Electronics14 (2014), 842–851.

Gillon

and Brochet

, Shape optimization of a permanent magnet motor using the experimental design method, IEEE Trans. on Magn.35(3) (1999), 1278–1281.

Alihalli

and Ilyas Bayindir

, Time-energy optimal control of vector controlled induction motor, COMPEL: Int. J. Computation and Mathematics in Electrical and Electronic Engineering21 (2002), 235–251.

Caillard

Gillon

Randi

and Janiaud

, Mono and bi-level optimization architectures for powertrain design, The International Journal for Computation and Mathematics in Electrical and Electronic Engineering (COMPEL)35(3) (2016), 847–859.

Mohammadi

M.H.

Silva

R.C.P.

and Lowther

D.A.

, Incorporating control strategies into the optimization of synchronous AC machines: A comparison of methodologies, IEEE Transaction on Magnetics54(3) (2018), 1–4, doi:10.1109/TMAG.2017.2748099.

Trélat

, Controle Optimal: Théorie et Applications, Mathématiques concrètes, 2005.

Geering

, Optimal control with engineering application, IEEE Control Systems31(5) (2011), 115–117.

Zhu

Trélat

and Cerf

M.P.

, Geometric optimal control and applications to aerospace, J. Math. Ind.9 (2017), 2198–4115.

10.

Uebel

Murgovski

Tempelhahn

and Bker

, Optimal energy management and velocity control of hybrid electric vehicles, IEEE Transactions on Vehicular Technology67(1) (2018), 327–337.

11.

Kim

Cha

S.-W.

and Peng

, Optimal control of hybrid electric vehicles based on Pontryagin’s minimum principle, IEEE Trans. Contr. Sys. Technol.19 (2011), 1279–1287, doi:10.1109/TCST.2010.2061232.

12.

Dagdougui

, Optimal control of a network of power microgrids using the Pontryagin’s minimum principle, IEEE Transactions on Control Systems Technology22(5) (2014), 1942–1948.

13.

Gillet

David

M.P.

and Messine

, Optimization of the control of a doubly fed induction machine, in: 2013 IEEE 11th International Workshop of Electronics, Control, Measurement, Signals and their Application to Mechatronics, Toulouse, 2013, pp. 1–5.

14.

Matsuri

Nakamura

and Kosaka

, Instantaneous torque analysis of hybrid stepping motor, IEEE Transactions on Industry Applications34(5) (1996), 1176–1182.

15.

Kang

S.G.

and Lieu

, Torque analysis of combined 2D FEM and lumped parameter method for a hybrid stepping motor, in: Proc. IEEE Int. Conf. Elect. Mach. Drives, 2005, pp. 1199–1203.

16.

Kuert

Jufer

and Perriard

, New method for dynamic modeling of hybrid stepping motors, in: Conf. Rec. 37th IEEE IAS Annu. Meeting, Vol. 1, 2002, pp. 6–12.

17.

Kosaka

Pollock

and Matsui

, 3 dimensional finite element analysis of hybrid stepping motors taking inter-lamination gap into account, in: Proc. End Int. Conf. PEMD, Vol. 2, 2004, pp. 534–539.

18.

Stuebig

and Ponick

, Comparision of calculation methods for hybrid stepping motors, IEEE Transaction on Industry Applications48(6) (2012), 2182–2189.

19.

Jiang

Gao

Wang

and Xiao

, A moving time node method for solving robot optimal control problem, in: 2019 Chinese Control Conference (CCC), Guangzhou, China, 2019, pp. 1899–1904.