2.1. Mecanum Robot Configuration

This study used a four-wheeled Mecanum carrier capable robot manufactured by Tosa Denshi Company Limited. A direct-current (DC) motor was used to drive the wheels of the robot. The dimensions of the robot are shown in Fig. 1. All units are in millimeters.

The body of the robot was rectangular. All four wheels were mounted parallel to one side of the robot, and the rollers on the wheels were oriented at an angle of 45° to the wheels, as shown in Fig. 1. The robot could move in all directions by adjusting the speed and direction of the wheels, which were controlled independently. The robot was capable of performing various motions, including straight-line movement, 360° rotation, and turning.

The robot received a control signal input every 0.1 s. The control signal consisted of three parts: target wheel velocity \(V\) [rpm], target velocity direction \(\theta\) [°], and target turning speed \(VC\) [rpm]. The resolutions of \(V\) and \(VC\) were 1 rpm, and that of \(\theta\) was 1° ¹⁸.

2.2. Modeling of Mecanum Robot

The robot model consisted of two parts: the kinematic equation between the wheel rotational velocity and the overall motion of the robot and the relationship between the input driving signal and the wheel rotational velocity, which was represented as a difference equation.

The inputs of this model were defined as the target wheel velocity \(V\), target velocity direction \(\theta\), and target turning speed \(VC\). The outputs were defined as the actual outputs of the robot model: wheel velocity \(v\), velocity direction \(\alpha\), and turning speed \(vc\).

The system configuration of the robot is described as follows: the host computer converted the control signal into the target velocity of each wheel and sent it to the four motor drivers that drove the motors. The internal structure of this robot is shown in Fig. 2. In Fig. 2, the abbreviations DAC and ADC stand for digital-to-analog converter and analog-to-digital converter, respectively.

All wheel models were identical, with each model comprising a PID controller, DC motor, reduction gearbox, and wheel.

The command-level relationship between the velocity of the robot and target rotational velocity of each wheel is given by Eq. (1). This equation was provided in the robot manual.

3.1. Wheel Control System Modeling

Based on the analysis of the system components (Section 2), the wheel system was modeled as follows.

Equations (3) and (4) are obtained based on Kirchhoff’s current law and Ampère’s force law ²¹.

Substituting \(i\) in Eq. (4) into Eq. (3) yields Eq. (5), which is a differential equation for the DC motor.

The DC motors and wheels are represented by continuous values, which must be discretized for modeling. Eqs. (6) and (7) were derived using the Taylor series expansion.

Substituting Eqs. (6) and (7) into Eq. (5) yields Eq. (8), which is a discrete equation for the DC motor.

The PID controller, DC motor, and wheel, together with its reduction gearbox, were parameterized as an entire system. This is because the robot is like a black box which cannot be disassembled; therefore, it was not possible to input the motor parameters separately to obtain its response for system identification. Only the input signal was sent to the robot, and the rotation of the wheel was measured as the system output. The discrete formula for the PID controller is as follows:

To derive Eq. (11), the difference form \(\Delta u(k)=u(k)-u(k-1)\) was first derived based on Eqs. (9) and (10), which eliminated the integral term. A similar form of \(\Delta u(k)\) was also derived from Eq. (8). By equating the two expressions of \(\Delta u(k)\) and isolating \(\omega(k+1)\), a difference equation was obtained as follows:

3.2. Gradient Descent Based Parameter Estimation

Parameter estimation was conducted by measuring the rotational velocity of robot wheel \(n(k)\) and then applying the same signal into the wheel model. The model output rotational velocity was denoted as \(\omega(k)\). The average squared error between \(n(k)\) and \(\omega(k)\) at each moment was used as the objective function, and the parameters of the model were estimated to minimize the value of the objective function. That is,

However, depending on the initial value selection, the gradient descent method may potentially converge to a local optimum. Therefore, a method is required to address this issue.

The step response data for parameter estimation cannot cover all dynamic or static characteristics of the robot. Moreover, considering the numerical values of mean squared error (MSE) as the only criterion may not reflect the physical characteristics of the robot.

