The main objective and contribution of this work is the development of a decentralized RMPC for virtual coupling in railways, considering the effects of adhesion loss and positioning errors. The research in [19] serves as a starting point for this work, from which the design of a nominal MPC is taken. Unlike this previous work, the implemented robust MPC considers the external disturbances with uncertainties that consider the previously mentioned positioning errors and a loss of the adhesion during braking. In this way, the proposed controller is compared with a nominal MPC in different simulations under the different applied disturbances. The results will be studied from an operational point of view, i.e., by analysing the distance between the trains and the overall behaviour of a two-train VCTS.

Ultimately, the objective of this work is to improve the line capacity under the condition of virtual coupling. In this sense, we consider that the VCTS fulfils its function and that all the components of the VCTS work as connected when all of them arrive and stop at the same time at the same station with a minimum delay time between them. In this way, the convoy operationally behaves as a single train since, from the users' point of view, all the VCTSs' components arrive at the station at the same time and behave as if they were a single complete train. Although when the convoy is running the distance between the components increases due to the safety conditions which have been set, we consider that the convoy preserves its integrity, operating as a single train, when the communication between trains is maintained and the above-mentioned condition of all the components arriving at the station at the same time is respected. The figures included in the simulations section illustrate this phenomenon.

The remainder of this paper is organized as follows. Section 2 presents the dynamic model. Section 3 formulates the design of the RMPC. Section 4 describes the simulation cases developed and presents their corresponding results. Finally, Section 5 includes the conclusions of this work.

### **2. System Dynamics**

The model defining the train motion of this work is based on the principles of longitudinal train dynamics (LTD). Then, the train is considered as a point mass with one degree of freedom, where the traction/brake system, rolling resistances, air intake, aerodynamic

drag, and slope and curve resistances are applied. The dynamic equations considered are as follows: . ̇ = ( − ) ⁄ (3) where *s (m)* and *v (m/s)* denote the position and train speed, *u (N)* is the controlled

The model defining the train motion of this work is based on the principles of longitudinal train dynamics (LTD). Then, the train is considered as a point mass with one degree of freedom, where the traction/brake system, rolling resistances, air intake, aerodynamic drag, and slope and curve resistances are applied. The dynamic equations

*Actuators* **2022**, *11*, x FOR PEER REVIEW 4 of 18

̇ = (− − −

**2. System Dynamics**

considered are as follows:

$$
\dot{s} = v \tag{1}
$$

̇ = (1)

2 − + )/ + (2)

$$\dot{\upsilon} = \left(-a - b\upsilon - c\upsilon^2 - F\_\varepsilon + F\right) / M + w^a \tag{2}$$

$$
\dot{\mathbf{F}} = (\mathbf{u} - \mathbf{F}) / \tau \tag{3}
$$

where *s* (*m*) and *v* (*m*/*s*) denote the position and train speed, *u* (*N*) is the controlled driving/braking force, *F* (*N*) is the integrated driving/braking force, *Fe* (*N*) is the resistance force due to the track, *τ* is the inertial lag of the longitudinal dynamics, *M* (*kg*) is the train's mass, *a* (*N*) is a term that includes the rolling resistance plus the bearing resistance, *b* (*Ns*/*m*) is a coefficient related to the air intake, *c* (*Ns*2/*m*2 ) is the aerodynamic coefficient, and *w a* (*m*/*s* 2 ) represents the uncertainty in the acceleration contemplated in the robust control. coefficient, and *(m/s2 )* represents the uncertainty in the acceleration contemplated in the robust control. The resistance force includes two components, and , which are defined as follows: = + (4)

The resistance force *Fe* includes two components, *Fg* and *FR*, which are defined as follows: = − × (5)

$$F\_{\varepsilon} = F\_{\mathcal{S}} + F\_{\mathcal{R}} \tag{4}$$

$$F\_{\mathcal{S}} = -M\mathfrak{g} \times slope \tag{5}$$

$$F\_{\mathbb{R}} = -M \times 6/\mathbb{R} \tag{6}$$

where *Fg* (*N*) is the component of the gravity force due to the slope of the track, the *slope* (m⁄m) is the slope of the track, g *m*/*s* 2 is the acceleration of gravity, *FR* (*N*) is the resistance in the curve, and *R* (*m*) is the radius of the curve. (m⁄m) is the slope of the track, g (/ ) is the acceleration of gravity, () is the resistance in the curve, and () is the radius of the curve. The values of the *slope* and *R* depend on the line profile and the position *s* of the train

2

The values of the *slope* and *R* depend on the line profile and the position *s* of the train on the line and are, therefore, known at each time. Figure 1 shows as an example of the value of the *slope* and *R* used in this work. on the line and are, therefore, known at each time. Figure 1 shows as an example of the value of the *slope* and *R* used in this work.

**Figure 1.** Slope and radius considered in the simulation scenario. **Figure 1.** Slope and radius considered in the simulation scenario.

defined as follows:

Starting from the dynamic equations, a linearized set of equations at the operating Starting from the dynamic equations, a linearized set of equations at the operating point will be used to handle uncertainties in the RMPC optimization algorithm [20].

point will be used to handle uncertainties in the RMPC optimization algorithm [20]. Thus, the vector of states *X* and outputs *Y*, that will be used for a linearization, are Thus, the vector of states *X* and outputs *Y*, that will be used for a linearization, are defined as follows:

$$\mathbf{X} = \begin{bmatrix} s \\ v \\ F \end{bmatrix} \tag{7}$$

$$\mathcal{Y} = \begin{bmatrix} \mu \cdot v \\ v^2 / 2a^e \end{bmatrix} \tag{8}$$

where *a e* (*m*/*s* 2 ) is the maximum possible deceleration of the controlled train.

By linearizing Equations (1)–(3), these equations are transformed into (9)–(11), thereby formulating the linearized dynamics as the first-order Taylor polynomial of the nonlinear model as follows:

$$\mathcal{U} = \begin{bmatrix} u \end{bmatrix}, \quad \mathcal{W} = \begin{bmatrix} 0 \\ w^a \\ 0 \end{bmatrix} \tag{9}$$

−1 ()

$$\dot{X} = \dot{X}\_0 + A(X - X\_0) + B(\mathcal{U} - \mathcal{U}\_0) + \mathcal{W} \tag{10}$$

$$Y = Y\_0 + \mathcal{C}(X - X\_0) + D(\mathcal{U} - \mathcal{U}\_0) \tag{11}$$

and

$$\begin{aligned} A &= \begin{bmatrix} 0 & 1 & 0 \\ 0 & -(B + 2Cv\_0)/M & 1/M \\ 0 & 0 & -1/\tau \end{bmatrix} & B &= \begin{bmatrix} 0 \\ 0 \\ 1/\tau \end{bmatrix} \\\\ \mathbf{C} &= \begin{bmatrix} 0 & u\_0 & 0 \\ 0 & v\_0/a\_f & 0 \end{bmatrix} & D &= \begin{bmatrix} v\_0 \\ 0 \end{bmatrix} \end{aligned}$$

where *v*0 (*m/s*) and *u*0 (*N*) corresponds to the state *v* and the control force *u* at the operating point, respectively; *X*0, *U*0, and *Y*0 are the state vector *X*, the decision vector *U*, and the output vector *Y* at the operating point, respectively; and . *X*0 is Equations (1)–(3) evaluated at the operating point. *Actuators* **2022**, *11*, x FOR PEER REVIEW 6 of 18

> In the above Equations (1)–(11), the states *s*, *v*, *F*; the controlled driving/braking force *u*; and the uncertainty *w* are time dependent. Therefore, *X* and *Y* are also time dependent. The matrices *A*, *B*, *C*, *D* are constant, but *v*0 and *u*0 depend on the operating point at which the equations are linearised. In addition, the values of *slope* and *R* are known at each time and dependent on the line profile and train location on the line, i.e., the train position *s*, as shown in Figure 1. The remaining variables in (1)–(11) are constant parameters that depend on the train's characteristics. described in Section 3.5. For the followers, we design an MPC robust controller that ensures a safety and control efficiency while considering the parameter uncertainties. With the proposed control strategy, the leader will track a given speed curve, and the followers will guarantee a safe minimum distance between trains and the string stability

#### **3. MPC Controller Design** of the VCTS with parameter uncertainties. Consequently, as seen in Figure 2, the follower

#### *3.1. Control Architecture* needs information about the preceding train to ensure its safety. This information (

We consider a convoy (Figure 2) composed of a leader and *n* followers, all of them with length *L*. The superscript indicates the train, with 0 for the leader and *i,* where *i* = 1, . . . , *n*, for the followers. and ) is used by the follower's controller for its calculations because the predicted states for the preceding train − 1 are calculated using the PTP embedded in each follower based on the information received from the preceding train.

Data communication

**Figure 2.** Control architecture. **Figure 2.** Control architecture.

ETCS.

*3.2. Leader Controller* When using a decentralized control architecture, the leader can operate under any control and signalling method. Therefore, the leader can use any conventional control We consider a decentralized VCTS control problem, with independent controllers for the leader and each follower. Figure 2 represents this control architecture. In this figure, index 0 represents the leader; therefore, *X* 0 *t* represents the current measured states for

method based on Automatic Train Control, Communication-Based Train Control, or the

Automatic Train Control system that tracks a given speed curve. This speed curve is

Because the role of the leader is to set the VCTSs' movement policy, we use the DP approach to precompute the reference behaviour of the leader. Then, the result of the DP establishes the general policy followed by the train's convoy. This general control policy can have different objective functions such as minimizing the energy consumption or maximizing the convoy speed. In this paper, we use the second policy with an optimal speed profile that finds the maximum velocity permitted by the speed limitations imposed

**Figure 3.** Maximum driving speed in compliance with speed limits.

obtained using a Dynamic Programming (DP) approach (Figure 3).

by the line operation, thereby satisfying the speed constraints at all times.

the leader; *v* 0 represents the leader's speed; *s* 0 is the leader's current position; and *U*0 represents the vector containing the driving/braking force for the leader. *X i t* represents the current measured states for the follower *i* in the platoon, which is virtually coupled to its previous train; *v i* and *s i* represent, respectively, the speed and current position of the follower *i*; and *Ui* represents the driving/braking force for the follower *i*. As seen in Section 2 for the linearized control model, *Ui o* is the value of *Ui* at the operating point. PTP 1 Train controller 0 Train controller 2 Train controller 1 PTP 2 *Xt 0 U 0 Xt 0 , Xt 1 , U0 1 , L Xt 1 , X 0 , L U 1 Xt 1 , Xt 2 , U0 2 , L Xt 2 , X 1 , L 1*

described in Section 3.5. For the followers, we design an MPC robust controller that ensures a safety and control efficiency while considering the parameter uncertainties.

and ) is used by the follower's controller for its calculations because the predicted states for the preceding train − 1 are calculated using the PTP embedded in each follower