One way to narrow the region of feasible solutions is to add constraints into the gradient descent process. Introducing the physical or system characteristics of the robot as constraints into the gradient descent process can ensure that the model has the same physical or system characteristics as that of the robot when obtaining the parameters. As the wheel of the robot is a stable system, system stability is a satisfactory constraint for optimizing the gradient descent process.

The Routh–Hurwitz stability criterion ²², which determines system stability by the change in the sign of the coefficients in the Routh arrays, was used as a condition for the stability constraints of the system.

The Routh–Hurwitz stability criterion is only applicable to continuous systems. Therefore, a bilinear transformation was used to map the unit circle stability region of the discrete system to the left half-plane of the continuous system. The model equation was first written in the \(z\)-domain form as follows:

A bilinear transformation was then applied to \(G(z)\). Substituting equation \(z=(1+(\tau/2)\omega)/(1-(\tau/2)\omega)\) into Eq. (13) yields Eq. (14):

According to the Routh–Hurwitz stability criterion, the Routh array was obtained as follows:

Thus, \(b_3\), \(b_2\), \((b_2b_1-b_3b_0)/b_2\), \(b_0>0\) are the stability conditions of the wheel part system.

In summary, the parameter estimation problem can be transformed into the following problem.

In the Routh–Hurwitz stability criterion, the values of the coefficients are restricted by inequality constraints. Therefore, the Karush–Kuhn–Tucker (KKT) condition ²³ can be used to impose inequality constraints on the descent process.

In this problem, the KKT conditions are as follows:

1. (a)

   Feasibility conditions: \(g_i(x)\ge 0\),

2. (b)

   Stationarity conditions:

   \begin{align*} \nabla \mathcal{L}(x,\lambda)=\nabla {H} ({x})+\sum_{i=1}^4 \lambda_i \nabla g_i(x), \end{align*}
3. (c)

   Non-negativity of Lagrange multiplier: \(\lambda_i\ge 0\),

4. (d)

   Complementary slackness condition: \(\lambda_i g_i(x)=0\).

The procedure for solving the target parameters \(x=(A,B,C,D,E,F)\) is as follows:

1. Step 1.

   Initialize the target parameter \(x^0\) and the Lagrange multiplier \(\lambda^0\).

2. Step 2.

   Calculate the gradient of the Lagrange function:

   \begin{align*} \nabla_x \mathcal{L}(x,\lambda) &=\dfrac{\partial \mathcal{L}\left(x^k,\lambda^k\right)}{\partial x}\\ &=\nabla H\left(x^k\right)+\sum_{i=1}^4 \lambda_i^k\nabla g_i\left(x^k\right). \end{align*}
3. Step 3.

   Update the target parameters using the gradient descent method:

   \begin{equation*} x^{k+1}=x^k-\alpha\nabla_x \mathcal{L}\left(x^k,\lambda^k\right). \end{equation*}
4. Step 4.

   Update the Lagrange multipliers:

   \begin{equation*} \lambda_i^{k+1}=\max \left(0,\lambda_i^k+\alpha g_i\left(x^k\right)\right). \end{equation*}
5. Step 5.

   Repeat until convergence. Terminate if the magnitude of the gradient of the Lagrange function is close to zero or the change in \(H(x^k)\) is sufficiently small. Otherwise, return to Step 2 until the convergence condition is satisfied.

3.3. Convergence and Optimality Analysis of the Proposed Method

This section presents the convergence and optimality properties of the proposed stability-constrained parameter estimation method and discusses their applicability.

The objective function defined in Eq. (12) is strictly convex with respect to the parameter \(x\) (\(H(x_k)\)). Without constraints, the gradient descent algorithm converges. The convergence point \(x^{\ast}\), at which \(H(x^{\ast})=H^{\ast}\), corresponds to the global optimal solution.

The gradient of the strictly convex objective function defined in Eq. (12) is a Lipschitz constant \(L_h\). According to the descent lemma ²⁴, the following inequality can be derived.

Whenever the unconstrained update \(x'_k\) falls outside the feasible region, the proposed algorithm projects \(x'_k\) back onto the feasible region using the KKT criterion, yielding the next iteration \(x_{k+1}\).

According to the iterative formula for gradient descent, the following equation can be obtained.

Equation (18) gives the feasible solution region for the problem defined in Eq. (15) with respect to the constraints.

The KKT condition for pulling the outbound point back to the feasible region can be considered as a projection operation, which is mathematically expressed as follows:

The projection theorem is applicable because \(S\) is a closed and convex set (as defined by \(g_1\), \(g_2\), and \(g_3\), as well as by the local convexity of \(g_4\)). According to the projection theorem, the projection operation on the feasible region \(S\) follows the below inequality ²⁵.

The following inequality is obtained by substituting Eq. (17) into Eqs. (16) and (20).

The above inequalities indicate that, even when the feasible region \(S\) is integrated with system stability constraints, the proposed method can still converge to \(H^{\ast}\), as long as the learning rate \(\alpha\) remains within specific thresholds: \(0\le \alpha \le 2/L_h\).

As mentioned previously, \(g_1\), \(g_2\), and \(g_3\) are convex polynomials; therefore, the value of the objective function at the convergence point (\(H(x^{\ast})=H^{\ast}\)) is globally optimal if the location of \(x^{\ast}\) is within \(S\) or on any one of \(g_1\), \(g_2\), and \(g_3\). When \(x\) is on \(g_4\), multiple boundary-point experiments must be conducted on the \(g_4\) boundary to determine whether the function value is the minimum value.

4.1. Parameters Estimation Results

To conduct parameter estimation and evaluate the effects of gradient descent with stability constraints using the Routh–Hurwitz criterion method and gradient descent without stability constraints, the input responses of the target system for four sets of step signals at 10, 40, 80, and 130 rpm were used as modeling data.

Two methods were used to estimate the parameters with the same initial parameters of 10 sets of random values.

Additionally, the gradient descent method was employed with stability constraints using the pole discrimination method for comparison with other methods.

Parameter estimation using the gradient descent method without stability constraints achieved the highest success rate (90%); however, only 30% yielded stable models.

Using the Routh–Hurwitz criterion as a stability constraint, the gradient descent method yielded a success rate of 60% with all successful models being stable.

Introducing the pole discrimination method as the stability constraint led to failure in every trial, rendering the average MSE too high to be calculated.

Figure 3 shows three sets of step inputs and outputs of the target system and successful cases with the same input response of the models obtained using the three methods. The first column shows the robot outputs, and the remaining three columns show the model outputs for the method proposed by Nakagawa et al., the standard gradient descent method, and the proposed method. The first three rows in Fig. 3 show the low-, medium-, and high-speed responses at 10, 80, and 130 rpm, respectively.

Columns (left to right): robot outputs and those obtained using the method proposed by Nakagawa et al., standard gradient descent, and proposed method. Horizontal axis: time [s], vertical axis: rotational speed [rpm].

The settling time of the model obtained using the proposed method was approximately 0.25 s, and the model exhibited an average relative steady-state error of 0.045%. The average relative steady-state error was calculated using Eq. (22).

The optimal value of the objective function \(H(x)\) obtained from the identification experiment was 21.80. The corresponding optimal parameters were as follows:

\(A=0.881\), \(B=-0.151\), \(C=0.122\), \(D=0.0744\), \(E=0.149\), \(F=-0.0748\).

For clarity, the parameter values were rounded off to three significant digits. These parameters were used for the subsequent experiments described in Sections 4.3 and 4.4.

4.2. Comparative Analysis of Method Accuracy

The response data obtained using 13 sets of step signals ranging from 10 to 130 rpm in steps of 10 rpm were used as the test data. For each set, the experiment was repeated five times. For each data group, the following three metrics were used to evaluate the accuracy of the models: relative steady-state error, rise time, and settling time.

The relative steady-state error was calculated using Eq. (22). The rise time was defined as the difference in the time required for the system output to increase from 10% to 90% of the final steady-state value. Settling time is the minimum time required for the system output to enter and remain consistently within a bandwidth of \(\pm 2\)% of the final steady-state value.