With the proposed control strategy, the leader will track a given speed curve, and the followers will guarantee a safe minimum distance between trains and the string stability of the VCTS with parameter uncertainties. Consequently, as seen in Figure 2, the follower needs information about the preceding train to ensure its safety. This information (

*Actuators* **2022**, *11*, x FOR PEER REVIEW 6 of 18

based on the information received from the preceding train.

*L*

*s*

The end-front distance *dt i* between two consecutive trains at any time is calculated using (12), where *L* is the train length: *v 0 v 1 Xt 0 Xt , L , L*

$$d\_t^{\dot{i}} = s\_t^{\dot{i}-1} - s\_t^{\dot{i}} - L. \tag{12}$$

*L*

−1 ()

Furthermore, *X i*−1 represents the predicted states for the preceding train (*i* − 1) calculated by the follower's Preceding Train Predictor (PTP) module, which will be described in Section 3.5. For the followers, we design an MPC robust controller that ensures a safety and control efficiency while considering the parameter uncertainties. **Figure 2.** Control architecture. Data communication

With the proposed control strategy, the leader will track a given speed curve, and the followers will guarantee a safe minimum distance between trains and the string stability of the VCTS with parameter uncertainties. Consequently, as seen in Figure 2, the follower needs information about the preceding train to ensure its safety. This information (*X i*−1 *t* (*t*) and *L*) is used by the follower's controller for its calculations because the predicted states for the preceding train *i* − 1 are calculated using the PTP embedded in each follower *i* based on the information received from the preceding train. *3.2. Leader Controller* When using a decentralized control architecture, the leader can operate under any control and signalling method. Therefore, the leader can use any conventional control method based on Automatic Train Control, Communication-Based Train Control, or the ETCS. For simplicity, for the leading train, we use a driving mechanism based on an Automatic Train Control system that tracks a given speed curve. This speed curve is

#### *3.2. Leader Controller* obtained using a Dynamic Programming (DP) approach (Figure 3). Because the role of the leader is to set the VCTSs' movement policy, we use the DP

*d1-2*

*1 s*

*v 2*

*U 2*

*L*

*2*

When using a decentralized control architecture, the leader can operate under any control and signalling method. Therefore, the leader can use any conventional control method based on Automatic Train Control, Communication-Based Train Control, or the ETCS. approach to precompute the reference behaviour of the leader. Then, the result of the DP establishes the general policy followed by the train's convoy. This general control policy can have different objective functions such as minimizing the energy consumption or maximizing the convoy speed. In this paper, we use the second policy with an optimal

For simplicity, for the leading train, we use a driving mechanism based on an Automatic Train Control system that tracks a given speed curve. This speed curve is obtained using a Dynamic Programming (DP) approach (Figure 3). speed profile that finds the maximum velocity permitted by the speed limitations imposed by the line operation, thereby satisfying the speed constraints at all times.

**Figure 3.** Maximum driving speed in compliance with speed limits. **Figure 3.** Maximum driving speed in compliance with speed limits.

Because the role of the leader is to set the VCTSs' movement policy, we use the DP approach to precompute the reference behaviour of the leader. Then, the result of the DP establishes the general policy followed by the train's convoy. This general control policy can have different objective functions such as minimizing the energy consumption or maximizing the convoy speed. In this paper, we use the second policy with an optimal speed profile that finds the maximum velocity permitted by the speed limitations imposed by the line operation, thereby satisfying the speed constraints at all times.

Figure 3 shows the speed profile obtained for the leader, corresponding to the convoy policy of the maximum possible velocity allowed by the speed limitations established for the line's operation. In this figure, the curve labelled "Max speed" represents the maximum limitations set in the line design. In this figure, it can be seen where the stations are located and that the speed at these stations is limited to 12.5 m/s (45 km/h). The other speed limitations depend on the design of the railway line.

On the other hand, the curve labelled "DP speed" represents the maximum speed at which it is possible to run without exceeding at any time the speed allowed on the line and represents the result of the DP calculation. On this curve, the stopping point of the leading train at the end of each station can also be seen. A detailed explanation of the implementation of this DP approach can be found in [19].

### *3.3. Follower Controller*

For the controller design, a model predictive control (MPC) approach is used. The MPC optimizes over a finite time horizon but implements only the current time window of the finite horizon optimization problem solution.

For the MPC formulation, a prediction horizon - *t*, *t* + *Np* with origin at time *t* is considered. The notation *xt*+*k*|*t* represents the state vector at time *t* + *k* predicted at time *t*, obtained from the current state *xt*|*t* = *x*(*t*) ≡ *xt* , while *u*·|*t* = h *ut*|*t* , . . . , *ut*+*Np*−1|*t* i denotes the unknown input variables (or inputs) to be optimized.

The RMPC for virtual coupling is obtained by introducing two uncertain variables.

The first uncertain variable is *w a* (*m*/*s* 2 ) and is added to the acceleration equation of the dynamic model, as seen in (9).

The second uncertain variable *w p* (*m*) is considered in the virtual coupling specific constraints and the follower cost functions, and it is used to define the uncertainty in the positioning, as will be seen later.

As a result of introducing uncertain variables, the optimization problem for the follower *i* is formulated as the *min–max* problem given in (13):

$$\min\_{\mathbf{u}^i, \mathbf{u}^i, \mathbf{w}^{a,i}, \mathbf{w}^{p,i}} \max\_{\mathbf{w}^{p,i}\_{\cdot \mid t}} J^i \left( \mathbf{X}^i{}\_{k|t}, \mathbf{u}^i{}\_{k|t}, \mathbf{w}^{a,i}\_{k|t'}, \mathbf{w}^{p,i}\_{k|t} \right). \tag{13}$$

The controller's objective for each follower is to try to keep the trains running as close together as possible by minimising the difference in the distance between trains *dk*|*t* and a desired separation *ddes* while maintaining safe conditions. Therefore, we can formulate the following optimization problem:

> *X i*

$$J^i = \sum\_{k=1}^{k=N\_p} K\_D || \frac{d^i\_{k+1|t} + w^{p,i}\_{k+1|t} - d\_{\text{des}}}{d\_{\text{des}}} ||\_1 \tag{14}$$

$$d^i\_{k|t} = \mathfrak{s}^{i-1}\_{k|t} - s^i\_{k|t} - L$$

subject to:

$$X^{i}\_{\ k+1|t} = X^{i}\_{\ k|t} + t\_{s} \; \dot{X}^{i}\_{\ k|t} \tag{15}$$

*t*|*t* = *X i t* (16)

$$0 \le v\_{k|t}^i \le v\_{lim} \left( s\_{k|t}^i \right) \tag{17}$$

$$-j\_{\max} \le j\_{k|t}^i \le j\_{\max} \tag{18}$$

$$j\_{k|t}^i = \left(u\_{k+1|t}^i - u\_{k|t}^i\right) / M / \Delta t$$

$$-M a\_{br}^i \le u\_{k|t}^i \le M a\_{dr}^i \tag{19}$$

$$-P^i\_{br} \le \mathcal{Y}^i\_{1k|t} \le P^i\_{dr} \tag{20}$$

$$\forall k = t\_\prime \dots \iota, t + N\_p - 1$$

$$0 \le v\_{t+N\_p|t}^i \le v\_{DP}^{i^\*} \left( s\_{t+N\_p|t}^i \right) \tag{21}$$

$$d\_{\min} \le w\_{k|t}^{p,i} + d\_{k|t}^i \tag{22}$$

$$\forall k = t, \dots, t + N\_p - 1$$

$$d\_{\rm min} \le w\_{k|t}^{p,i} + d\_{t+N\_p|t}^i \tag{23}$$

$$d\_{\min} \le w\_{k|t}^{p,i} + d\_{t+Np|t}^i + \frac{\left(\overline{v}\_{t+Np|t}^{i-1}\right)^2}{2a^s} - Y\_{2k|t}^i \tag{24}$$

$$w\_{a,\min} \le w\_{k|t}^{a,i} \le w\_{a,\max} \tag{25}$$

$$w\_{p\_\prime \min} \le w\_{k|t}^{p,i} \le w\_{p\_\prime \max}.\tag{26}$$

Superscript *i* is related to the corresponding follower.

In the cost function *J i* (14), *KD* ≥ 0 is a dimensionless coefficient representing the weight that penalizes the deviation from the desired distance *ddes*, and *v i*−1 *t*+*Np*|*t* is the predicted speed of the preceding train, which is pre-calculated according to the PTP module that will be described later in Section 3.5.

The equations in (15) represent the train dynamics updates for the model obtained in (9)–(11) according to the *ts* integration time step. The initial state is set in (16).

Equation (17) corresponds to the velocity constraint *vlim* (Max speed in Figure 3).

Equation (18) represents the jerk constraint established for *u i* .

Equations (19) and (20) represent the input constraints, including the maximum driving/braking force and the maximum power on traction/braking, with *Y i j k*|*t* representing the estimated *j* component of vector *Y* at instant k calculated at time *t* according to (11).

Equation (21) represents a terminal constraint. The upper speed limit *v i* ∗ *DP*(*m*/*s*) is obtained by pre-calculating the maximum speed at which the train can run while respecting the speed limits set for the line (DP speed in Figure 3). This terminal constraint is used because the train must always respect the speed limit, and there are areas above the blue line in Figure 3 where if the train runs under these conditions, it may not respect the speed limits at future times.

Equation (22) represents the safety condition, stating that the distance between the trains must be greater than *dmin* at any time.

Equations (23) and (24) are the terminal constraints. These constraints are imposed to ensure that the controller is recursively feasible and safe (see Section 3.6). These equations guarantee that the follower can come to a complete stop without collision by applying a maximum service deceleration *a s* when the preceding train performs emergency braking with the maximum deceleration *a e* . Note that the terminal constraint (24) functions as an operational constraint, while the terminal constraint (23) works as a safety constraint. The safe train-to-train distance *dmin* is used to define a constraint to avoid a collision if the leader brakes. The distance *ddes* is used in the cost function as a target to maintain a safe desired driving distance greater than the minimum safe distance *dmin*. Finally, constraint (25) bounds the acceleration uncertainty *w a* and (26) bounds the second uncertain variable *w p* (m), which addresses the distance error from the preceding train. Both constraints are handled through the traditional method of robust optimization based on the duality theory (for more details, see [20]).

The first uncertain variable *w a* (m/s2 ) is added to the dynamic model's acceleration equation, as seen in (2). By means of this uncertainty, we can consider different phenomena. For instance, using this variable, an additional resistance due to the start resistance can be taken into account, thus solving the lack of accuracy of the Davis formula at speeds below 3 m/s. In addition, this variable can counteract the additional resistance introduced

by the tunnel factor and the longitudinal component of lateral wind's impact. Moreover, this variable can avoid the risky situations produced by perturbation, such as a loss of adhesion during braking. This variable also enables the correction of the modelling errors produced by incorrect estimations of the Davis coefficients. The specific values that limit this uncertain variable will be discussed in Section 4.1.

Finally, the second uncertain variable *w p* (m) addresses the positioning error that can influence the estimation of the follower's distance from the preceding train. The specific values that limit this uncertain variable will be discussed in Section 4.2.

#### *3.4. Control Loop*

For each train, the resulting optimal states and inputs of (13)–(26) are denoted as follows:

$$\begin{array}{ccccc} \mathbf{X}\_{t}^{\mathrm{i}^\*} = \left( \begin{array}{cccc} \mathbf{X}\_{t|t}^{\mathrm{i}^\*} & \mathbf{X}\_{t+1|t}^{\mathrm{i}^\*} & \dots & \mathbf{X}\_{t+N\_p|t}^{\mathrm{i}^\*} \\ \mathbf{u}\_{t}^{\mathrm{i}^\*} = \left( \begin{array}{cccc} u\_{t|t}^{\mathrm{i}^\*} & u\_{t+1|t}^{\mathrm{i}^\*} & \dots & u\_{t+N\_p-1|t}^{\mathrm{i}^\*} \end{array} \right)^T . \end{array} \tag{27}$$

Since constraints (16), (17), and (21) are space-dependent, to avoid numerical problems and reduce the computational burden, they can be estimated a priori by *s i t*+*k*|*t* :

$$
\overline{X}^i\_{\;k+1\vert t} = \overline{X}^i\_{\;k\vert t} + t\_s \stackrel{\cdot}{\overline{X}}^i\_{\;k\vert t} \tag{28}
$$

$$\forall k = t, \dots, t + N\_p - 1$$

$$
\overline{u}\_{k|t}^i = u\_{k|t-1}^{i^\*} \tag{29}
$$

$$\forall k = t, \dots, t + N\_p - 2$$

$$\overline{u}\_{t + N\_p - 1 \mid t}^i = u\_{t + N\_p - 2 \mid t - 1}^{i^\*} \tag{30}$$

$$\overline{X}^i\_{t|t} = X^i\_{t|t} \tag{31}$$

where *u i* ∗ *k*|*t*−1 is the input predicted from (27) at the previous time step *t* − 1 with *u i* ∗ *k*|*t*−1 = 0 ∀*k* = 0, . . . , *Np* − 2.

To close the loop, the first input is applied to system (9)–(11) during the time interval [*t*, *t* + 1 ):

$$
\mu^i\_t = \mu^{i^\*}\_{t|t}.\tag{32}
$$

At the next time step *t* + 1, a new optimal problem in the form of (13)–(26) based on a new state measurement will be solved over a shifted horizon.

### *3.5. Preceding Train Predictor*

To establish the control loop, it is necessary to estimate a prediction of the preceding train's movement (train *i* − 1). This task is performed by the module, called the Preceding Train Predictor (PTP), which is described in this section.

Since both vehicles are equipped with a communication module, and intervehicular communication is available, the prediction task is executed by the follower based on the information received from the preceding train. For the prediction, the preceding train *i* − 1 communicates to the follower its position *s i*−1 *t*|*t* , speed *v i*−1 *t*|*t* , and the driving/braking force *F i*−1 ∗ *t*|*t*−1 . In this way, the distance between trains *d i k*|*t* is determined by the difference in the positions of the two trains. On the other hand, and as a redundancy, the follower has a sensor that can measure the distance to the front train. Therefore, for the estimation of the position of the preceding train, the worst case is considered, which corresponds to the lower of the two calculated values above.

Then, the MPC assumes that the preceding train's information can be predicted over a short horizon *Np*. The corresponding expressions for the prediction of *s i*−1 *k*|*t* and *v i*−1 *k*|*t* are

$$\overline{X}\_{t|t}^{i-1} = \begin{pmatrix} s\_{t|t}^{i-1} & v\_{t|t}^{i-1} & F\_{t|t-1}^{i-1^\*} \end{pmatrix} \tag{33}$$

$$\overline{X}^{i-1}{\overline{X}^{i-1}}\_{k+1|t} = \overline{X}^{i-1}{\overline{X}^{i-1}}\_{k|t} + t\_s \stackrel{\cdot}{\overline{X}}^{i-1}{\overline{X}^{i-1}}\_{k|t} \tag{34}$$

$$\forall k = t+1, \dots, t+N\_p-1,$$

where

$$
\overline{u}\_{k|t}^{i-1} = u\_{k|t-1}^{i-1^\*} \tag{35}
$$

$$\forall k = t, \dots, t + N\_p - 2$$

$$\overline{u}\_{t + N\_p - 1 \mid t}^{i - 1} = u\_{t + N\_p - 2 \mid t - 1}^{i - 1^\*} \tag{36}$$

Here, *s i*−1 *t*|*t* and *v i*−1 *t*|*t* are the measured position and velocity of the leader at *t*, respectively; *F i*−1 ∗ *t*|*t*−1 is the force at time *t* predicted from (27) at the previous time step *t* − 1 initialized with *F i*−1 ∗ 0|*t*−1 = 0; and *u i*−1 ∗ *k*|*t*−1 is the input predicted from (27) at the previous time step *t* − 1, with *u i*−1 ∗ *k*|*t*−1 = 0 ∀*k* = 0, . . . , *Np* − 2.

In (34), . *X i k*|*t* is calculated from the expressions of (9)–(11), particularly (10), and the disturbance *w a* is considered to be the lower limit *wa*,*min*, which is the most unfavourable situation corresponding to the highest risk of a collision between the two consecutive trains.

#### *3.6. Min-Max Recursive Feasibility*

To ensure the recursive feasibility of the RMPC, we define set C as the robust control invariant set. According to [21], set C is said to be a robust control invariant set for the system (9)–(11) subject to constraints (13)–(25) if

$$X(t) \in \mathbb{C} \implies \exists [u] \in \mathbb{U} \text{ such that } system \ (4) \in \mathbb{C}, \ \forall w(t) \in \mathbb{W} \ \forall t \ge 0 \tag{37}$$

where U is the set of the valid decision variable values given by (18)–(20), and W is the set that bounds the uncertain variables, i.e., the set given by (25) and (26).

As in [22], set C can be introduced by the constraints (21)–(23), which define a robust control invariant set because (21) ensures a proper decision of the train's controller according to the track profile, while (22)–(24) ensure a safe reaction to the emergency braking of the preceding train. In other words, these constraints ensure that there exists an output value for the control variables that satisfies (37).

#### **4. Simulations**

This section compares the results obtained with the nominal controller developed in [19] with those obtained with the robust controller developed in this paper. To test the robust MPC for the VCTS, two simulation cases were used for a metro line. The first simulation case analyses the influence of an uncertainty in acceleration *w a* , while the second one studies the uncertainty in the train's positioning *w p* .

The parameters considered for the rolling stock can be found in Table A1 of Appendix A. The profile of the selected line can be found in Figure 1, and its corresponding speed profile can be found in Figure 3.

The simulations were performed by comparing the results between an NMPC (without uncertainties) and an RMPC, where uncertainties are considered.

Simulations were performed in MATLAB and Yalmip [23] on a computer with an i7-1165G7 processor at 2.8 GHz and with 8 GB of RAM. For the different simulations, we determined that the average CPU computation time of the solver was 0.17 s for the RMPC. Moreover, the total CPU computation time was found to be 302 s, which is lower than

the total simulated time (345 s). Therefore, the RMPC presents promising results for a hardware implementation that needs to be solved in real time.

#### *4.1. Simulation 1: Uncertainty in Acceleration*

At the end of Section 3.3, we introduced what the first uncertain variable *w a* (m/s2 ) can represent. Now, from a practical point of view in the train operation, we bound the values of the acceleration uncertainty (25) according to Section 3.3.

The start resistance of a train on a track with a null slope is 50 N/t, i.e., 0.05 m/s2 . Since this resistance is applied against the movement (negative), −0.05 m/s2 is established as the lower limit *wa*,*min*. Even though this value is lower because the Davis formula (2) is still applied in our controller, it allows for the simultaneous consideration of additional small resistances. For instance, the maximum errors in the estimation of the Davis coefficients for passenger trains are small (below 0.01 m/s2 for this train) and would be equally compensated by the controller thanks to this setting. Therefore, this value corrects the modelling errors introduced by the Davis formula: both above 3 m/s (Davis coefficients) and lower speeds (start resistance and Davis coefficients).

Additional resistances under 0.05 m/s2 will also be correctly compensated at speeds above 3 m/s, which is the case for an external perturbation, such as the resistance introduced by either the tunnel factor or the longitudinal component of a lateral wind's impact.

Conversely, *wa*,*max* is established to counteract an extreme loss of adhesion of 10% during braking. If a friction factor of 0.1 is considered, this loss will produce a reduction in the deceleration of 0.1 m/s2 . This reduction works as an additional traction force during braking. Therefore, 0.15 m/s2 is established as the upper bound *wa*,*max* as a conservative approximation.

Among all the external perturbations, only the extreme perturbation of a 10% loss of adhesion will be studied in the simulation cases of this paper, but these results can be extrapolated to any of the acceleration uncertainty effects mentioned above.

In Simulation 1, the behaviour of the uncertain variable *w a* has been assessed and corresponds to a disturbance involving a loss of adhesion at the entrance of the third station. This situation is extremely unfavourable since the disturbance was applied only to the follower, and the train in front was considered to continue braking without any loss of adhesion.

Figures 4–6 represent the behaviour of the convoy when the follower experiences a 10% loss of adhesion during braking.

The upper plot in Figure 4 represents the time–velocity curve and shows the speed of the leader, which is the same for both the nominal and robust controllers. Moreover, in acceleration (i.e., the increasing parts of the curves), the leader has a slightly higher velocity than the followers, thus increasing the distance between the trains. However, in braking (i.e., the decreasing part of the curve), when the velocity of the leader decreases, the follower becomes closer as it has a higher speed, always maintaining the safe distance imposed by the constraints. It is also possible to observe the stop at the stations.

In Figure 4, it can also be seen how the integrity of the convoy is maintained, as the two trains stop at the station at almost the same time. This effect can be seen in the time–velocity plot, where both trains in the convoy stop at the stations at almost the same time. The plot below shows the distance between the leader and the follower. Here, the most critical moment occurs at the entrance of the stations, where the leader stops, and the follower must also stop maintaining a safe distance.

The curves labelled "no adhesion loss" include the distance between the leader and the follower for the NMPC and RMPC when there is no disturbance, i.e., no adhesion loss. These curves show a very similar behaviour in both cases.

If the behaviour of the robust controller is compared with that of the nominal controller without disturbances, it can always be observed that the distance between two consecutive trains is always larger in the case of the robust controller. This result is logical since the controller contemplates the presence of uncertainty in its own design. In this case, the minimum distance between the trains is obtained at the station entrance, when the leader has stopped. It can also be seen that in the NMPC, this distance is 5 m, while the RMPC is more conservative and maintains a distance of approximately 8 m.

*Actuators* **2022**, *11*, x FOR PEER REVIEW 12 of 18

**Figure 4.** Results of simulation 1: 10% adhesion loss. Variables versus time.

**Figure 4.** Results of simulation 1: 10% adhesion loss. Variables versus time.

controller without disturbances, it can always be observed that the distance between two

consecutive trains is always larger in the case of the robust controller. This result is logical

since the controller contemplates the presence of uncertainty in its own design. In this

case, the minimum distance between the trains is obtained at the station entrance, when

the leader has stopped. It can also be seen that in the NMPC, this distance is 5 m, while

a 10% loss of adhesion is considered in the entrance area of the third station. In this case,

while the RMPC can maintain a minimum distance of 8 m between the trains at the

entrance to the station, a distance of 0.5 m is obtained with the NMPC. This situation

seriously compromises safety because a collision is about to occur, and the safety

However, in the case of a loss of adhesion, differences are observed. In the simulation,

the RMPC is more conservative and maintains a distance of approximately 8 m.

conditions for which the nominal control was designed are not respected.

**Figure 5.** Results of simulation 1: 10% adhesion loss. Variables versus distance travelled. **Figure 5.** Results of simulation 1: 10% adhesion loss. Variables versus distance travelled.

Figure 5 represents the speed and distance between the trains versus the distance travelled. How and when the trains stop at the station is accurately visualized in this figure. It can also be seen that the front train stops at the end of the station, and the rear train stops somewhat further back at the desired distance. Additionally, Figure 6 represents the driving/force obtained as a controlled variable. However, in the case of a loss of adhesion, differences are observed. In the simulation, a 10% loss of adhesion is considered in the entrance area of the third station. In this case, while the RMPC can maintain a minimum distance of 8 m between the trains at the entrance to the station, a distance of 0.5 m is obtained with the NMPC. This situation seriously compromises safety because a collision is about to occur, and the safety conditions for which the nominal control was designed are not respected.

Figure 5 represents the speed and distance between the trains versus the distance travelled. How and when the trains stop at the station is accurately visualized in this figure. It can also be seen that the front train stops at the end of the station, and the rear train stops somewhat further back at the desired distance. Additionally, Figure 6 represents the driving/force obtained as a controlled variable.

**Figure 6.** Results of simulation 1: 10% adhesion loss. Force vs. time. **Figure 6.** Results of simulation 1: 10% adhesion loss. Force vs. time.

#### *4.2. Simulation 2: Uncertainty in Positioning Error 4.2. Simulation 2: Uncertainty in Positioning Error*

At the end of Section 3.3, we introduced the second uncertain variable (m), which represents the positioning error that can influence the estimation of the follower's distance from the preceding train. Now, from a practical point of view in the train operation, we bound the values of the positional uncertainty (26) according to Section 3.3. At the end of Section 3.3, we introduced the second uncertain variable *w p* (m), which represents the positioning error that can influence the estimation of the follower's distance from the preceding train. Now, from a practical point of view in the train operation, we bound the values of the positional uncertainty (26) according to Section 3.3.

The positioning error is produced by two components. To begin with, ATO (Automatic Train Operation) controllers have a specified maximum positioning error of ±0.5 m [24]. Additionally, a velocity sensor usually has a maximum error of 1% of the real speed. If a metro line with a 110 km/h maximum speed is considered, the velocity The positioning error is produced by two components. To begin with, ATO (Automatic Train Operation) controllers have a specified maximum positioning error *ep* of ±0.5 m [24]. Additionally, a velocity sensor usually has a maximum error *ev* of 1% of the real speed. If a metro line with a 110 km/h maximum speed is considered, the velocity error will be ±0.3 m/s.

error will be ±0.3 m/s. Hence, the maximum distance error produced by the follower predictor appears at the end of the prediction horizon (38), where the error is more critical in the constraints. In (38), this result is multiplied by two to consider an extreme situation in which the preceding train is nearer, and the follower is farther. Since this distance error implies that the preceding train might be nearer, this value is considered in the lower bound. Therefore, −3.5 m is estimated as the lower bound of , as a conservative approximation of (38). Nevertheless, 0 m is set as the upper bound , in order to Hence, the maximum distance error produced by the follower predictor appears at the end of the prediction horizon (38), where the error is more critical in the constraints. In (38), this result is multiplied by two to consider an extreme situation in which the preceding train is nearer, and the follower is farther. Since this distance error implies that the preceding train might be nearer, this value is considered in the lower bound. Therefore, −3.5 m is estimated as the lower bound of *wp*, *min* as a conservative approximation of (38). Nevertheless, 0 m is set as the upper bound *wp*, *max* in order to avoid a distance reduction above virtual coupling conditions, which could end in a collision between the trains:

avoid a distance reduction above virtual coupling conditions, which could end in a

$$d\_{error} = \mathcal{D}\left(e\_p + t\_s N\_p e\_v\right) = 3.4 \text{ m.} \tag{38}$$

 = 2( + ) = 3.4 . (38) This simulation focuses on testing the uncertain variable , which is applied to test the performance of the system when there is noise in the positioning information used by the follower. In this simulation, two sinusoidal signals are introduced in the positioning and speed information of the leading train ( and ). These signals have a period of 90 s and amplitudes of 0.8 m and 0.6 m/s, respectively. In addition, a random noise is superimposed on the previous wave. This noise is uniformly distributed and has an This simulation focuses on testing the uncertain variable *w p* , which is applied to test the performance of the system when there is noise in the positioning information used by the follower. In this simulation, two sinusoidal signals are introduced in the positioning and speed information of the leading train (*s l* and *v l* ). These signals have a period of 90 s and amplitudes of 0.8 m and 0.6 m/s, respectively. In addition, a random noise is superimposed on the previous wave. This noise is uniformly distributed and has an amplitude of 0.001 in each signal (m and m/s, respectively).

amplitude of 0.001 in each signal (m and m/s, respectively). The simulation was performed in both the nominal MPC and robust MPC. The results presented in Figure 5 under the label "errors" correspond to the simulation with disturbances, while the results without the label "errors" correspond to the simulation The simulation was performed in both the nominal MPC and robust MPC. The results presented in Figure 5 under the label "errors" correspond to the simulation with disturbances, while the results without the label "errors" correspond to the simulation without disturbances.

without disturbances.

Figure 7 shows that an error in the estimation of the position of the preceding train, e.g., due to a failure in the odometry or noise in the position information of the preceding train, can cause a collision with the NMPC (green dashed line in Figure 7). This collision can occur at the end of braking from high speeds, i.e., near stations. Figure 7 shows that an error in the estimation of the position of the preceding train, e.g., due to a failure in the odometry or noise in the position information of the preceding train, can cause a collision with the NMPC (green dashed line in Figure 7). This collision can occur at the end of braking from high speeds, i.e., near stations.

*Actuators* **2022**, *11*, x FOR PEER REVIEW 15 of 18

**Figure 7.** Simulation results 2: Errors in the preceding train information. **Figure 7.** Simulation results 2: Errors in the preceding train information.

The results of the simulations show that in the case of the NMPC, the minimum distance reached with this disturbance was found to be 0.8 m near the 300 s instant of the The results of the simulations show that in the case of the NMPC, the minimum distance reached with this disturbance was found to be 0.8 m near the 300 s instant of the simulation time.

simulation time. In contrast, in the case of the RMPC, these risk situations are avoided (the solid magenta line in Figure 7). Additionally, Figure 8 represents the driving/force obtained as a

controlled variable. In this simulation, the RMPC was found to achieve a minimum safe distance of 5.1 m, whereas the NMPC presented 0.8 m in the same situation. a controlled variable. In this simulation, the RMPC was found to achieve a minimum safe distance of 5.1 m, whereas the NMPC presented 0.8 m in the same situation.

In contrast, in the case of the RMPC, these risk situations are avoided (the solid magenta line in Figure 7). Additionally, Figure 8 represents the driving/force obtained as

*Actuators* **2022**, *11*, x FOR PEER REVIEW 16 of 18

**Figure 8.** Results of simulation 2: errors in the preceding train information. Force vs. time. **Figure 8.** Results of simulation 2: errors in the preceding train information. Force vs. time.

The most important result of this simulation is that the nominal MPC control system tends to lead to risky situations when the environment in which it operates is not a nominal one. This statement is corroborated in Figure 7, which illustrates that, at the station entrance, trains approach up to an absolute distance of 0.8 m, which is an The most important result of this simulation is that the nominal MPC control system tends to lead to risky situations when the environment in which it operates is not a nominal one. This statement is corroborated in Figure 7, which illustrates that, at the station entrance, trains approach up to an absolute distance of 0.8 m, which is an undesirably risky situation.

undesirably risky situation. However, the robust controller results in a minimum distance of 5.1 m compared to 7.5 m for an operation without disturbances. Nevertheless, despite the unfavourable situation, a distance greater than the specified odometry error (3.5 m) is maintained at 48.6%, thereby ensuring a safe traffic flow. Furthermore, with respect to the NMPC, the minimum distance achieved by the RMPC in this unfavourable situation is almost six times greater. In other words, the NMPC reduces the distance in this situation by 83% due However, the robust controller results in a minimum distance of 5.1 m compared to 7.5 m for an operation without disturbances. Nevertheless, despite the unfavourable situation, a distance greater than the specified odometry error (3.5 m) is maintained at 48.6%, thereby ensuring a safe traffic flow. Furthermore, with respect to the NMPC, the minimum distance achieved by the RMPC in this unfavourable situation is almost six times greater. In other words, the NMPC reduces the distance in this situation by 83% due to its lack of robustness, which leads to risky situations.

#### to its lack of robustness, which leads to risky situations. **5. Conclusions**

**5. Conclusions** In this paper, a decentralized RMPC for the VCTS in railroads was developed. Two types of uncertainties were included in the controller. We compared the proposed In this paper, a decentralized RMPC for the VCTS in railroads was developed. Two types of uncertainties were included in the controller. We compared the proposed RMPC controller with an NMPC in different simulations that consider the positioning errors and a loss of adhesion during braking.

RMPC controller with an NMPC in different simulations that consider the positioning errors and a loss of adhesion during braking. The result is a more conservative controller than a nominal controller, which is nevertheless capable of guaranteeing collision safety in situations which the nominal The result is a more conservative controller than a nominal controller, which is nevertheless capable of guaranteeing collision safety in situations which the nominal controller is not able to overcome, such as under a 10% loss of adhesion during braking and in the case of errors in the position information of the front train.

controller is not able to overcome, such as under a 10% loss of adhesion during braking and in the case of errors in the position information of the front train. However, it is possible to improve the estimation of the virtual coupling predictor. In future research, fewer conservative solutions could be proposed based on the leader's However, it is possible to improve the estimation of the virtual coupling predictor. In future research, fewer conservative solutions could be proposed based on the leader's braking estimates calculated from the last known information instead of using only the last known information directly in the controller.

braking estimates calculated from the last known information instead of using only the last known information directly in the controller. On the other hand, the decentralized control architecture is acceptable as an advanced railway signalling system that improves the current train operations but not from the perspective of the VCTSs' concept, since in the proposed solution, the leader sets the VCTSs' movement policy but is not affected by the behaviour of the followers. Therefore, in the

near future, we will use these results to design a robust centralized MPC for the VCTSs' strategy where all the elements of the VCTS can interact.

**Author Contributions:** Conceptualization, J.F., M.A.V.-S. and J.d.D.S.; methodology, J.F. and M.A.V.-S.; software, J.F. and M.A.V.-S.; investigation, J.F., M.A.V.-S. and J.d.D.S.; writing—original draft preparation, J.F. and M.A.V.-S.; writing—review and editing, J.F. and M.A.V.-S.; supervision, J.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the H2020 ECSEL EU Project Intelligent Secure Trustable Things—InSecTT (www.insectt.eu, (accessed on 7 December 2022)) receiving funding from the ECSEL Joint Undertaking (JU) under grant agreement No 876038, and in part by the Spanish Science and Innovation Ministry—State Research Agency under Grant PCI2020-112126/AEI/10.13039/<PHONE_NUMBER>33. The JU receives support from the European Union's Horizon 2020 research and innovation programme and Austria, Sweden, Spain, Italy, France, Portugal, Ireland, Finland, Slovenia, Poland, Netherlands, Turkey. The document reflects only the author's view and Commission is not responsible for any use that may be made of the information it contains.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

### **Appendix A**

**Table A1.** Parameters used in the simulation.

#### **References**

## *Article* **Time-Optimal Current Control of Synchronous Motor Drives**

**Václav Šmídl \* , Antonín Glac and Zdenˇek Peroutka**

Research and Innovation Centre for Electrical Engineering (RICE), University of West Bohemia, 30100 Plze ˇn, Czech Republic

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** We are concerned with the problem of fast and accurate tracking of currents in the general synchronous drive. The problem becomes complicated with decreasing available voltage, which is common in high-speed and field weakening regimes. The existing time-optimal controllers rely on a simplified model, ignoring stator resistance and differences in inductances. We derive a solution for the general model considering all parameters and show how the parameters affect the current trajectory. One simplifying assumption had to be made, but we show in simulation that it has a negligible impact on accuracy. The simplification allowed for the design of a feed-forward controller that has a low computational cost and can be easily implemented in realtime. We provide experimental validation of the controller on the developed IPMSM drive prototype of the rated power of 4.5 kW using conventional industrial DSP. The controller is compared to conventional PI and deadbeat solutions, demonstrating that the time-optimal controller can reach the required setpoint four times faster than the competitors at the field weakening regime of the drive. The proposed feed-forward control can be seen as a universal building block that can be combined with existing feedback controllers and observers and thus incorporated into existing control solutions.

**Keywords:** interior permanent magnet synchronous motor (IPMSM); maximum torque per ampere (MTPA); maximum torque per current (MTPC); dead beat control, predictive control

### **1. Introduction**

Control of synchronous motor drives is traditionally decomposed into nested loops where the current loop is the fastest of them. A prominent example is the conventional fieldoriented control of using a cascade of PID controllers [1]. Since the time constants of the current loops are relatively short, the optimal operating points of the drive are usually given as inputs to the current control loop [2]. The optimal current setpoints are computed for steady-state operations, using well-known results such as the maximum torque per ampere MTPA curve [3]. Even in the calculation of the MTPA, the stator resistance of the drive is often neglected to achieve simpler solutions. This restriction has been recently relaxed by designing optimal steady-state references for the full electric model of a synchronous motor drive [4].

Optimization of the steady state operation is typically considered for permanent operation of the drive, with application to pumps or ventilation, hence the most common concern of optimization is efficiency [4,5]. However, this optimization is incomplete for highly dynamic operations such as robotic actuators or manipulators, where high precision of set-point tracking is required. In such applications, the set-point trajectory is often known in advance, allowing optimization of the feed-forward part of the controller. In this paper, we focus on feed-forward optimization of the current loop. All other aspects of the control, such as set-point design [4], feed-back part of the control [6], and disturbance rejection [7,8] will be briefly commented at the end of the paper.

Computation of the optimal feed-forward trajectories for transients can be achieved numerically, using model predictive control [9]. However, even many predictive control solutions often use steady-state solutions [10,11] or decomposition into speed and current

**Citation:** Šmídl, V.; Glac, A.; Peroutka, Z. Time-Optimal Current Control of Synchronous Motor Drives. *Actuators* **2023**, *12*, 15. https://doi.org/10.3390/ act12010015

Academic Editor: Shuxiang Dong

Received: 30 November 2022 Revised: 16 December 2022 Accepted: 21 December 2022 Published: 29 December 2022

**Copyright:** © 2022 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 (https:// creativecommons.org/licenses/by/ 4.0/).

loop [12,13]. While MPC can also be utilized on a higher level with beneficial results [14,15], its application to demanding current control transients is not successful. The reason is the hard constraint on the input voltage. Due to this constraint, a very long prediction horizon is required to obtain the optimal current trajectory with a completely different trajectory for shorter horizons [16]. The results of long-horizon optimization are counter-intuitive, since the optimal current profile is increasing the tracking error at the beginning of the transient. The reason is that it will reach the set-point in a shorter time than if it tried to decrease the error. This solution can be obtained numerically with very high computational cost. While it is possible to transform the numerical solution to explicit model predictive control [12,17], any change in motor parameters requires recomputing of the expensive solution.

In this paper, we address this problem using time-optimal control (TOC) methods [18]. Minimizing the time to convergence can be treated using general-purpose methods [19]; however, some problems allow for an analytical solution. Specifically, we will be formalizing the current control problem in continuous time, where it can be addressed using Pontryagin's principle of maxima. This technique has been applied to the current control of PMSM drive problem, which has been studied in [20,21] using highly simplifying assumptions, such as neglected resistance. The solution of the full model problem is conceptually known [22,23]; however, its evaluation for the full model is non-trivial. It can be found, e.g., by dynamic programming [24], which is computationally demanding. An extension considering resistance has been presented in [25]; however, without considering different inductances in the direct and quadrature axis. In essence, the TOC serves to generate current setpoints for low level controllers such as deadbeat [2,26].

The contributions of our paper are as follows:

The paper is organized as follows. The review of previous approaches is summarized in Section 2. The general formula of time-optimal control is derived in Section 3. A sensitivity study of the proposed control to parameters of the drive is performed in Section 4 using simulation. The experimental comparison with a deadbeat controller and the PI controller is reported in Section 5.

#### **2. Review of Existing Approaches**

#### *2.1. Mathematical Model of Synchronous Motor Drive*

Mathematical model of a synchronous motor drive in the *dq* reference frame linked to a rotor flux linkage vector is generally described by stator flux dynamics [4]

$$
\psi\_s(t) = -R\_s i\_s(t) - \omega(t) \mathbf{J} \psi\_s(t) + \mathfrak{u}\_s(t), \tag{1}
$$

where *ψs* denotes the vector of stator flux in the *dq* reference frame, *ω*(*t*) is the electrical rotor speed, *Rs* is the stator winding resistance, *J* is the rotation matrix *J* = [0, −1; 1, 0] and *us*(*t*) is the stator voltage vector in the *dq* reference frame. We will assume that the flux can be (at least locally) approximated by linear formula

$$
\psi\_s(t) = L\_s i\_s(t) + \psi\_{\rm pm} \tag{2}
$$

where we assume a diagonal matrix of inductances *Ls* = [*Ls*,*d* , 0; 0, *Ls*,*q*] and flux excited by the permanent magnets on the rotor *ψ*pm = [*ψ*pm,*d* ; *ψ*pm,*q*]. In this work, we consider *ψ*pm to be a general constant vector, which can be specialized for particular drives by setting its elements to zero [4]. For example, *ψ*pm,*q* = 0 for PMSM and PM-enhanced RSM, *ψ*pm,*d* = 0 for PM-assisted RSM, and *ψ*pm = [0; 0] for RSM drive. We will derive all equations for the general synchronous motor model.

Since we are concerned with drives that have a much faster current response than mechanical dynamics, we will consider the rotor speed *ω* to be constant during the transient. The state variable is thus the vector of fluxes *x* = *ψ* = [*ψs*,*d* ; *ψs*,*q*], with the control variable *u* = *us* and with dynamics

$$
\dot{\mathfrak{x}}(t) = \mathbf{A}\mathfrak{x}(t) + \mathfrak{u}(t) + \mathfrak{q},\tag{3}
$$

where

$$\mathbf{A} = -\mathbf{R}\_{\mathbf{s}} \mathbf{L}\_{\mathbf{s}}^{-1} - \omega \mathbf{J}\_{\prime} \qquad \qquad \qquad \qquad \boldsymbol{\eta} = \mathbf{R}\_{\mathbf{s}} \mathbf{L}\_{\mathbf{s}}^{-1} \boldsymbol{\Psi}\_{\text{pm}}.$$

The trajectory of the system is thus

$$\mathbf{x}(t) = e^{tA} \left( \mathbf{x}(0) + \int\_0^t e^{-sA} (\mathbf{u}(s) + \mathbf{q}) ds \right). \tag{4}$$

#### *2.2. Time-Optimal Control (TOC)*

Time-optimal control is a special case of optimal control where the aim is to reach the desired state *x*des from the original state *x*0 in the minimum possible time *τ*. The wellknown Pontryagin principle of maxima [18] states necessary conditions for the optimal trajectory *x* ∗ and optimal control *u* ∗ via two dynamic equations

$$\begin{aligned} \dot{\mathfrak{x}}^\*(t) &= \frac{\partial H}{\partial \mathfrak{p}}(t) = \mathbf{A}\mathfrak{x}^\*(t) + \mathfrak{u}^\*(t) + \mathfrak{q}(t),\\ \dot{\mathfrak{p}}(t) &= -\frac{\partial H}{\partial \mathfrak{x}}(t) = -\mathbf{A}^T \mathfrak{p}(t), \end{aligned}$$

and two boundary conditions *x*(0) = *x*0 and *x*(*τ*) = *x*des. Here, *H*(*p*(*t*), *x*(*t*), *u*(*t*)) = 1 + *p*(*t*) *T* (*Ax*(*t*) + *u*(*t*) + *q*) is the Hamiltonian of the system, and *p* is known as the costate (adjoint) variable satisfying

$$\mathfrak{p}(t) = e^{-t\mathbf{A}^T}\mathfrak{p}\_0. \tag{5}$$

The optimal control *u*(*t*) has to satisfy for almost every *t* ∈ [0, *τ*] the maximum principle

$$H(\mathfrak{p}(t), \mathfrak{x}^\*(t), \mathfrak{u}^\*(t)) = \max(H(\mathfrak{p}(t), \mathfrak{x}^\*(t), \mathfrak{u}(t)))$$

$$\text{s.t. } ||\mathfrak{u}(t)|| \le \overline{\mathcal{U}}.\tag{6}$$

Under the circular voltage constraint (6), the optimal control has the form

$$\mathfrak{u}(t) = \overline{\mathcal{U}} \frac{\mathfrak{p}(t)}{||\mathfrak{p}(t)||}. \tag{7}$$

Substituting (7) into (4), we obtain

$$\mathbf{x}(t) = e^{tA} \left( \mathbf{x}\_0 + \int\_0^t \left[ \overline{\mathbf{U}} \frac{e^{-s\mathbf{A}} e^{-s\mathbf{A}^T} \mathbf{p}\_0}{\sqrt{\mathbf{p}\_0 e^{-s\mathbf{A}} e^{-s\mathbf{A}^T} \mathbf{p}\_0}} + e^{-s\mathbf{A}} q \right] ds \right). \tag{8}$$

Solving (8) together with constraint *x*(*τ*) = *x*des for *τ* and *p*0 yields the required solution. Note that we may normalize the costate by setting ||*p*0|| = 1. The complexity of the solution depends highly on the model matrix *A*.

#### *2.3. TOC with Neglected Stator Resistance*

To our best knowledge, all previous approaches such as [20,21,24], were designed with neglected stator resistance *Rs* ≈ 0. It is understandable because this assumption is very convenient and allows us to greatly simplify the problem, since *A* = −*ωJ* and *q* = 0. Then,

$$e^{tA} = \begin{bmatrix} \cos\omega t & \sin\omega t \\ -\sin\omega t & \cos\omega t \end{bmatrix} \text{ , } e^{-sA}e^{-sA^T} = I\_A$$

Using assumption ||*p*0|| = 1, Equation (8) becomes

$$\mathbf{x}(t) = e^{t\mathbf{A}} \left( \mathbf{x}\_0 + \overline{\mathbf{U}} \mathbf{p}\_0 t \right). \tag{9}$$

Substituting *t* = *τ* and using the terminal condition *x*(*τ*) = *x*des in (9) yields *p*0 = 1 *Uτ e* −*τAx*des − *x*0 . Since *p*0 has to satisfy the normalization condition ||*p*0|| = 1, we obtain an implicit equation:

$$||e^{-\tau A} \mathfrak{x}\_{\text{des}} - \mathfrak{x}\_0||^2 = (\overline{\mathcal{U}}\tau)^2,\tag{10}$$

where the only free variable is *τ*. Finding the solution of (10) can be done using, e.g., the bisection method.

#### **3. Time-Optimal Flux Control Considering Stator Resistance**

Since the assumption of negligible resistance is not valid in many cases, we now use the study solution of the time-optimal control (7) for more general cases. First, we will analyze the case of equal inductances *Ls*,*d* ≈ *Ls*,*q*, which is a reasonable assumption for surface mounted PMSM. The more general case is studied subsequently.

### *3.1. TOC with Equal Stator Inductances*

Under the assumption of equal stator inductances in the *d* and *q* axes, *Ls* = *Ls*,*d* = *Ls*,*q* the system (3) has the form

$$A = -\rho I - \omega J, \ \mathfrak{q} = \rho \mathfrak{y}\_{\text{pm}}.$$

where *ρ* = *Rs*/*Ls* , yielding

$$e^{tA} = e^{-\rho t} \begin{bmatrix} \cos \omega t & -\sin \omega t \\ \sin \omega t & \cos \omega t \end{bmatrix} \; \; e^{-sA} e^{-sA^T} = e^{2\rho t} I \,\tag{11}$$

then, Equation (8) becomes

$$\mathbf{x}(t) = e^{tA} \left( \mathbf{x}\_0 + \int\_0^t \left[ e^{\rho s} \overline{\mathbf{U}} \mathbf{p}\_0 + e^{-sA} \mathbf{q} \right] ds \right) \tag{12}$$

$$\mathbf{r} = e^{tA} \left( \mathbf{x}\_0 + \mathbf{A}^{-1} (\mathbf{I} - e^{-tA}) \mathbf{q} + \frac{\overline{\mathbf{U}}}{\rho} (e^{\rho t} - \mathbf{1}) \mathbf{p}\_0 \right). \tag{13}$$

This yields an explicit form for the initial condition

$$p\_0 = \frac{\rho}{\overline{\mathcal{U}}(e^{\rho t} - 1)} \left( e^{-tA} \mathbf{x}(t) - \mathbf{x}\_0 - A^{-1} (I - e^{-tA}) \mathbf{q} \right) \tag{14}$$

which once again has to be normalized and together with *x*(*τ*) = *x*des yield

$$||e^{-tA}\mathbf{x}\_{\rm des} - \mathbf{x}\_0 - A^{-1}(I - e^{-At})q||^2 = \left(\frac{\overline{\mathbf{U}}}{\rho}(e^{\rho t} - 1)\right)^2. \tag{15}$$

When multiplying this formula by *e* −*ρτ*, the vector on the left-hand side can be simplified and the whole formula results in

$$0 = \left(\frac{\overline{\mathbf{U}}}{\rho} (\mathbf{1} - e^{-\rho \tau})\right)^2 - \left(-a\_1 \cos \omega \tau + a\_2 \sin \omega \tau + \beta\_1 e^{-\rho \tau}\right)^2$$

$$- \left(-a\_2 \cos \omega \tau - a\_1 \sin \omega \tau + \beta\_2 e^{-\rho \tau}\right)^2. \tag{16}$$

where

$$\begin{aligned} a\_1 &= \psi\_{\rm des,d} - \frac{\rho \psi\_{\rm pm,d} + \omega \psi\_{\rm pm,q}}{\rho^2 + \omega^2} \rho, & \quad \beta\_1 &= \psi\_{0,d} - \frac{\omega \psi\_{\rm pm,q} + \rho \psi\_{\rm pm,d}}{\rho^2 + \omega^2} \rho, \\\ a\_2 &= \psi\_{\rm des,q} + \frac{\omega \psi\_{\rm pm,d} - \rho \psi\_{\rm pm,q}}{\rho^2 + \omega^2} \rho, & \quad \beta\_2 &= \psi\_{0,q} + \frac{\omega \psi\_{\rm pm,d} - \rho \psi\_{\rm pm,q}}{\rho^2 + \omega^2} \rho. \end{aligned} \tag{17}$$

This is one equation in the time variable and thus of complexity comparable to that of (10). The only increase is due to the need to evaluate the exponential function. However, accurate approximations for its evaluation are available [27].

#### *3.2. Different d, q Inductances*

The assumption of the previous section is no longer valid in drives that have significantly different inductances. Since we use flux as the main variable, different inductances influence only the stator resistance dependent terms of Equation (3) via *RsL* −1 *s*,*d* and *RsL* −1 *s*,*q* . For analytical convenience, we denote *ρ* = 1 2 *Rs*(*L* −1 *s*,*d* + *L* −1 *s*,*q* ) and *δ* = 1 2 *Rs*(*L* −1 *s*,*d* − *L* −1 *s*,*q* ). Then, (3) has the form

$$\mathbf{A} = -\begin{bmatrix} \rho + \delta & 0\\ 0 & \rho - \delta \end{bmatrix} - \omega \mathbf{J}. \tag{18}$$

The form of the system trajectory via the matrix exponential (4) then depends on the relation between *δ* and *ω*. Specifically, if *ω* > *δ*, the trajectory is periodic in the base of sin and cos functions just like (11), but with a different frequency √ *ω*2 − *δ* 2 . When *δ* > *ω*, the trajectory becomes aperiodic with parameter √ *δ* 2 − *ω*2 . When *ω* = *δ*, the trajectory becomes a straight line. To achieve compact notation for all cases, we introduce auxiliary variables:

$$\begin{aligned} c\_1 &= \sqrt{|\omega^2 - \delta^2|}, \\ c\_2 &= \delta^2 - \omega^2 - \rho^2, \end{aligned}$$

and auxiliary functions *σ*(*t*) and *µ*(*t*), defined as:

$$\begin{cases} \sigma(t) = \frac{\sin(c\_1 t)}{c\_1} & \mu(t) = \cos(c\_1 t) & \text{for } |\delta| < |\omega|, \\\sigma(t) = \frac{\sinh(c\_1 t)}{c\_1} & \mu(t) = \cosh(c\_1 t) & \text{for } |\delta| > |\omega|, \\\sigma(t) = t & \mu(t) = 1 & \text{for } |\delta| = |\omega|. \end{cases} \tag{19}$$

Note that the last equation is a limit of both former cases for *c*1 → 0. Using auxiliary functions (19), the system dynamics can be written as

$$e^{tA} = e^{-\rho t} \begin{bmatrix} \mu(t) - \delta \sigma(t) & \omega \sigma(t) \\ -\omega \sigma(t) & \mu(t) + \delta \sigma(t) \end{bmatrix}. \tag{20}$$

In this case, the norm of the costate, ||*p*0||, required in (7) is rather complex and analytical expression of (5) becomes intractable. Therefore, we will proceed with the same normalization as in the case of PMSM, which is now an approximation:

$$e^{-sA}e^{-sA^T} \approx e^{2\rho t},\tag{21}$$

validity of this approximation will be tested in simulation.

Under this approximation, the solution of the costate equation becomes formally equal to (14), giving an implicit equation very similar to (16):

$$\begin{split} 0 = \left(\frac{\overline{\mathcal{U}}}{\rho}(1 - e^{-\rho\tau})c\_2\right)^2 - \left(a\_1\mu(\tau) + a\_2\sigma(\tau) + \beta\_1 e^{-\rho\tau}\right)^2 \\ &- \left(a\_3\mu(\tau) + a\_4\sigma(\tau) + \beta\_2 e^{-\rho\tau}\right)^2. \end{split} \tag{22}$$

where the only difference is in different constants

$$\begin{aligned} \alpha\_1 &= c\_2 \psi\_{\rm des,d} - \delta q\_1 + \omega \eta\_2 + q\_1 \rho\_{\prime} \\ \alpha\_2 &= q\_1 \omega^2 - q\_1 \delta^2 + \delta c\_2 \psi\_{\rm des,d} + \delta q\_1 \rho - \omega c\_2 \psi\_{\rm des,q} - \omega q\_2 \rho\_{\prime} \\ \alpha\_3 &= c\_2 \psi\_{\rm des,q} + \delta q\_2 - \omega q\_1 + q\_2 \rho\_{\prime} \\ \alpha\_4 &= q\_2 \omega^2 - q\_2 \delta^2 - \delta c\_2 \psi\_{\rm des,q} - \delta q\_2 \rho + \omega c\_2 \psi\_{\rm des,d} + \omega q\_1 \rho\_{\prime} \\ \beta\_1 &= -\omega q\_2 - c\_2 \psi\_{0,d} + \delta q\_1 - q\_1 \rho\_{\prime} \\ \beta\_2 &= -\delta q\_2 - c\_2 \psi\_{0,q} + \omega q\_1 - q\_2 \rho. \end{aligned} \tag{23}$$

Note that since these constants are independent of the time, they can be computed only once before the bisection method.

When the time of the transient *τ* is known, it is substituted to (14) to obtain *p*0, which is then substituted to (7) to obtain *u* ∗ (0):

$$
\mu\_d^\* = \frac{\rho}{(1 - e^{-\rho \tau}) c\_2} (\alpha\_1 \mu(\tau) + \alpha\_2 \sigma(\tau) + \beta\_1 e^{-\rho \tau}),
\tag{24}
$$

$$
\mu\_q^\* = \frac{\rho}{(1 - e^{-\rho\tau})c\_2} \left( a\_3 \mu(\tau) + a\_4 \sigma(\tau) + \beta\_2 e^{-\rho\tau} \right). \tag{25}
$$

#### *3.3. Implementation Details*

#### 3.3.1. Time Discretization

The proposed controller is designed in continuous time; however, its implementation will be done in discrete time. Let us denote the time instant at which we will apply the controller by time index *ut* . The time required to compute the solution prevents the use of the measurement at the same time. Hence, we assume that we have only one-step delayed measurement, *it*−1. Therefore, we have to perform delay compensation, i.e., the initial state *x*0 is calculated using one-step ahead prediction (2). This value and the requested state *x*req are substituted into the implicit Equation (16). The solution of this equation using the bisection method is a scalar value *τ* ∗ , which is substituted into (24) and (25), which is the time-optimal control law.

#### 3.3.2. Multiple Extremes in Bisection

To use the bisection method, we need to define the minimum and maximum of the search interval such that the function in these points have opposite signs. A suitable choice of the lower bound is zero, where the implicit function (22) is always positive. The upper bound is more problematic. One possible choice is setting the upper limit to the maximum expected time of a transient (e.g., 512 sampling periods, ∆*t*). This setting becomes problematic at a very high speed regime in the field weakening operation. At this region, the current trajectory starts circulating and may cross the requested value multiple

times in a short window. The implicit function then has multiple roots, see Figure 1, and we need to select the one closest to zero. An optimal solution would be to derive an analytic upper bound on it. However, in the experimental evaluation, we use a simple heuristic solution. Specifically, we choose *tmax* = 256∆*t*, but the first evaluation is not made in the middle of the interval but at point 10∆*t*. All subsequent evaluations are done using the standard bisection algorithm. We found that 20 steps of the bisection algorithm for solving Equation (22) are sufficient.

**Figure 1.** An example of multiple roots in the implicit Equation (22).

#### 3.3.3. Small Step Solution

Another problem arises when the trajectory becomes too close to the requested value and *τ* < ∆*t*. In such a case, we compute the reference value for the PWM using the conventional deadbeat solution, which arises from the discrete-time version of model (3)

$$\frac{\mathfrak{x}\_{t+1} - \mathfrak{x}\_{t}}{\Delta t} = \mathbf{A}\mathfrak{x}\_{t} + \mathfrak{u}\_{t} + \mathfrak{q}\_{t}$$

where it is assumed that the requested value *x* ∗ can be reached in one step, yielding

$$
\mu\_t^{\rm DB} = \frac{\mathbf{x}^\* - \mathbf{x}\_l}{\Delta t} - A\mathbf{x}\_l - q. \tag{26}
$$

When the state *xt* is close to the requested value, the deadbeat solution is within the voltage limit. However, it becomes suboptimal when the amplitude of the solution |*u* DB *t* | is higher than *U*. This motivates the proposed algorithm for evaluation of the current control *u* opt *t* , Algorithm 1.

The deadbeat solution alone is often used for current control [2]. However, in such a case, it has to be truncated to the voltage limit. Typical solution is

$$\mu\_t^{\rm DBt} = \begin{cases} \mu\_t^{\rm DB} & \text{if } |\mu\_t^{\rm DB}| < \overline{\mathcal{U}},\\ \frac{\overline{\mathcal{U}}}{|u\_t^{\rm DB}|} \mu\_t^{\rm DB} & \text{otherwise.} \end{cases} \tag{27}$$

which we will use for comparison with the proposed TOC. The key difference between DB and TOC is the situation when the available voltage cannot satisfy the requirement in one step (i.e., the "else" part of step 4 in Algorithm 1).

Note that both algorithms are feed-forward controllers, relying on model correctness. Potentially this may lead to a small steady-state error, which can be compensated by an additional feedback controller. However, we have not experienced any steady state error in our experiments.

#### **4. Simulations**

In this section, we study the difference between the proposed solution and previous solutions that neglect stator resistances and different inductances on the time-optimal trajectory. First, we use the parameters of our experimental rig and then we study the sensitivity of the solution to parameter variations. All simulations were done in Matlab, without the use of any toolbox, since the simulation is fully determined by elementary arithmetic operations.

#### *4.1. Nominal Parameters*

Our experimental rig has the parameters

$$\begin{aligned} R\_s &= 1.8 \,\Omega, & \quad \psi\_{pm,d} = 0.438 \,\text{Wb}, \\ L\_{s,d} &= 14.0 \,\text{mH}, & \Delta t = 100 \,\text{\mu s}, \\ L\_{s,q} &= 19.3 \,\text{mH}, & \quad \mathcal{U}\_{dc} = 450 \,\text{V}. \end{aligned} \tag{28}$$

We compare three different methods for time-optimal flux control: (i) control with neglected resistance [20,23], (ii) solution with equal inductances, where the inductance is set to mean value of the true machine parameters, *Ls* = 0.5(*Ls*,*d* + *Ls*,*q*), and (iii) the proposed method with full model parameters. The methods are compared on a step change of the requested stator current from zero to *i* ∗ *d* = −3A, *i* ∗ *q* = 14A. The current trajectories of different control strategies are displayed in Figure 2.

In two modes, open-loop (left) and closed-loop (right). The open loop trajectory applies the control strategy designed at the origin. In the closed loop, the strategy is recomputed at every sampling period. For comparison, the deadbeat control (27) is also computed and displayed in the right column of Figure 2. Note that all time-optimal control strategies are relatively close to each other. This indicates that the effect of resistance on the trajectory is low. In contrast, the deadbeat controller (27) takes a different path. The difference is in the time to reach the setpoints. While for *ω* = 10 rad/s all controllers reach the setpoint in 16 sampling periods, for *ω* = 400 rad/s, all time-optimal controllers reach the setpoint in 46 sampling periods (46∆*t*), but the deadbeat reaches the setpoint in 131∆*t*.

#### *4.2. Sensitivity Study*

The difference between all versions of the time-optimal control becomes more obvious when the ratio *ρ* between the resistance and inductance is higher. This is simulated by setting smaller values of the inductances, with *Ls*,*d* > *Ls*,*q*, namely *Ls*,*d* = 5 mH, *Ls*,*q* = 3 mH. The resulting current trajectories for all tested algorithms and step change of the requested stator current from zero to [5, 30]A are displayed in Figure 3.

**Figure 2.** Comparison of studied versions of time-optimal control and deadbeat control on stator current reference step change from zero to [−3, 14]A for *ω* = 10 rad/s (**top** row), *ω* = 120 rad/s (**middle** row), and *ω* = 400 rad/s (**bottom** row). The difference of the solutions in the origin is visible in the open-loop strategy (**left**), its impact on the receding horizon reevaluation in each sampling period (**right**), sampling times are denoted by dots on the full lines.

**Figure 3.** Comparison of studied versions of time-optimal control and deadbeat control on current reference step change from zero to [5, 30]A for modified parameters of the system with lower inductance at *ω* = 10 rad/s (**top** row), and *ω* = 250 rad/s (**bottom** row). The difference of the solutions in the origin is visible in open-loop strategy (**left**), its impact on the receding horizon reevaluation in each sampling period (**right**), sampling times are denoted by dots on the full lines.

In this case, simplifications of the time-optimal trajectory differ from the proposed solution. This is most obvious on the open-loop results, where only the proposed solution reaches the requested value, but the strategy with neglected resistance goes far to the left (negative *id* ) and the one with averaged inductances far to the right (positive *id* ). This tendency is corrected in the closed-loop due to recalculations of the trajectory in each time step. However, even with a correctly reached target, the strategy with neglected resistance generates a trajectory with very low *id* currents, and consequently reaches the setpoint in the longest time (15∆*t*). The strategy with averaged inductances takes a path with positive *id* currents and reaches the setpoint in 15∆*t*. The proposed strategy is at low speed almost equivalent to the deadbeat controller, and reaches the setpoint in 14∆*t* equally with the deadbeat.

An important conclusion is that even for these parameters, the open-loop trajectory of the proposed control reaches the target setpoint. This implies that the approximation proposed in (21) is sufficiently accurate, in contrast to the simplified solutions.

#### *4.3. Discussion of Results*

The results of the controllers significantly differ for different parameters. When the stator resistance is low, it can be neglected in the time-optimal controllers, and the resulting controller yields in high speeds significantly faster transients than the deadbeat controller. However, when the effect of the resistance is higher, the TOC with neglected resistance provides a poor solution that is actually slower than the deadbeat solution and yields higher Joule losses. In such a case, the proposed TOC with resistance yields a better solution. The

only computationally slow operation is the sinh function, which needs to be computed for *ω* < *δ*. If computational cost is a concern, it is safe to use the deadbeat controller in this operating regime and switch to TOC only for higher speeds.

### **5. Experimental Results for IPMSM**

A laboratory prototype of the IPMSM drive with the same parameters as in the simulation (28) was used to verify the approach experimentally. The test rig is displayed in Figure 4.

**Figure 4.** Photo of the test rig with the controlled IPMSM and loading induction machine (**left**), and the controlled converter, current sensors and the control board (**right**).

The rated power of the IPMSM drive is 4.5 kW, rated voltage 400 V, rated current 12.47 A rms, and rated speed 1500 rpm. The IPMSM drive is equipped with 12 bit absolute angular position encoder LARM ARC 405, torque sensor Burster 8661, voltage transducer LEM LV 25-P for converter dc-link voltage measurement, and current transducers LA 55-P for measurement of the stator phase currents. The switching frequency of the voltage-source converter supplying IPMSM is 10 kHz.

The current control strategies will be used to follow setpoint designed by optimal steady-state solution that respects stator resistance and different inductances as described in [4]. The full block diagram is displayed in Figure 5.

**Figure 5.** Block diagram of the tested closed loop controller.

This follows the conventional cascade structure, where the speed of the drive is controlled by a PI controller (or any other controller) yielding the requested torque, *T* ∗ *e* . The requested stator current vector in the *dq* reference frame with elements *i* ∗ *d* and *i* ∗ *q* is calculated from the requested torque by an optimization scheme that minimizes Joule losses in the steady-state [4]. Specifically, the optimal setpoint is calculated as the intersection of the torque curve with the MTPA curve or with the field weakening (FW) curve. The method presented in [4] extends previous approaches (e.g., [2]) by explicit consideration of stator resistance. The effect of the stator resistance on the FW curve is visualized in Figure 6.

**Figure 6.** Visualization of optimal setpoints for rotor speed *ω* = 400 rad/s. The optimal setpoint lies on the intersection of the FW curve with the iso-torque curve (T = 38.6 Nm). FW curve with correct resistance (*Rs* = 1.8 Ω) is compared to curve with neglected resistance (*Rs* = 0 Ω).

Note that the FW curve that does not consider that the resistance is symmetric around the *x*-axis, while the curve that considers that it is not. This means that the *id* current in the real drive has to be lower than predicted using the simplified FW curve. This mismatch has been achieved by using "safety" coefficient *ζ*, [2], *U* = *ζUdc*/2 where *Udc* is the dclink voltage. In real experiments, the safety coefficient for the simplified approach has to be set to *ζ* = 0.7 to achieve good performance. The improved solution of [4] allows us to use higher safety coefficient, *ζ* = 0.9 was found to be sufficient in our case. This improvement is achieved at a higher computational cost due to the need to solve roots of fourth-order polynomials.

In the experiments, the current control loop is either the proposed time-optimal controller (TOC), the standard truncated deadbeat controller (DB) from (27), or the FOC using the conventional PI controllers

$$u\_{d,t}^{\*} = k\_p e\_{d,t} + k\_i \sum\_{i=1} e\_{t-i\nu} \quad e\_{d,t} = (i\_{d,t} - i\_{d,t}^{\*})\_{\prime\prime}$$

and analogously for the *q* axis. Both the TOC and the deadbeat controllers are tuning free, but the performance of the PI controller depends on the choice of coefficients *kp* and *ki* (tuning). We have tuned the PI controller to yield the best overall performance. It is possible to re-tune the controller to obtain better results at one operating point, however, at the cost of deteriorating performance at another. However, the behavior of the control remains very similar for different tuning due to the fact that it is a pure feedback controller. The computational times of the control loop building-blocks are given in Table 1.

**Table 1.** Computational times of the control algorithm on DSP TMS320F28377S.

The execution time of the TOC depends on the rotor speed. For lower speeds, the base functions are the hyperbolic functions (sinh,cosh), which are more expensive, while for higher speeds, the base functions are the sin and cos functions, which are cheaper.

Since the TOC controller is based on the assumption of the perfect state knowledge, we use the Kalman filter (block KF in Figure 5) to reconstruct the speed, position, and stator currents from the position and current measurements.

The performance of the controller is tested on a testing profile displayed in Figure 7.

**Figure 7.** Speed profile of the experimental evaluation. Speed request is displayed in a dashed blue line, the actual speed in a red full line. Letters a, b, c and d indicate moments that are analyzed in detail in Figure 8.

Which is composed of multiple step changes of the requested rotor speed. Since the step change is high, the PI controller of the speed generates a request for the maximum possible torque. The optimization routine provides current setpoints corresponding to the best option at the given operating point. For low-speed operations, it selects points on the maximum torque per ampere curve. In the field weakening regime, it computes optimal currents for field weakening. References on the torque and currents remain almost constant for many sampling periods since the time-constant of the speed is much longer than that of the currents. To visualize the differences between the proposed controllers, we provide details of the current transient at the time of the step change at multiple operating points. These points are at different speeds *ω* = 0, 280, 400, 480 for acceleration request and *ω* = 600 rad/s for speed reversal request, see Figure 7 for illustration.

The resulting torque and current trajectories at selected operating points are displayed in Figure 8 in two modes: (i) time trajectories visualize the temporal convergence of the trajectories to setpoints displayed by dashed lines in the right subplots, and (ii) current *d*-*q* plane trajectories to the setpoints displayed by crosses in the left subplots. Note that at zero speed, the current trajectories of all controllers are very similar and equally fast. This is understandable since the TOC trajectory is almost a straight line at very low speeds. It starts deviating from a straight line with a higher speed, as can be noted for case (b) in Figure 8.

**Figure 8.** Details of torque and current trajectories at selected moments—(**a**–**d**)—of the experimental testing profile. For each transient, we plot the resulting trajectories for TOC controller (**top** row), PI controller (**middle** row), and deadbeat controller (**bottom** row). Trajectories of the currents in the *d*-*q* coordinates are displayed in the left subfigures with the requested current values denoted by crosses. Trajectories in the time domain are displayed in the right subfigures. The requested values of the torque and current are marked by dashed lines.

At *ω* = 280 rad/s, the field weakening limit crosses the MTPA line, and the reference of the *d* current moves to negative values. At this moment, the TOC controller is able to reach the requested torque in 2.5 ms and remain stable at this value. The PI controller reaches the requested torque around the same time, but only due to the fact that it reaches the torque curve at different current values than requested. In an attempt to reach the requested currents, the torque request overshot, and it takes a while to settle it. The DB controller reaches the requested torque later, around 3 ms, but it is also able to remain stable

at this value. Even more obvious is the difference at *ω* = 400 rad/s and *ω* = 480 rad/s, where the TOC controller is able to reach the requested torque in 4 ms, contrary to the other controllers, that are not able to reach it in twice that time. Note that the TOC controller starts the transient with decreasing *id* current at the cost of decreasing the *iq* current. This strategy yields a decrease of the torque for a very short time. However, this short time decrease allows a much faster increase of the torque in the second part of the transient. This demonstrates the ability of the TOC controller to optimize on a very long horizon. Performance of all controllers at the speed reversal moment at *ω* = 600 rad/s is again comparable.

A summary of the settling times for all transients is provided in Table 2. The settling time was measured as the time after which the measured current is within 5% of the requested value.

**Table 2.** Comparison of settling time of different controllers for the transients in Figure 8. All times are given in milliseconds.

As expected, the TOC has the most visible benefits at high speeds. While minimization of the speed of the transient is its primary objective, we would like to point out that the current trajectories of the TOC controller also exhibit lower oscillations. This is propagated into the torque trajectory, which is followed by TOC with lower oscillations than those provided by the PI and deadbeat controllers.

#### **6. Conclusions**

We have proposed a time-optimal control strategy for current control of the general synchronous motor drive, that considers also stator resistance and differences in stator inductances in the direct and quadrature axis. We have shown that the base functions of the optimal trajectory differ from those of the previous simplified solutions. We have derived a bisection-based optimization algorithm for the evaluation of the more accurate control strategy, which has only a minor increase in computational demands compared to previous simplified solutions.

The proposed strategy was able to reach the requested current in the minimum possible time. The difference to conventional PI or DB controllers is negligible at lower speeds; however, it provides significantly faster transients at higher speeds. On the testing prototype, TOC achieved a four times faster settling time than the PI and DB controllers at rotor speed of 480 rad/s. This is caused by the fact that TOC is a feed-forward controller increasing the tracking error at the beginning of the transient, which is not natural for a feedback controller.

While the controller is designed to minimize the time of the transient, we have found that it also provides more stable tracking of the required torque. The use of this control may thus benefit applications that require fast and accurate torque tracking such as high dynamic servo drives for robotics or manipulators.

Note that the proposed feed-forward strategy can be combined with any feed-back strategy, even with the conventional PI cascade. We believe that the proposed controller can be used as a universal building block complementing existing solutions. Exploring all potential benefits of the approach in combinations with different feedback solutions and different observers is left for future study.

**Author Contributions:** Methodology, V.Š. and L.A.; software, L.A. and A.G.; validation, A.G.; data curation, A.G.; writing—original draft preparation, V.Š., L.A. and A.G.; writing—review and editing, Z.P. and L.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Education, Youth and Sports of the Czech Republic under the project OP VVV Electrical Engineering Technologies with High-Level of Embedded Intelligence, CZ.02.1.01/0.0/0.0/18 069/0009855, project OP VVV Research Center for Informatics, CZ.02.1.01/0.0/0.0/16 019/0000765, and by UWB Student Grant Project no. SGS-2021-021.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Implementation and Control of a Wheeled Bipedal Robot Using a Fuzzy Logic Approach**

**Chun-Fei Hsu 1,\* , Bo-Rui Chen 2 and Zi-Ling Lin 1**

**Abstract:** This study designs and implements a wheeled bipedal robot (WBR) that combines the mobility of wheeled robots and the dexterity of legged robots. The designed WBR has extra knee joints to maintain body balance when encountering uneven terrain. Because of the robot's highly nonlinear, dynamic, unstable, and under-actuated nature, an intelligent motion and balance controller (IMBC) based on a fuzzy logic approach is proposed to maintain the balance of the WBR while it is standing and moving on the ground. It should be emphasized that the proposed IMBC system does not require prior knowledge of system dynamics and the controller parameters are tuned using the qualitative aspects of human knowledge. Furthermore, a 32-bit microcontroller that has memory, programmable I/O peripherals, and a processor core is used to implement the IMBC method. Finally, moving and rotating, height-changing, posture-keeping, and "one leg on slope" movement scenarios are tested to demonstrate the feasibility of the proposed IMBC system. The experimental results show that, by using the proposed IMBC system, the WBR can not only balance and move well both on flat ground and in complex terrain but also extend each leg independently to maintain body balance.

**Keywords:** fuzzy control; balance control; movement control; wheeled bipedal robot

### **1. Introduction**

Automation has become an integral part of all kinds of work, touting the benefits of not only reducing personnel costs but also maintaining the quality of the work done. In the past, industrial robots and human workers worked independently of each other in the workspace; however, robots and human workers will need to interact more directly in the future. Mechanisms for robot motion can be divided into two main domains: wheel-based approaches [1–3] and leg-based approaches [4–6]. Wheeled robots can change their position easily and quickly and their cost is lower; however, they cannot handle uneven terrain. Legged robots can travel almost anywhere but require heavy computational processing to actuate such complex movements. As they can achieve the advantages of both wheeland leg-based designs, wheel-legged robots have attracted much attention [7–12]. The wheel-legged robot can move fast on flat ground and pass over uneven terrain, thereby improving the robot's adaptability in application scenarios. It also allows the robot's height to be changed to avoid collisions when it needs to pass under low obstacles.

Both academia and industry have conducted extensive research on self-balancing two-wheeled robots [13–15]. These have become an important member of the family of autonomous service robots, thanks to their good movement speed, load capacity, and terrain adaptability, and the concept has been successfully applied to transport robots [16] and hospitality robots [17,18]. Considering the diversity of locations, the wheeled bipedal robot (WBR) [19,20], which has two legs and two wheels (with the wheels installed on each foot), has more degrees of freedom than self-balancing two-wheeled robots [13–18] in the vertical direction. By achieving a higher movement speed and a larger movement

**Citation:** Hsu, C.-F.; Chen, B.-R.; Lin, Z.-L. Implementation and Control of a Wheeled Bipedal Robot Using a Fuzzy Logic Approach. *Actuators* **2022**, *11*, 357. https://doi.org/ 10.3390/act11120357

Academic Editor: Matteo Cianchetti

Received: 29 October 2022 Accepted: 30 November 2022 Published: 2 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 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 (https:// creativecommons.org/licenses/by/ 4.0/).

range, it is designed to be flexible even in the variable terrain of indoor environments. The most famous WBRs are "Handle" from Boston Dynamics and "Ascento" from ETH Zurich, which can move quickly on flat terrain, go down stairs, and jump through obstacles. However, because there are only two contact points between the robot and the ground, the ability to balance plays a key role in the practical application of WBRs.

Recently, there is increasing interest in studying the design, modeling, and control of WBRs. Traditional control schemes have been employed to solve the balance and height control issues [21–32]. A linear quadratic regulator (LQR) controller was designed to implement stabilization and driving control, and a proportional-integral-derivative (PID) controller was designed to synchronize the two leg motors in order to change the height of the robot [21]. However, these control gains needed to be pre-constructed using a tuning procedure. A cascade PID controller was proposed as the movement control system of a WBR [22]. Although it could achieve satisfactory balancing motion performance, it did not take into account how the WBR could pass over rough terrain. An easy-to-implement LQR-based control system was proposed [23–26]. However, the control performance degrades in the case of WBR operating point changes. An inverse dynamics controller was designed to consider the full dynamics of the WBR without neglecting any nonlinear dynamics terms [27]. However, these system uncertainties are unavoidable in a real WBR. A whole-body control frame was proposed for a WBR so that it could not only accelerate and decelerate along the forward direction and balance dynamically but also withstand disturbances more efficiently [28–32]. However, its control performance was highly dependent on the accuracy of the WBR's dynamic model.

The parameters of the model-based controller in recent studies [21–32] have had to be manually adjusted to ensure effective control when the physical parameters of the WBR change, such as robot height and payload weight. To overcome this shortcoming as well as solve the WBR balance problem, intelligent control schemes have been proposed [33–35]. A cascade PID controller was proposed for the path tracking of a WBR in which a reinforcement learning method automatically obtained optimal gains for the control law in a simulation [33]. Although the obtained controller gains could be directly applied in real-world scenarios without any re-training, the effect of robot height variation on the control response was not considered. A learning-based solution was proposed based on adaptive optimal control, where the controller parameters could be learned directly from the input-state data in the Gazebo robotic simulator [34]. However, it required heavy computational processing. An intelligent adaptive sliding-mode controller was proposed for the trajectory tracking and stability of a WBR, where a fuzzy-basis-function-network was used to approximate the system dynamics online [35]. However, the controller did not consider the effect of WBR height variation on the system response.

Fuzzy control (FC) using linguistic information can model the qualitative aspects of human knowledge and possesses several advantages, such as robustness, being model-free, and being a rule-based algorithm [36–38]. Most of the operations in an FC use error and change-of-error as the fuzzy input variables, where the rules table is constructed in a two-dimensional input space, which cause difficulties in implementation. It can be found that the rule tables have the skew–symmetry property, where the absolute magnitude of the control input is proportional to the distance from its main diagonal line in the normalized input space. To reduce the complexity of the implementation, a single-input fuzzy controller (SFC) was proposed based on the sliding-mode control scheme, where a new fuzzy input variable (called the signed distance) was derived [39–42]. Compared to FC system, the total number of fuzzy rules used in SFC system is greatly reduced, so the control rules can be easy to generate and adjust. However, due to the linguistic expression of the FC and SFC systems, it has not been easy to guarantee the stability and robustness of the closed-loop control systems.

This study proposes a WBR that consists of a body and two legs with wheels where each leg can be independently extended and retracted by driving a motor mounted on the hip of the body. It has the advantages of both leg and wheel movement, with good maneuverability and a high load capacity. The height of the proposed WBR can be adjusted through the realization of a four-bar linkage. Since the four-bar linkage approximates the linear motion of the wheels perpendicular to the ground, the WBR can traverse simple terrain in the same manner as wheeled robots and adjust its posture to the shape of uneven ground.

To overcome the highly dynamic and under-actuated nature of the proposed WBR's behavior, this study proposes an intelligent motion and balance control (IMBC) system that comprises fuzzy movement balancing control (FMBC), fuzzy yaw steering control (FYSC), and fuzzy roll balancing control (FRBC). The FMBC, based on a dual-loop control structure, is designed to cope with the under-actuated nature problem of the WBR and allow the WBR to move at the desired velocity command while maintaining balance. The FYSC is a compensated speed control signal for the WBR that can rotate about its vertical axis, and the FRBC extends each leg independently to change the robot's height and altitude when driving over challenging terrain. Since the proposed IMBC method was designed based on the SFC method, it does not require complex analysis or heavy computational processing. A low-cost 32-bit STM32F446RE microcontroller is used to implement the control method for the WBR. Finally, four experimental scenarios are tested to explore the stability performance of the proposed IMBC system. The experimental results show that the proposed IMBC system is feasible and effective for robot motion on both flat and complex terrain.

The contribution of this study can be summarized as follows: (1) The mechanical structure of the WBR, which has two legs and two wheels (with the wheels installed on each foot), is designed and implemented. (2) An experimental setup is presented for the WBR using a low-cost 32-bit STM32F446RE microcontroller. (3) The IMBC system is shown to be able to maintain the balance of the WBR while standing and moving on the ground. (4) The IMBC system is shown to be able to independently control two leg motors to change the robot's height and attitude and allow the WBR to move over challenging terrain.

The remainder of this study is outlined as follows. In Section 2, the design and dynamic model of the WBR is presented. In Section 3, the structure of the IMBC system is presented. Experimental tests are presented in Section 4, and conclusions are given in Section 5.

### **2. System Structure and Modeling of the WBR**

#### *2.1. Mechanical Design*

Figure 1 shows the mechanical structure of the WBR: Figure 1a shows the leg components and Figure 1b shows the kinematic models. The developed WBR, which consists of a body and two legs with wheels, is a promising way to improve the performance of bipedal robots in terms of flexibility, movement speed, and energy efficiency. By adding wheels to the ends of the legs, each leg can be independently extended and retracted by driving a leg motor mounted on the hip of the body. The torsion springs installed in the inner joints (see Figure 1a) are used to counteract the weight of the WBR itself, which can reduce the control effort of the leg motors when driving or standing to improve the overall efficiency of the robot. The total height of the WBR can be adjusted between 26 cm and 33 cm through the realization of a four-bar linkage.

In order to illustrate the relationship between the leg length of the WBR and the angle of the leg motor, part of the right leg is considered as shown in Figure 1b. The length between the pin joint and the inner joint, *LkR*, can be found according to the Cosine law as follows:

$$L\_{\rm kR} = \sqrt{L\_1^2 + L\_2^2 - 2L\_1L\_2 \cos \theta\_{\rm kR}} \tag{1}$$

where *L*1 is the length between the pin joint and the knee joint, *L*2 is the length between the pin joint and the leg joint, and *θkR* is the leg motor angle of the right leg. The torsion spring angle *θtR* and the length *LR* of the right leg of the WBR can be found as follows:

$$\theta\_{lR} = \pi - \cos^{-1}\left(\frac{L\_1^2 + L\_{kR}^2 - L\_2^2}{2L\_1L\_{kR}}\right) - \cos^{-1}\left(\frac{L\_3^2 + L\_{kR}^2 - L\_1^2}{2L\_3L\_{kR}}\right) \tag{2}$$

$$L\_R = \sqrt{L\_1^2 + L\_4^2 - 2L\_1L\_4\cos(\theta\_{lR})}\tag{3}$$

where *L*3 is the length between the knee joint and the inner joint and *L*4 is the length between the inner joint and the wheel motor. The height of the right side of the WBR is given as:

$$h\_R = L\_R + R \tag{4}$$

where *R* is the wheel radius. Similarly, the height of the left side of the WBR, *hL*, can be obtained using the leg motor angle of the left leg *θkL* as follows:

$$h\_L = \sqrt{L\_1^2 + L\_4^2 - 2L\_1L\_4\cos(\theta\_{lL})} + R \tag{5}$$

in which

$$\theta\_{tL} = \pi - \cos^{-1}\left(\frac{L\_1^2 + L\_{kL}^2 - L\_2^2}{2L\_1L\_{kL}}\right) - \cos^{-1}\left(\frac{L\_3^2 + L\_{kL}^2 - L\_1^2}{2L\_3L\_{kL}}\right) \tag{6}$$

$$L\_{kL} = \sqrt{L\_1^2 + L\_2^2 - 2L\_1L\_2\cos\theta\_{kL}}\tag{7}$$

where *θtL* is the torsion spring angle and *LkL* is the length between the pin joint and the inner joint of the left leg of the WBR, respectively. Thus, the joint angles are adjusted to produce a change in the lengths of both sides of the robot (*hR* and *hL*) when the two leg motors are driven (*θkR* and *θkL*). *Actuators* **2022**, *11*, x FOR PEER REVIEW 4 of 18

**Figure 1.** Mechanical structure. (**a**) Leg components. (**b**) Kinematic model of the right leg. **Figure 1.** Mechanical structure. (**a**) Leg components. (**b**) Kinematic model of the right leg.

#### In order to illustrate the relationship between the leg length of the WBR and the angle *2.2. Hardware and Software*

given as:

in which

where

θ

motors are driven (

*2.2. Hardware and Software*

of the leg motor, part of the right leg is considered as shown in Figure 1b. The length between the pin joint and the inner joint, *LkR* , can be found according to the Cosine law as follows: *LkR L L L L* θ *kR* 2 cos 1 2 2 2 2 = 1 + − (1) where *L*1 is the length between the pin joint and the knee joint, *L*2 is the length between the pin joint and the leg joint, and θ*kR* is the leg motor angle of the right leg. The torsion spring angle θ*tR* and the length *LR* of the right leg of the WBR can be found as follows: + − − + − = − − − *kR kR kR kR tR L L L L L L L L L L* 3 2 1 2 2 1 3 1 2 2 2 2 1 1 2 cos 2 θ π cos (2) Figure 2 shows the photographs of the WBR: Figure 2a shows the WBR platform and Figure 2b shows the hardware connection diagram. The hardware connection diagram is composed of a STM32F446RE microcontroller, a MPU 6050, two wheel motors with optical encoders, and two leg motors. The MPU 6050, consisting of a 3-axis accelerometer and a 3-axis gyroscope, is located at the center of the robot body, the STM32F446RE microcontroller is placed on the upper body, and the motor drives are located on the lower body. The leg motor is a Dynamixel motor (MX-64) made by ROBOTIS, which can provide enough power for a leg to adjust the height of the robot. The wheel motor is a geared DC motor (GM25-370) made by ChiHai MOTOR, which can provide smooth stabilization and counteract disturbances to the tilt angle of the system. An encoder is installed behind the wheel motor to provide information on the robot's translational velocity and position. The MPU 6050 is used to measure the roll and pitch angles of the robot tilt. The STM32F446RE microcontroller, which is coded by Keil C, is used as a platform

> 2 4 2 *LR* = *L*1 + *L* − *L L*

between the inner joint and the wheel motor. The height of the right side of the WBR is

where *R* is the wheel radius. Similarly, the height of the left side of the WBR, *hL* , can be

*hL* = *L* + *L* − 2*L*1*L*4 cos( *tL* ) + *R* 2 4 2 1

*kL kL tL L L*

> 2 2

inner joint of the left leg of the WBR, respectively. Thus, the joint angles are adjusted to produce a change in the lengths of both sides of the robot ( *hR* and *hL* ) when the two leg

*L L L L L*

1

*LkL L L L L*

*tL* is the torsion spring angle and *LkL* is the length between the pin joint and the

2 2 θ

2 cos( ) 1 4

θ

*kL* as follows:

θ

 + − −

θ*kL* 2 cos 1 2

cos

*tR* (3)

(5)

2 1

*kL kL*

2 = 1 + − (7)

*L L L*

2 2 1 3

2

cos (6)

3

*hR* = *LR* + *R* (4)

+ − = − − −

2

2 2 1 1

obtained using the leg motor angle of the left leg

θ

θ*kR* and π

θ*kL* ).

for the implementation of control algorithms; it has a Cortex-M4 core, a processor clock frequency of up to 180 MHz, and a rich set of IO ports for connecting to various sensors and motors. The microcontroller, leg motors, wheel motors, and remaining electronic devices are powered by a 3-cell LiPo Battery. control algorithms; it has a Cortex-M4 core, a processor clock frequency of up to 180 MHz, and a rich set of IO ports for connecting to various sensors and motors. The microcontroller, leg motors, wheel motors, and remaining electronic devices are powered by a 3-cell LiPo Battery.

Figure 2 shows the photographs of the WBR: Figure 2a shows the WBR platform and Figure 2b shows the hardware connection diagram. The hardware connection diagram is composed of a STM32F446RE microcontroller, a MPU 6050, two wheel motors with optical encoders, and two leg motors. The MPU 6050, consisting of a 3-axis accelerometer and a 3-axis gyroscope, is located at the center of the robot body, the STM32F446RE microcontroller is placed on the upper body, and the motor drives are located on the lower body. The leg motor is a Dynamixel motor (MX-64) made by ROBOTIS, which can provide enough power for a leg to adjust the height of the robot. The wheel motor is a geared DC motor (GM25-370) made by ChiHai MOTOR, which can provide smooth stabilization and counteract disturbances to the tilt angle of the system. An encoder is installed behind the wheel motor to provide information on the robot's translational velocity and position. The MPU 6050 is used to measure the roll and pitch angles of the robot tilt. The STM32F446RE microcontroller, which is coded by Keil C, is used as a platform for the implementation of

*Actuators* **2022**, *11*, x FOR PEER REVIEW 5 of 18

**Figure 2.** Photographs of the WBR. (**a**) WBR platform. (**b**) Hardware connection diagram. **Figure 2.** Photographs of the WBR. (**a**) WBR platform. (**b**) Hardware connection diagram.

#### *2.3. Modeling 2.3. Modeling*

In order to manipulate the WBR, it is necessary to build a mathematical model that allows the robot and the designed controller to successfully achieve the desired control goals. Figure 3 shows the coordinate systems of the WBR and Appendix A lists all the symbols used and their definitions. The robot's position can be measured as ( ) 2 1 θ *w* = θ *R* + θ *L* and the robot's body yaw angle can be given as: In order to manipulate the WBR, it is necessary to build a mathematical model that allows the robot and the designed controller to successfully achieve the desired control goals. Figure 3 shows the coordinate systems of the WBR and Appendix A lists all the symbols used and their definitions. The robot's position can be measured as *θw* = 1 2 (*θR* + *θL*) and the robot's body yaw angle *ψ* can be given as:

ψ

$$
\psi = \frac{R}{2W} (\theta\_R - \theta\_L) \tag{8}
$$

where *θR* is the rotary angle of the right wheel motor and *θL* is the rotary angle of the left wheel motor. The dynamic model of the WBR for the pitch axis can be described as [43]:

$$\mathbf{E}\begin{bmatrix}\ddot{\theta}\_w\\ \ddot{\theta}\end{bmatrix} + \mathbf{F}\begin{bmatrix}\dot{\theta}\_w\\ \dot{\theta}\end{bmatrix} + \mathbf{G}\begin{bmatrix}\theta\_w\\ \theta\end{bmatrix} = \mathbf{H}\begin{bmatrix}v\_L\\ v\_R\end{bmatrix} + \mathbf{D} \tag{9}$$

where **E** is the inertia matrix, **F** is the matrix of centripetal and Coriolis forces, **G** is the gravity matrix, **H** is the gain matrix, **D** is the lump of uncertainties matrices (including the external disturbances and parameter uncertainties), and *vL* and *vR* are the input voltages of the left and right wheel motors, respectively. These matrices can be found by considering the limits *θ* → 0, sin *θ* → *θ* and cos *θ* → 1 as:

$$\mathbf{E} = \begin{bmatrix} (2m+M)R^2 + 2J\_W + 2J\_m & MLR - 2J\_m \\ MLR - 2J\_m & ML^2 + J\_\Psi + 2J\_m \end{bmatrix} \tag{10}$$

$$\mathbf{F} = \mathbf{2} \begin{bmatrix} \frac{K\_l K\_b}{R\_{\text{ll}}} & -\frac{K\_l K\_b}{R\_{\text{ll}}} \\ -\frac{K\_l K\_b}{R\_{\text{m}}} & \frac{K\_l K\_b}{R\_{\text{m}}} \end{bmatrix} \tag{11}$$

$$\mathbf{G} = \begin{bmatrix} 0 & 0 \\ 0 & -MgL \end{bmatrix} \tag{12}$$

$$\mathbf{H} = \begin{bmatrix} \frac{K\_t}{R\_m} & \frac{K\_t}{R\_m} \\ -\frac{K\_t}{R\_m} & -\frac{K\_t}{R\_m} \end{bmatrix} \tag{13}$$

**Figure 3.** Generalized coordinates of the system. (**a**) Side view. (**b**) Top view. (**c**) Front view. **Figure 3.** Generalized coordinates of the system. (**a**) Side view. (**b**) Top view. (**c**) Front view.

**3. Controller Design for the WBR**  A dual-loop control structure was proposed to cope with the under-actuated nature problem of a self-balancing two-wheeled robot moving and turning while maintaining the robot's balance [43]. Extending the design ideas, this study proposes an IMBC system for the WBR, as shown in Figure 4. It comprises three controllers: FMBC, FYSC, and FRBC. Equation (9) can be divided into two subsystems: one is the velocity control subsystem and the other is the balancing control subsystem. The velocity control subsystem is used to allow the WBR to move at the desired velocity command, and the balancing control subsystem is used to maintain the WBR's balance. Meanwhile, the dynamic model of the WBR for the yaw axis can be expressed as [43]:

$$\bar{\psi} = -\frac{J}{I} + \left[ -\frac{K}{I} \frac{K}{I} \right] \left[ \begin{array}{c} v\_L \\ v\_R \end{array} \right] \tag{14}$$

design that independently controls two leg motors to change the robot's height and altitude and allows the WBR to travel over challenging terrain. where *I* = 2*mW*2 + *Jφ* + 2*W*2 *R*2 (*Jw* + *Jm*), *J* = 2*W*2 *R*2 *KtKb Rm* , and *K* = *W R Kt Rm* . The robot body roll angle can be given as:

$$\phi = \tan^{-1}(\frac{h\_L - h\_R}{2W}) \tag{15}$$

leg motor *b v* θ *kL* θ *kR* , FRBC φ *e* + θ intelligent motion and balance control *h* When the robot's two legs are on the ground. From (9), (14) and (15), it shows the nonlinear and under-actuated nature of the WBR system. Since these robot parameters and the uncertainties matrix are unknown or perturbed in real WBR applications, a model-free control scheme should be designed that not only accelerates and decelerates along the forward direction and balances dynamically but is also robust to uncertainties.

#### wheel motor transform *L R* FMBC *v* ,*v* **3. Controller Design for the WBR**

kinematic

*3.1. FMBC Design*

θ *wc*

> ψ*c*

φ*c*

−

*w* + *e* − θ *w*

−

ψ

is designed as [43]:

According to the velocity command

wheeled bipedal robot FYSC *y v* ψ ψ*e* + A dual-loop control structure was proposed to cope with the under-actuated nature problem of a self-balancing two-wheeled robot moving and turning while maintaining the robot's balance [43]. Extending the design ideas, this study proposes an IMBC system for the WBR, as shown in Figure 4. It comprises three controllers: FMBC, FYSC, and FRBC. It

θ

θ *w*

given by the users, an outer-loop sliding index

θ*wc*

The FMBC should control the WBR to move at the desired velocity command and remain balanced even in the height variation of the WBR. A dual-loop control scheme is utilized, where the outer-loop controller (for the velocity control subsystem) is used to calculate the desired robot pitch angle for tracking velocity commands and the inner-loop controller (for the balancing control subsystem) is applied to keep the body balanced with the desired robot pitch angle. Thus, the advantage of the FMBC scheme lies in its simplicity and robust performance, compared to some other existing schemes [23–32].

is similar to [43], where the FMBC designed using a dual-loop control structure allows the WBR to move at the desired velocity command while maintaining balance and the FYSC can control the WBR to rotate about its vertical axis. Unlike [43], the FRBC is a novel design that independently controls two leg motors to change the robot's height and altitude and allows the WBR to travel over challenging terrain. It is similar to [43], where the FMBC designed using a dual-loop control structure allows the WBR to move at the desired velocity command while maintaining balance and the FYSC can control the WBR to rotate about its vertical axis. Unlike [43], the FRBC is a novel design that independently controls two leg motors to change the robot's height and altitude and allows the WBR to travel over challenging terrain.

**Figure 3.** Generalized coordinates of the system. (**a**) Side view. (**b**) Top view. (**c**) Front view.

(b)

ψ

*D* θ*R* *hL*

*hL* − *hR*

(c)

φ

*W*

*hR*

*W*

φ

*W*

2*W*

*W*

A dual-loop control structure was proposed to cope with the under-actuated nature problem of a self-balancing two-wheeled robot moving and turning while maintaining the robot's balance [43]. Extending the design ideas, this study proposes an IMBC system for the WBR, as shown in Figure 4. It comprises three controllers: FMBC, FYSC, and FRBC.

*Actuators* **2022**, *11*, x FOR PEER REVIEW 7 of 18

θ*L*

**3. Controller Design for the WBR**

θ

*L*

(a)

θ *w* θ *L* θ*R* , ,

θ *kL* θ*kR* ,

**Figure 4.** The block diagram of the IMBC system for the WBR. **Figure 4.** The block diagram of the IMBC system for the WBR.

#### *3.1. FMBC Design 3.1. FMBC Design*

The FMBC should control the WBR to move at the desired velocity command and remain balanced even in the height variation of the WBR. A dual-loop control scheme is utilized, where the outer-loop controller (for the velocity control subsystem) is used to calculate the desired robot pitch angle for tracking velocity commands and the inner-loop controller (for the balancing control subsystem) is applied to keep the body balanced with the desired robot pitch angle. Thus, the advantage of the FMBC scheme lies in its simplicity and robust performance, compared to some other existing schemes [23–32]. According to the velocity command θ *wc* given by the users, an outer-loop sliding index is designed as [43]: The FMBC should control the WBR to move at the desired velocity command and remain balanced even in the height variation of the WBR. A dual-loop control scheme is utilized, where the outer-loop controller (for the velocity control subsystem) is used to calculate the desired robot pitch angle for tracking velocity commands and the inner-loop controller (for the balancing control subsystem) is applied to keep the body balanced with the desired robot pitch angle. Thus, the advantage of the FMBC scheme lies in its simplicity and robust performance, compared to some other existing schemes [23–32]. According to the velocity command . *θwc* given by the users, an outer-loop sliding index is designed as [43]:

$$s\_{\overline{w}} = e\_{\overline{w}} + \lambda\_1 \int\_0^t e\_{\overline{w}} d\tau \tag{16}$$

where *ew* = . *θwc* − . *θw* is the translational velocity error and *λ*1 is a positive constant. In order to converge the translational velocity error of the robot to zero quickly and have an excellent system response such as a faster rising time and smaller overshoot, each fuzzy rule of the outer-loop controller is given as follows:

$$\text{IF } s\_w \text{ is } F\_{sw}^i \text{ THEN } \theta\_c \text{ is } F\_{\theta c'}^i \text{ for } i = 1, 2, \dots, 7 \tag{17}$$

where *F i sw* are the labels of the fuzzy sets, *F i θc* are the singleton control actions, and *θc* is the desired robot pitch angle. The fuzzy rules base of the outer-loop controller is given in Table 1, where the fuzzy labels are negative big (NB), negative medium (NM), negative small (NS), zero (ZO), positive small (PS), positive medium (PM), and positive big (PB), as shown in Figure 5a, where the membership functions of the fuzzy sets are given in a triangular form. The fuzzy rules in Table 1 are constructed using the idea that the outerloop sliding index *sw* can quickly reach zero via a time-consuming trial-and-error tuning procedure. The adopted outer-loop sliding index reduces the number of fuzzy rules, thus reducing the complexity of adjustment. It shows that the translational velocity error *ew* will converge to zero when the outer-loop sliding index converges to zero. The outer-loop controller allows the WBR to move at the desired velocity command; thus, the velocity control subsystem can be asymptotically stable under the action of the outer-loop controller.

**Table 1.** The fuzzy rules for the FMBC. **Table 1.** The fuzzy rules for the FMBC. *ws* **NB NM NS ZO PS PM PB**

controller.

*i F*θ

where *w wc w e* θ

θ

*Actuators* **2022**, *11*, x FOR PEER REVIEW 8 of 18

rule of the outer-loop controller is given as follows:

where *i Fsw* are the labels of the fuzzy sets, *i F*

= − is the translational velocity error and

IF *ws* is *i Fsw* , THEN

**Figure 5.** Membership functions. (**a**) FMBC and FRBC; (**b**) FYSC. **Figure 5.** Membership functions. (**a**) FMBC and FRBC; (**b**) FYSC.

To cope with the system uncertainties caused by the changes in robot height, each fuzzy rule of the inner-loop controller is given as follows: To cope with the system uncertainties caused by the changes in robot height, each fuzzy rule of the inner-loop controller is given as follows:

*sw* **NB NM NS ZO PS PM PB**

*c* −15.0 −8.5 −5.4 0.0 5.4 8.5 15.0

 = + *t sw ew ewd* 0 1 λ

order to converge the translational velocity error of the robot to zero quickly and have an excellent system response such as a faster rising time and smaller overshoot, each fuzzy

> θ *c* is *i F*θ

θ

is the desired robot pitch angle. The fuzzy rules base of the outer-loop controller is given in Table 1, where the fuzzy labels are negative big (NB), negative medium (NM), negative small (NS), zero (ZO), positive small (PS), positive medium (PM), and positive big (PB), as shown in Figure 5a, where the membership functions of the fuzzy sets are given in a triangular form. The fuzzy rules in Table 1 are constructed using the idea that the outerloop sliding index *ws* can quickly reach zero via a time-consuming trial-and-error tuning procedure. The adopted outer-loop sliding index reduces the number of fuzzy rules, thus reducing the complexity of adjustment. It shows that the translational velocity error *we* will converge to zero when the outer-loop sliding index converges to zero. The outer-loop controller allows the WBR to move at the desired velocity command; thus, the velocity control subsystem can be asymptotically stable under the action of the outer-loop

τ

λ

(16)

*c* , for *i* = 1,2,..,7 (17)

*c* are the singleton control actions, and

1 is a positive constant. In

θ*c*

$$\text{IF } h \text{ is } F\_{h'}^{j} \text{ THEN } v\_b \text{ is } \mathbf{k}\_j^T \mathbf{e}\_{\Theta'} \text{ for } j = 1, 2, \dots, 3 \tag{18}$$

where the robot height can be measured as *h* = 1 2 (*hR* + *hL*), *F j h* are the labels of the fuzzy sets, *eθ* = *θ* − *θc* is the error of the body pitch angle, **e**θ = [*eθ* , . *eθ* ] *T* is the error vector, and **k***j* is the control gains vector of the inner-loop controller. The fuzzy labels used in this study are low (LO), medium (MD), and high (HI), as shown in Figure 5b, where the membership functions of the fuzzy sets are given in a triangular form. It is known that the outer-loop controller obtains the virtual angle command by defining the outer-loop sliding index and the inner-loop controller tracks the virtual angle command by defining the robot angle error. Thus, the balancing control subsystem can be asymptotically stable under the action of the inner-loop controller.

### *3.2. FYSC Design*

FYSC is designed to control the WBR so that it can rotate along its vertical axis. The speed control signals applied to the right or left driving wheels are obtained as follows:

$$
v\_L = v\_b + v\_y \tag{19}$$

$$
v\_R = v\_b - v\_y \tag{20}$$

where *vb* is the FMBC output and *vy* is the FYSC output. According to the velocity command . *ψc* given by the users, a yawing sliding index is designed as [43]:

$$s\_{\psi} = e\_{\psi} + \lambda\_2 \int\_0^t e\_{\psi} d\tau \tag{21}$$

where the yaw speed error *eψ* = . *ψc* − . *ψ* and *λ*2 are positive constants. Each fuzzy rule of the FYSC is given as follows:

$$\text{IF } s\_{\psi} \text{ is } F\_{s\psi}^{i} \text{ THEN } v\_{y} \text{ is } F\_{vy}^{i} \text{ for } i = 1, 2, \dots, 7 \tag{22}$$

where *F j sψ* are the labels of the fuzzy sets as shown in Figure 5a and *F i vy* are the singleton control actions. The fuzzy rules base of the FYSC is shown in Table 2, where it is constructed using the idea that the yawing sliding index *sψ* can quickly reach zero via a time-consuming trial-and-error tuning procedure. The adopted yaw sliding index reduces the number of fuzzy rules, thus reducing the complexity of adjustment. It shows that the yaw speed error *eψ* will converge to zero when the yaw sliding index converges to zero. The FYSC allows the WBR to rotate along its vertical axis at the desired velocity command. Thus, it can be asymptotically stable under the action of the FYSC.

**Table 2.** The fuzzy rules for the FYSC.

#### *3.3. FRBC Design*