The results of the relative steady-state error, rise time, and settling time obtained for the successful models are shown in Fig. 4. The first column shows the values of the three metrics for the robot output, and the next three columns show the three metrics obtained from the model proposed by Nakagawa et al., the standard gradient descent model, and the proposed model.

Columns (left to right): robot output; method of Nakagawa et al.; standard gradient descent method; proposed method. Rows (top to bottom): relative steady-state error [%], rise time [s], and settling time [s]. Horizontal axis: commanded speed [rpm], vertical axis: indicated metric value.

The first row of Fig. 4 compares the relative steady-state errors acquired on the test data.

The relative steady-state errors on robot output show an average of \(4.52\times 10^{-4}\). The model proposed by Nakagawa et al. and the standard gradient descent model exhibit averages of \(2.79\times 10^{-6}\) and 0.0109, respectively. By contrast, the proposed method exhibits an almost zero relative steady-state error of \(-4.41\times 10^{-8}\).

The second row in Fig. 4 shows the rise times. The actual rise time of the robot is approximately 0.11 s. The rise time of the model proposed by Nakagawa et al. and the standard gradient descent model is 0.11 s. The proposed method exhibits a rise time of 0.12 s.

As shown in the third row of Fig. 4, the settling time of the robot is approximately 0.25 s. The settling times of the models presented by Nakagawa et al. and this study are 0.23 s and that of the standard gradient model is 0.22 s. All three methods accurately yielded the step response of the robot.

Nevertheless, the proposed method reduced steady-state error to a near-zero value, while maintaining a dynamic performance comparable to that of the robot.

From left to right: model outputs from Nakagawa et al.’s method, standard gradient descent method, proposed method. Horizontal axis: time [s], vertical axis: rotational speed [rpm].

4.3. Input Tests Beyond Linear Operating Range

The input signal beyond the linear operating range and the response of the target system are shown in Fig. 5. The three response curves of the models for the successful case are shown in the same figure. As observed from the response curve of the target system, the maximum output speed is limited to 135 rpm when the command speed exceeds 135 rpm. This constraint can be incorporated into the model proposed by Nakagawa et al. and the proposed model via programming. The standard gradient descent model exhibits a significant fluctuation in the ascending part of the model response curve, which is also observed for command speeds higher than 110 rpm. This phenomenon leads to faster cumulative errors in response to rapidly changing inputs.

From top to bottom: model outputs from Nakagawa et al.’s method; standard gradient descent method; proposed method. Horizontal axis: time [s], vertical axis: rotational speed [rpm].

4.4. Random Inputs Response

An illustrative five-trial random-input response is shown in Fig. 6. This experiment was conducted to observe the dynamic and steady-state responses of the models. The values of all the inputs were randomly generated. In each experiment, the input design included five cases intended to capture steady-state responses lasting longer than the settling time and five cases intended to capture dynamic responses lasting less than the settling time. Random experiments were conducted five times, with each trial lasting 5.2 s. Table 1 shows the MSE of the model obtained using the three methods in five trials.

As shown in Table 1 and Fig. 6, the average MSE of the model under stability constraints is only 1.077 rpm\(^2\) larger than that of the model without stability constraints. However, the stability constraints ensure the stability of the obtained model without compromising accuracy. Hence, the obtained model is closer to the actual system in theory. This helps in preventing the identification of unstable models, which may produce gradually diverging outputs even under small initial state perturbations and near-zero inputs over a long simulation time. In addition, the MSE of the proposed model is smaller than that of the model proposed by Nakagawa et al.

As evident from the aforementioned results, the gradient descent method incorporating stability constraints using the Routh–Hurwitz criterion can stabilize the model obtained after parameter estimation by ensuring gradient descent accuracy. However, this can lead to a descent failure owing to the addition of constraints. In successful trials, the models satisfied both the stability conditions and exhibited acceptable accuracy performance. By contrast, the step responses of the models obtained in the unsuccessful cases exhibited highly oscillatory responses or nearly flat outputs with singular spikes. These models exhibited extremely large MSE with respect to the actual outputs and lacked physically meaningful behavior, rendering them unsuitable for control design.

The 60% success rate (Section 4.1) was obtained using different initial values. This indicates that these failures appear to be strongly related to the choice of initial values. This is likely because the feasible regions defined by the Routh–Hurwitz criterion may have nonconvex portions, which may induce the gradient descent algorithm to converge to suboptimal regions. These characteristics can hinder the convergence process and cause the descent process to stagnate or diverge even when the input signals are appropriate. This suggests that the success of the method depends not only on constraint formulation but also on proper initialization strategy. Additionally, repeated trials with different initializations can increase the likelihood of obtaining a valid model. From an engineering perspective, a viable approach is to perform multiple runs and select the converged results for deployment. Future studies should consider improved initialization schemes to enhance the convergence rate of this method.

4.5. Analysis of Identified Model in Robot Motion Experiments

To evaluate the applicability of the model, a series of speed–response experiments were conducted on a four-wheeled Mecanum robot platform, and its motion characteristics was measured using an MPU9250 sensor. Specifically, the robot was programmed to start from rest and move in a straight line toward the target velocity or rotate in place during the acceleration, steady-state, and deceleration phases. The target velocities were set at 10, 40, 80, and 130 rpm for linear motion and 1.484 rad/s for rotational motion.

Figure 7 shows the proposed model and the model presented by Nakagawa et al. The figure also shows the measured responses of the robot platform during forward motion experiments at different speeds and stationary rotation experiments. As shown in Fig. 7, the overall trends predicted by the proposed model and that presented by Nakagawa et al. align well with the measured data.

However, during the rising phase, the proposed model and that presented by Nakagawa et al. exhibited rising times of 0.12 s and 0.11 s, respectively, both of which were lower than the rising time of the robot platform (0.19 s). During the steady-state phase, the relative steady-state errors at input speeds of 10, 40, 80, and 130 rpm were 0.186, 0.0689, 0.156, and 0.310, respectively. The relative steady-state error at the input rotational speed of 1.484 rad/s was 0.0306. The model output speeds consistently exceeded the measured values.

The maximum relative steady-state error of 0.31 was observed at the speed of 130 rpm, which can be attributed to the output saturation limits of the robot system, which prevent the actual response from scaling further with higher inputs, as mentioned in Section 4.3. The errors at 10, 40, and 80 rpm likely resulted from the system load, rolling friction, and sensor response latency, which were not explicitly modeled at the robot level.

While the proposed model shows satisfactory performance at the wheel level, deviations were observed between the model outputs and robot-level responses, particularly at the steady-state speed and rise time. These results indicate that the models established under no-load conditions exhibit systematic deviations under actual operating conditions. As the focus of this study was to present a system identification method that is feasible under low command-update frequency conditions, our primary objective was to model the core motor dynamics rather than the full robot behavior.

In future studies, we plan to incorporate load modeling, nonlinear friction compensation, and dynamic sensor error correction to enhance the applicability of the model in real-world control systems. Additionally, by collecting response data under various load conditions and integrating real-time modeling with online update mechanisms, the model parameters can be dynamically optimized. This approach improves control reliability and prediction accuracy in complex environments.

4.6. Computational Efficiency Assessment

The computational efficiencies of the proposed method, gradient descent method, and the method proposed by Nakagawa et al. were compared. All methods used the same programming language, compiler version, operating system, computer, and modeling data. The average computation times and mean memory usages of the three methods are listed in Table 2. The proposed and gradient descent methods were performed 10 times with different initial parameters. The method presented by Nakagawa et al. was used only five times because of its long computation time, and the mean memory usage of this method was sampled during operation.

As shown in Table 2, the computation time and memory usage of the gradient descent methods are markedly shorter than those of the “exhaustive search method” proposed by Nakagawa et al. Owing to the addition of constraints, the average memory usage of the proposed method is higher than that of the gradient descent method without constraints. Constraint incorporation necessitates additional computations and updating of the Lagrange multipliers during each computation. However, reasonable constraints help narrow the region of feasible solutions, thereby enabling the algorithm to locate the optimal solution faster and reducing the number of unnecessary iterations. Therefore, the average computation time of the proposed method was less than that of gradient descent without constraints.