**1. Introduction and Related Works**

Robotics has been an area of research that has attracted the attention of scientists around the world for many decades, but in recent years has enjoyed enormous interest, especially in the safe interaction between humans and robots [1]. The robotic system may be assigned tasks that a person cannot perform in some respects, e.g., if the operations are physically demanding for a person, or if the activity requires high accuracy of movement, or last but not least, the operating and working conditions are not suitable for human work in the vicinity of such a mechanism. All these facts give researchers from different fields many actuations to progress in the research and development of robots and manipulators, whether in the process of designing [2,3], production or the actual application of robots [4–7]. Ref. [4] describes the application of robots, and thus, in particular, the topic of self-employed robots so that its interaction is safe in a dynamic environment. In [5], are described, several types of robots, robot drives, sensors, control methods, as well as industrial applications are described in several sections. Ref. [6], using Lab Activities, helps to understand the parallel between the theory and the application of industrial robotics. Advances in robotics are reflected in every single application area of robotics. This is the case with rehabilitation robotics too, which is described in [7]. Current advances in robotics (with a focus on explaining new concepts) are also described in more detail in [8].

**Citation:** Trojanová, M.; Cakurda, T.; ˇ Hošovský, A.; Krenický, T. Estimation of Grey-Box Dynamic Model of 2-DOF Pneumatic Actuator Robotic Arm Using Gravity Tests. *Appl. Sci.* **2021**, *11*, 4490. https://doi.org/ 10.3390/app11104490

Academic Editor: Angelo Luongo

Received: 23 April 2021 Accepted: 12 May 2021 Published: 14 May 2021

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

**Copyright:** © 2021 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/).

Behind the effective movement of the robot or manipulator are many aspects that have to be mastered in motion control (detailed characteristics of manipulator properties, complex process of motion simulation, correct derivation of kinematic and dynamic models, or last but not least, correct identification of control algorithm). Therefore, the manipulator control will not be correct until the system described in the form of a model is correct. In the publications [9–11], which are manuals in the field of robotics, it is in some chapters to the introduction of kinetic and dynamic methods of modeling such systems. For example, in [9], Chapter 3 is devoted to Inverse Kinematics, and Chapters 4 and 5 are oriented to Newton–Euler Dynamics and Lagrangian Dynamics. Ref. [10], in the introductory parts, describes the basic topics of robotics, and then, among other things, one part is devoted to various areas of application of robotics. Another more comprehensive publication is [11], which explains the constructions, drives, principles of operation of robots and applications of robots, so that the work helps not only researchers but also potential enthusiasts and students of robotics. Publication [12] summarizes the most basic information on the design, but also the control of mechatronic systems, focusing on the needs of engineers and designers who develop and apply their machines to the industrial environment. The part of the work is devoted to robotics, including the dynamic modeling of systems.

Modeling of kinematics and dynamics of manipulators plays an important role in the control accuracy of these devices. For example, in the publication [13], the kinematics and dynamics of the Hyundai HS165 robot are described, which is one of the industrial robots used for spot welding. A large part of the article is devoted to the modeling of kinematics and dynamics, where the load capacity of the manipulator is monitored. Another article [14], in turn, describes a mathematical model of a mobile manipulation robot, where a kinematic model is used for the draft regulation, and a dynamic model is used to monitor (compensate for errors) the speed of the robot. Additionally, the publication [15] describes dynamic modeling, but for a manipulator, which is designed for underwater work and is driven by an electric motor. The manipulator has heterogeneous arms, while each joint of the arm behaves differently within the kinematics. By identifying the minimum number of parameters, the model is experimentally validated on an underwater hybrid remotecontrolled vehicle Ifremer Ariane. In [16] is described the dynamic characteristics of the arm of a manipulator that is driven by pneumatic artificial muscles. In addition to the dynamics of the manipulator, the article also deals with the dynamics of the drive. It is validated the model using measured data from a real device. The created complex model of system dynamics is to be used for control.

Dynamic modeling is generally used to describe the basic relationship between two parameters of the system, namely between force and motion. The dynamic model of the system can be derived within the dynamic modeling based on two basic principles, which are described in the publication [17]. The first is the Euler–Lagrange formulation. This method assesses the system as a whole based on the total kinetic and potential energy (or the difference between these energies). On the contrary, the Newton–Euler formulation describes each member of the system separately. It is necessary to take into account when choosing a method, e.g., number of joints in the manipulator (number of degrees of freedom), location of joints and other mechanical parts, and others. In [18], the dynamic modeling of a robotic manipulator is described using the Euler–Lagrange formulation. The robot has four degrees of freedom, and the location of the base is oriented to the ground/floor (horizontal mounting upwards). In addition to analyzing the properties of this model, they suggest using a PID controller to control the position of the manipulator. In [19], a combination of the Lagrange equation and the Morison formula is used to derive a dynamic model of a manipulator with 3 DOF, which can be used in a hydrostatic environment, while during modeling, they had to take into, for example, other forces acting on the joints when working with the manipulator underwater. The next publication [20] used the Euler–Lagrange method to derive dynamic equations for the prototypes of a four-legged robot—Mini-Bot—where they tracked the position and angular velocity of each leg in the Cartesian coordinate system. They used the Matlab-Mupad software tool as

a tool for deriving dynamic equations. The Newton–Euler formulation was used in the publication [21] to derive object equations, while a recursive Gibbs–Appell formulation was used to obtain the dynamic equations of a flexible, cooperative manipulator. The manipulator consists of two arms, which have two flexible links. Two limitations had to be taken into account in the simulations, namely dynamic and kinematic limitations resulting from the placement of the object on a non-holonomic mobile platform. In the article [22], the research of the derivation of a dynamic model using a Newton–Euler formulation for a rigid manipulator is described, which consists of two links.

Some papers also present new methods in their publications that can be used for the dynamic modeling of manipulators. For example, [23] proposes a new method for derivate a dynamic model of a manipulator with 6 DOF. The hysteresis effect of friction in the joints is described using a centrosymmetric static friction model. The dynamic parameters of the six joints were identified using a hybrid whale optimization algorithm and a genetic algorithm. One of the most recent publications [24], contains another new approach to deriving equations of motion, called Pseudo-Symbolic Dynamic Modeling. It is a numerical algorithm whose task is to create functions that form the basis for the inverse dynamics of the manipulator. The result is simplified dynamic models, while this algorithm has been implemented in Matlab and can be applied to manipulators up to 7 degrees of freedom of movement.

To achieve better results of the modeling, models can be optimized using numerical iteration methods, for example, algorithms as, e.g., the Levenberg–Marquardt method or the Trust Region Algorithm. In [25], the Levenberg–Marquardt method for the unconstrained optimization problem is studied, whereas the method of convergence technique (using line search) is proposed. In one of the publication's chapters [26], the Levenberg– Marquardt method is used to determine the minimum volume of the longitudinal cooling fin. The Trust Region Algorithm is also used to solve nonlinear problems. At present, its popularity is higher than with the previously mentioned method. Many researchers describe algorithms that are derived from the Trust Region Algorithm. The aim is to increase its reliability and convergence properties. For example, in [27], two trust-region algorithms are proposed that examine trust outside the specified area. In [28] is described the application of the algorithm to the online estimation of tire/road friction parameters in real-time. In addition, a new LuGre model-based nonlinear least-squares parameter estimation algorithm is presented in this work.

The aim of the research described in the article is to use an optimization method to estimate parameters of the model (improve the results of validation) and the comparison of the model with the parameters that come from CAD and its improvement by estimating the parameters based on gravitational measurements. The article describes the dynamics of a manipulator with two degrees of freedom, while the dynamic model of the manipulator's arm is derived using Lagrangian formalism. The compiled model was implemented in Matlab, taking into account the physical parameters of the manipulator. A scheme (model) was compiled in the Simulink, which was used for the subsequent validation process. The article is divided some parts, where the introduction is followed by the dynamic model of the kinematic structure of a planar arm. Next follows the characteristics of the experimental manipulator. Data measured from the system using gravity test are described in chapter Measured Data. A simulation model was developed for the validation process, and the results of the validation with CAD parameters are described in the fifth chapter. To achieve better results were estimated the parameters of the model. The sixth chapter described the optimization algorithm used to estimate parameters of the model (Trust Region Algorithm for Nonlinear Least Squares), and in the next chapter are summarized results of validation with estimated parameters. The last part consists of a discussion, in which the results of validation and estimation of parameters, based on statistical indicators MAE and NRMSE.

#### **2. Dynamic Model of the System**

The experimental system represents one of the most basic kinematic structures (a planar arm with two rotating joints moving in the x-y plane). The schematic diagram of the kinematic structure of the manipulator shown in Figure 1 serves to show the basic elements, their location within the whole system, but also the basic parameters. Two rigid arms, lengths *l*<sup>1</sup> and *l*2, have a mass of *m*<sup>1</sup> and *m*2. The distance from the axis of rotation of the joint to the center of gravity of the arm is denoted as *c*<sup>1</sup> and *c*2. *I*<sup>1</sup> and *I*<sup>2</sup> represent the moments of inertia of the bodies (arms) concerning an axis that passes through the center of gravity of the body and is parallel to the horizontal axis *x*. The deflection of the arm from the zero position (the angle of rotation of the joint is not equal to 0◦) is indicated by the angles *ϕ*1, *ϕ*2, which represent the degrees of freedom of the manipulator. The deflection of the arm in the positive direction from the vertical *y*-axis is expressed by the angle of rotation of the joint *ϕ*1. The angle *ϕ*<sup>2</sup> indicates the angle of rotation of the second joint, but concerning the axis of the first deflected link.

**Figure 1.** Scheme of kinematic structure of experimental system.

Lagrange's equations were used to derivate a dynamic model of the described system, which has two degrees of freedom [29]. The general formulation for a robot with 2 degrees of freedom is given in Equation (1), where *L* is Lagrangian, *ϕ<sup>i</sup>* is joint positions, derivation of the joint position is joint velocities, and *τ<sup>i</sup>* is force term.

$$\frac{d}{dt}\left[\frac{\partial L}{\partial \dot{q}\_i}\right] - \frac{\partial L}{\partial q\_i} = \tau\_i \qquad \qquad i = 1, 2\tag{1}$$

The Lagrangian is given by the difference of kinetic and potential energy, and therefore Equation (2) applies to its determination. Member *Ek* (*ϕ*, . *ϕ*) expresses the kinetic energy of the manipulator and is formed by the sum of individual kinetic energies of members 1 and 2 (*Ek*<sup>1</sup> and *Ek*2). Member *Ep (ϕ)* characterizes the energy of the manipulator concerning its position and is given by the sum of individual potential energies of arms 1 and 2 (*Ep*<sup>1</sup> and *Ep*2).

$$L\left(\boldsymbol{\varrho}, \dot{\boldsymbol{\varrho}}\right) = E\_k\left(\boldsymbol{\varrho}, \dot{\boldsymbol{\varrho}}\right) - E\_{\mathcal{P}}\left(\boldsymbol{\varrho}\right) \tag{2}$$

#### *Determination Lagrange's Equations*

By calculating the Lagrangian and applying Equation (1), two nonlinear differential equations are obtained, which represent the equations of the manipulator dynamics. For link 1, the force component *τ*<sup>1</sup> is determined in the final form:

$$\begin{array}{l} \pi\_{1} = \begin{array}{l} \left[ m\_{1} \cdot c\_{1} \right]^{2} + m\_{2} \cdot c\_{2} \, ^{2} + m\_{2} \cdot l\_{1} ^{2} + 2m\_{2} \cdot c\_{2} \cdot l\_{1} \cdot \cos \, \varphi\_{2} + \, l\_{1} + l\_{2} \right] \cdot \ddot{\boldsymbol{\varphi}}\_{1} \\ + \left[ m\_{2} \cdot c\_{2} \, ^{2} + m\_{2} \cdot c\_{2} \cdot l\_{1} \cdot \cos \, \varphi\_{2} + \, l\_{2} \right] \cdot \ddot{\boldsymbol{\varphi}}\_{2} - 2m\_{2} \cdot c\_{2} \cdot l\_{1} \cdot \sin \, \varphi\_{2} \cdot \dot{\boldsymbol{\varphi}}\_{1} \cdot \dot{\boldsymbol{\varphi}}\_{2} \\ - m\_{2} \cdot c\_{2} \cdot l\_{1} \cdot \sin \, \varphi\_{2} \cdot \dot{\boldsymbol{\varphi}}\_{2}^{2} + m\_{2} \cdot \boldsymbol{g} \cdot c\_{2} \cdot \sin \, (\varphi\_{1} + \varphi\_{2}) \\ + \left[ m\_{1} \cdot c\_{1} + m\_{2} \cdot l\_{1} \right] \cdot \boldsymbol{\varphi} \cdot \sin (\varphi\_{1}) \end{array} \tag{3}$$

and for link 2, the force component *τ*<sup>2</sup> is defined:

$$\begin{aligned} \tau\_2 &= \left[ m\_2 \cdot c\_2^2 + m\_2 \cdot c\_2 \cdot l\_1 \cdot \cos \,\varphi\_2 + l\_2 \right] \cdot \ddot{\varphi}\_1 + \left[ m\_2 \cdot c\_2^2 + l\_2 \right] \cdot \ddot{\varphi}\_2 \\ &+ m\_2 \cdot c\_2 \cdot l\_1 \cdot \sin \,\varphi\_2 \cdot \dot{\varphi}\_1^2 + m\_2 \cdot \dot{\varphi} \cdot c\_2 \cdot \sin \left( \varphi\_1 + \varphi\_2 \right) \end{aligned} \tag{4}$$

The Lagrangian dynamic model expressed by Equations (3) and (4) can be written in a compact form, and when extended by the friction component *F* then has the form of Equation (5).

$$M(\varphi)\ddot{\varphi} + \mathcal{C}\left(\varphi, \dot{\varphi}\right)\dot{\varphi} + G(\varphi) + F(\dot{\varphi}) = \pi \tag{5}$$

The friction component *F* can be defined by Equation (6), where *ω* represents the angular velocity in [rad/s]; *Tc* is the Coulomb friction magnitude is given in [N·m]; the Viscous friction coefficient *β* has given in units [N·m·s/rad]. This whole expression on the right side of the equation is also called Coulomb and Viscous friction [30].

$$F = \text{sgn}(\omega) \cdot T\mathbf{c} + \beta \cdot \omega \tag{6}$$

#### **3. Experimental System—Planar Robotic Arm**

The object of the research was a planar robotic arm of the manipulator, which is driven by an unconventional type of actuators (pneumatic artificial muscles) and has two degrees of freedom. The manipulator is shown in Figure 2. The construction of the manipulator consists of metal profiles, which ensure the stability of the manipulator and the attachment of the upper arm of the manipulator to the upper base of the structure. A load weighing 3.34 kg is attached to the other arm of the manipulator. The type of pneumatic drive was chosen for the manipulator drive-fluid muscles from the manufacturer FESTO. If these muscles are pressurized, their activation occurs, as a muscle contraction, which causes the translational tensile force of the muscle. However, since the muscles within the manipulator are antagonistically involved, and also the chain-sprocket and joint subsystems are present within the system, the resulting movement is rotational (the subsystems ensure the transmission of torque). During the research, the angle of rotation of the joint was monitored as an output parameter. An incremental sensor was installed in the manipulator to be able to read this parameter. To achieve the required outputs with the manipulator, the pressure in the muscles must be regulated—electronic pressure regulators (one for each muscle) are used for this. In addition to the possibility of pressure coordination, the regulators contain built-in pressure sensors. An AMD Athlon X2 computer with a 5 GB RAM processor was used as a control unit to control the manipulator. The software support of the stand consisted of the use of the Matlab and Matlab Simulink programs, while the control signals were operated with a Humusoft MF624 card.

**Figure 2.** Experimental system—planar robotic arm with 2-DOF and actuated with PAMs.

Table 1 contains an overview of the main components used in the experimental system. In addition to the name of the component, its basic parameters are also given. A total of 4 fluid muscles of the same size and diameter were used. Each muscle had one electronic pressure regulator (a total of four were used). An incremental sensor was used to read the output variable. The mechanical parts of the stand (except the structure) consist of two arms, two joints, two chains and two sprockets with a diameter of 0.07 m.


#### **4. Measured Data**

Measurements were performed by using gravity experimental tests, the task of which was to obtain the waveforms of the output monitored parameter—the joint angle of the manipulator's arm. The first measurement started at a joint rotation angle of 10◦, i.e., the arm was moved from the zero position by 10◦. Then, it was released, and the values of the joint rotation angle were recorded by incremental sensors until the moment when the arm stopped. Thus, measurements were subsequently performed for the joint angle 20◦, 30◦ and 40◦ too. The total number of measurements for one initial angle was 20, with the tests being performed in both directions. During the processing, all waveforms were converted to one side based on a simplified assumption of symmetry of the arm dynamics. The tests were performed independently for link 1 and link 2. The total measurement time was different for each initial joint angle, but the sampling was the same at 0.005 s.

#### *4.1. Measured Data—Link 1*

Figure 3 shows 4 boxplots (when viewed from left to right for an initial joint angle of 10◦–40◦), which show how the initial values of the joint angle have changed in 4 different measurement cases. Each boxplot was processed based on the measured initial joint angles in 20 measurements. It can be seen from the graphs that the mean of the 20 measured initial values (red line) approaches the required input at most in the case of an angle of 20◦ (the second graph from the left with mean value 20.016◦), where it is not possible to see the 25th and 75th percentiles (blue rectangle) is present. Only one different value is also a remote value of 19.872◦ (shown in the second graph on the left by a red cross). It is possible also to observe a remote measured value of 40.464◦ in the fourth graph from the left when the initial angle was set to the reference value of 40◦. The largest difference between the minimum and maximum set input value (black horizontal line) can be observed for an initial joint angle of 10◦ (first graph from the left).

**Figure 3.** Boxplots of initial values of joint angles—link 1.

From the 20 measured time courses of the joint angle, where the initial joint angle is 10◦, the angles were set for individual measurements in the interval <9.216◦; 10.08◦>. The maximum value of the interval represents the statistic parameter Mode. Figure 4 represents the mean response (red curve) and standard deviation (the grey area) from measured data.

**Figure 4.** The mean response and standard deviation for the joint angle 10◦—link 1.

The interval of the initial angles of rotation of the joint with a reference value of 20◦ was <19.872◦; 20.016◦>. Such a small range was caused by the fact that out of 20 measurements, the value 20.016◦ was set 19 times. Waveforms of standard deviation and mean response of data are shown in Figure 5.

**Figure 5.** The mean response and standard deviation for the joint angle 20◦—link 1.

Figure 6 shows the waveform of mean response at a joint angle of 30◦ for link 1, where the initial values of the joint angle were set in the range <29.808◦; 30.24◦>. More than half of the measurements started from 30.096◦. Figure 6 shows a more pronounced standard deviation also (compared to the deviation for angles of 10 ◦ and 20 ◦).

**Figure 6.** The mean response and standard deviation for the joint angle 30◦—link 1.

The last group within link 1 is closed by the curve mean response and standard deviation if the initial angle was 40◦, as shown in Figure 7 (respectively, in reality, the initial angles were set in the interval <39.888◦; 40.464◦> during the measurements). Half of the measurements were started at a joint rotation angle of 40.032◦.

**Figure 7.** The mean response and standard deviation for the joint angle 40◦—link 1.

#### *4.2. Measured Data—Link 2*

Figure 8 shows four boxplots for link 2 (when viewed from left to right for an initial joint angle of 10◦ to 40◦). In these graphs, in comparison with the graphs for link 1, it can be seen that the minimum and maximum values of the set initial angles are much further apart, which results in a more significant shift of the mean value from the desired initial angle (especially for the angle −40◦, fourth graph from the left). The remote value was measured only for an angle of 20◦ (red crosses, second graph from the left), two times (−17.28◦ and −14.832◦).

**Figure 8.** The boxplots of initial values of joint angles—link 2.

While for link 1, the initial angles were set in the range of up to 1◦, the setting of the initial joint angle for the measurements of link 2 was more difficult in terms of constraints in the construction of the system. Figure 9 graphically presents the mean response (red curve) and standard deviation (grey area) when the initial angle was set to −10◦. The total interval of the initial values of the angle of rotation was <−7.92◦; −10.08◦>, the Mode from the twenty initially set angles was −8.928◦.

**Figure 9.** The mean response and standard deviation for the joint angle 10◦—link 2.

As the joint angle for link 2 increases, the range of the intervals of the initial values of the joint angle also increases (it is also visible in the area of the standard deviation parameter in Figure 10). The interval for the initial value of the joint rotation angle −20◦ is <−14.832◦; −21.888◦>, where the difference between the minimum and maximum value of the interval already exceeds 7◦. Nevertheless, the value of the Mode parameter (−20.736◦) is close to the reference value.

**Figure 10.** The mean response and standard deviation for the joint angle 20◦—link 2.

The interval for the reference value −30◦ for link 2 was <−23.904◦; −31.248◦>, where from the range of initial values, it is possible to monitor a more significant value of the parameter Mode (−24.048◦) in comparison with the reference value. The shift of the beginning of the mean response curve (Figure 11) and the standard deviation area itself is caused by one of the measurements when the initially set value was −23.904 ◦.

**Figure 11.** The mean response and standard deviation for the joint angle 30◦—link 2.

The last measured data was for the initial joint angle was −40◦. The total interval of the initial values of the joint angle was <−32.976◦; −42.192◦>, the Mode of these twenty values was −38.736◦. In Figure 12, it is possible to see the fact that the variance in the measurements is the most expressive, which of course, was reflected in the standard deviation.

**Figure 12.** The mean response and standard deviation for the joint angle 40◦—link 2.

Gravity tests showed that within link 1, the time course of mean response of the joint angle was more uniform with gradual stabilization to the value of the zero joint angle. For link 2, it is clear from the graphic waveforms that the course is different compared to link 1. A more pronounced stabilization for link 2 for each initial joint angle begins only after a certain amplitude has been exceeded. Additionally, in the case of graphical dependencies, link 2 can see more significant friction, which can affect the validation process. Further, from the sample mean responses and standard deviations for link 1 and link 2, it can be presented that for initial joint angles of 10◦ and 20◦, the standard deviation range is narrower than for 30◦ and 40◦ for both axes, which means that at higher initial joint angles, the deviation from the mean response is more pronounced. This could be caused by the stronger manifestation of unmodeled effects, such as additional oscillations of the structure, which was not tightly fastened to the floor (itself slightly uneven). In addition, the larger dispersion of the initial angles had an impact on the wider standard deviation band for link 2 compared to link 1. At the same time, it can be seen from the mean response that the arm stabilization was not completed at 0◦ in any case due to the effect of friction being more pronounced for this axis. The most visible deviation is for a joint rotation of 20◦ and 30◦ for link 1 (Figures 5 and 6). After the phase of implementation of gravity tests and analysis of measured data, the stage of validation in the Matlab and Simulink followed.

In research were used only experiments (not simulation) obtained through gravity tests as a validating agent. However, in the first case, are used parameters values that could be obtained from CAD, and thus could directly be used for simulation without the need for experiments (nevertheless, even in this case, this was validated using the measurements—Chapter 5). In the second case, these parameter values were used as initial values for parameter optimization to obtain an improved model that includes the effects not incorporated in the model with CAD parameters—Chapter 7 (such as fluidic muscles, chains, etc.).

#### **5. Validation**

Based on the dynamic model presented in the section Dynamic Model of The System, a simulation scheme in the Simulink was created for the validation process. The basic physical parameters of the system (mass, length, a moment of inertia and distance from the axis of rotation of the joint to the center of gravity for link 1 and link 2) were exported from the CAD software, which was used to make a 3D model of the system. The coefficients viscous friction coefficient and Coulomb friction magnitude were determined experimentally. An overview of these parameters (with parameter designation, value and unit) is summarized in Table 2.


**Table 2.** The list of input parameters of the Simulink model.

The validation outputs were evaluated based on the MAE (Mean Absolute Error) and NRMSE (Normal Root Mean Square Error) statistical indicators. The indicators are used to compare the output of the model with the outputs of the system. The MAE indicator by a numerical value represents the error that occurred during the validation compared to the measured data. The general mathematical formulation of MAE is given in Equation (7), where yj represents the output of the system on the k-th sample and у¯<sup>j</sup> represents the output of the model on the k-th sample. The NRMSE, unlike the MAE, compares the degree of agreement of the simulated output with the measured data (Equation (8)). The desired value is for NRMSE = 100% when the model output corresponds 100% to the system output.

$$MAE = \frac{1}{n} \sum\_{j=1}^{n} |y\_j - \hat{y}\_j| \tag{7}$$

$$F\_{\mathcal{G}of} = \left( 1 - \frac{\sqrt{\sum\_{j=1}^{n} \left[ y\_j - \mathcal{G}\_j \right]^2}}{\sqrt{\sum\_{j=1}^{n} \left[ y\_j - \frac{1}{n} \sum\_{j=1}^{n} y\_j \right]^2}} \right) \times 100\% \tag{8}$$

#### *5.1. Results of Validation—Link 1*

The validation results based on the selected MAE and NRMSE criteria for link 1 are summarized in Table 3. The model output for link 1 achieved the best result based on both statistical indicators at an initial joint angle of 10◦ (MAE = 1.0932 and NRMSE = 35.47%). However, even with the best model, the results represent that the agreement of the model output with the measured data is low and insufficient.


**Table 3.** Results of validation for link 1.

Figures 13–16 show the model outputs for the joint 10◦–40◦ for link 1. It can be seen from the waveforms that the model outputs (blue solid line) are largely unable to simulate an expected output represent by mean responses (red dashed line) at the specified input parameters. The error (black dotted line) for a joint angle of 10◦ has a decreasing tendency with a gradual approach to zero (Figure 13). The error waveforms for the joint angles 20–40◦ (Figures 14–16) are initially smaller but then increase over time, and at the moment when the mean response and thus the output model approaches zero, the value of the error is minimized.

**Figure 13.** Model output of validation for the joint angle 10◦—link 1.

**Figure 14.** Model output of validation for the joint angle 20◦—link 1.

**Figure 15.** Model output of validation for the joint angle 30◦—link 1.

**Figure 16.** Model output of validation for the joint angle 40◦—link 1.

#### *5.2. Results of Validation—Link 2*

From the analysis of the validation results for link 2 (Table 4) based on the above indicators, it is clear that the values of the indicators achieve even worse results than in the case of link 1 (this may be due to oscillations that originate from the interconnection of individual parts of the system for link 2 as for link 1). The highest value of the indicator NRMSE = 18.14% was achieved for the model when the input joint angle was 30◦, but MAE achieved the best result at the initial joint angle of 10◦ (MAE = 0.5655).

**Table 4.** Results of validation for link 2.


Since the course of the measured data for link 2 has a different character compared to the courses of link 1, and at the same time, a larger oscillation is present, these facts are also reflected in the actual outputs of the model. In Figures 17–20, where the model outputs of validation for the joint angle 10–40◦ for link 2 are presented, it is possible to see to what extent over time the model output (blue solid line) was able to approach an average of the measured values (red dashed line). Since the output of the model moves many times close to the zero value of the joint angle, in all cases, the shape of the error (black dotted line) is similar to the course of the mean response. Similarly, for both link 1 and link 2, at the moment when the mean response and thus also the output model approaches zero, the error value is minimized.

**Figure 17.** Model output of validation for the joint angle 10◦—link 2.

**Figure 18.** Model output of validation for the joint angle 20◦—link 2.

**Figure 19.** Model output of validation for the joint angle 30◦—link 2.

**Figure 20.** Model output of validation for the joint angle 40◦—link 2.

Based on the insufficient results obtained in the validation phase and presented either in the form of MAE and NRMSE criteria parameters or by time dependences of the model outputs alone, compared to the average of measured values and the error between model output and real setpoint, was as next step selected estimation of model parameters. Parameter estimation is an efficient tool that uses mathematical models to accurately describe the behavior of systems. This phase is described in the following two sections.

#### **6. Parameter Estimation of 2-DOF Planar Arm**

To estimate the parameters of the model given in Equations (9) and (10), one has to formulate this as an optimization problem. Based on [31], the parameter estimation problem can be formulated as follows: assuming it is a vector of measured values of joint angle *<sup>λ</sup>* ∈ *n<sup>λ</sup>* consisting of K samples, a vector of model parameters *<sup>θ</sup>* ∈ *n<sup>θ</sup>* has to be found using an optimization algorithm and properly defined objective function. The objective function used in the experiments was a sum-squared type of error function defined as:

$$\theta(\lambda) = \sum\_{k=1}^{K} \left[ \overline{\theta}(k) - \theta(k, \lambda) \right]^T \left[ \overline{\theta}(k) - \theta(k, \lambda) \right] \tag{9}$$

The purpose of the optimization algorithm is to find the values of parameter vector *<sup>λ</sup>*<sup>ˆ</sup> ∈ *n<sup>λ</sup>* so that

$$
\hat{\lambda} \in \arg\min\_{\lambda \in L}(\lambda) \tag{10}
$$

where *<sup>L</sup>* ⊂ *n<sup>λ</sup>* is the parameter space.

#### *Trust Region Algorithm for Nonlinear Least Squares*

Nonlinear optimization is the minimization of a function according to the following formulation:

$$\min\_{\mathbf{x}\in\mathbb{R}^n} f(\mathbf{x})\tag{11}$$

subject to *ci*(*x*) = 0 for *i* = 1,2, ... ,*me* and *ci*(*x*) ≥ 0 for *i* = *me* + 1, ... ,*m*. Functions *f*(*x*) and *ci*(*x*) hold that they are real, defined in *Rn*; at least one of them is nonlinear; for the coefficients *m* and *me*, it holds that are nonnegative integers, and *m* ≥ *me*. It holds that, if *m* = *me* = 0, the function (Equation (11)) presents an unconstrained optimization problem. If this condition does not apply, then Equation (11) is a constrained optimization problem.

Optimization problems are solved by numerical iteration methods, which means that the approximate solution of the function at the point *xk* is in the k-th iteration. The new point (approach to the solution) *xk+*<sup>1</sup> is determined computationally based on an algorithm. Such approximation to the solution is applied until the solution is acceptable. An important aspect in the direction of the search, whereby solving subtasks it is verified whether the points are approaching, resp. do not come close to the solution. One of the algorithms that can be applied is the Trust Region Algorithm [32–35]. The Trust Region Algorithm for Nonlinear Least Squares is similar to the method Levenberg–Marquardt [25,36] used to solve nonlinear functions. Although it is a newer algorithm (compared to Levenberg– Marquardt), it brings advantages such as robustness, reliability and convergence properties. In this algorithm, the model is said to be trusted if the area is close to the solution. Such an area is called a trust region and is modified at each subsequent iteration. This means that if in a certain iteration the approximate model is not good enough (it, therefore, contains a bad point), the area of confidence will not be increased but, on the contrary, will decrease.

After a nonlinear optimization problem was defined using Equation (12), if this problem were unconstrained, the algorithm model could be constructed under the condition that the trial step in the k-th iteration would be determined as [34]:

$$\min\_{d \in \mathbb{R}^n} g\_k^T d + \frac{1}{2} d^T b\_k d = \varrho\_k(d) \tag{12}$$

Subject to *d*<sup>2</sup> ≤ Δ*<sup>k</sup>* (13)

where *bk* is a symmetric matrix of dimension *n*x*n* approaching the Hessian of the function *f*(*x*); *gk* represents the gradient at the current iteration *xk*; Δ*<sup>k</sup>* is greater than zero and represents the trust-region radius.

If Equation (13) is subtracted from Equation (12), then the difference is *Sk*. Subsequently, if the difference *ϕk*(*0*)−*ϕk*(*Sk*) is determined, the result is *Pk*, which represents the predicted reduction. For *Pk*, holds it is positive if *xk* is a stationary point and *bk* is a positive semi-definite. Actual reduction *Ak* (reduction in objective function) can be determined as Algorithm 1 [34]:

$$A\_k = f(\mathbf{x}\_k) - f(\mathbf{x}\_k + \mathbf{S}\_k) \tag{14}$$

If the ratio is the actual reduction *Ak* and the predicted reduction *Pk*, it is obtained an important ratio *rk*, which firstly expresses whether the test step is acceptable and secondly determines the adaptation of the new trust-region radius.


*τ<sup>i</sup>* is chosen by the user (usually *τ<sup>0</sup>* = 0, *τ*<sup>1</sup> = 2, *τ*<sup>2</sup> *= τ<sup>3</sup>* = 0.25, *τ<sup>4</sup>* = 0.5). Normally, τ<sup>0</sup> is chosen to be zero, but if *τ<sup>0</sup>* >0, count of lim *k*→∞ *inf gk*<sup>2</sup> = 0 will be stronger. As a satisfactory optimal criterion for terminating the algorithm at point *xk* is *gk*<sup>2</sup> ≤ .

#### **7. Parameter Estimation of Simulink Model**

Within the estimation of model parameters, attention, in this case, was focused on the parameter estimation of the Simulink model, where the task is to process the input/output data used for the testing phase. In the validation phase, tests were performed separately for link 1 and link 2 and at the same time for each axis at four different initial joint angles. However, when estimating the model parameters, the estimation takes place simultaneously for link 1 and link 2 for the specified joint angle. To create the experiment, it was necessary to select the measured outputs—the matrix of the required measured outputs depending on time. Specifically, in this case, it was the time dependence of the mean response obtained from the measured data. Another important step in configuring the experiment was to choose which parameters to estimate. In the performed experiment, all physical parameters of the system were selected for estimation, which are listed in Table 5. At the same time, it was necessary to determine the limits of the set parameter. The limits were set to ±30% of the original parameter values obtained from the CAD software (excluding friction). The basis of parameter estimation is an optimization method and algorithm that determine the sequence of steps in which the estimation takes place and what criteria for parameter estimation are chosen. Based on what was mentioned in Section 4, the Trust Region Algorithm for Nonlinear Least Squares was chosen.


**Table 5.** The list of the number of samples and the parameters' spaces.

The result of parameter estimation is new estimated values of selected model parameters, which were reused for validation to obtain better model results. An overview of these parameters is summarized in the list of model parameters after parameter estimation in Table 6. The simulation scheme in Simulink was modified with new detected parameters, and the process of output validation based on measured data was performed as in section Validation. The results of the validation after estimating the parameters are divided into two parts—for link 1 and link 2.

**Table 6.** The list of model parameters after parameter estimation.


#### *7.1. Results of Validation after Parameter Estimation—Link 1*

The results of validation after the estimation of parameters were assessed based on the same statistical criteria as in the validation before estimation. Table 7 provides an overview of the MAE and NRMSE indicators for link 1. From the results, it can be seen that based on the MAE indicator, the model at the initial joint angle of 10◦ (MAE = 0.6076) is the best, while based on the NRMSE indicator, it is the model at an initial joint angle of 40◦, with NRMSE = 84.62%.

**Table 7.** Results of validation after parameter estimation for link 1.


The improvement of the model results is obvious also from the graphic curves. Figures 21–24 show the model outputs after parameter estimation for the joint 10◦–40◦ for link 1. The model output (blue solid line) approaches the required mean response (red dashed line) with increasing initial joint angle, which also indicates an increase in the value of the NRMSE indicator in Table 7. Before estimating the parameters, the model output had the problem of achieving mean response amplitudes and also the peaks of the model output amplitudes that were shifted for the time axis.

**Figure 21.** Model output after parameter estimation for the joint angle 10◦—link 1.

**Figure 22.** Model output after parameter estimation for the joint angle 20◦—link 1.

**Figure 23.** Model output after parameter estimation for the joint angle 30◦—link 1.

**Figure 24.** Model output after parameter estimation for the joint angle 40◦—link 1.

#### *7.2. Results of Validation after Parameter Estimation—Link 2*

The results of the validation before estimating the parameters were significantly worse for link 2 compared to the results for link 1. However, after estimating the parameters and performing the validation again, the models were able to produce much better model outputs, as evidenced by the MAE and NRMSE values in Table 8. The MAE indicator had a value of less than one for each initial joint angle, with the best output of the model obtained at an initial joint angle of 10◦ (MAE = 0.2407). The value of the NRMSE indicator was around 80% at all initial joint angles examined, with the model achieving the best result at an angle of 20◦ (NRMSE = 82.82%).

**Table 8.** Results of validation after parameter estimation for link 2.


This significant difference in the improvement of the model results after parameter estimation is especially visible in Figures 25–28. The model before parameter estimation had the biggest problem for link 2 with data in the time zone 0 to 2.5 s when there are significant fluctuations of the manipulator's arm. However, these jumps are followed by more pronounced damping. However, the model after estimation can simulate the output (blue solid line) at this time with high accuracy to the mean response (red dashed line) for each initial joint angle, which is of course also reflected by the change in the curve error (black dotted line).

**Figure 25.** Model output after parameter estimation for the joint angle 10◦—link 2.

**Figure 26.** Model output after parameter estimation for the joint angle 20◦—link 2.

**Figure 27.** Model output after parameter estimation for the joint angle 30◦—link 2.

**Figure 28.** Model output after parameter estimation for the joint angle 40◦—link 2.

#### **8. Discussion**

The estimation of model parameters using the Trust Region Algorithm for Nonlinear Least Squares optimization method for the known structure of the planar kinematic system model derived based on Lagrange formalism was intended to improve the validation results. As mentioned above, the limits of the estimated parameters were set to ±30% of the original value of each parameter (excluding friction). This value was set because deviations could occur within the parameters specified by the CAD software, e.g., the resulting weights can be skewed, as an element, such as a chain or muscles, is not included in the CAD. CAD, therefore, determined the parameters only based on the modeled construction. A summary and comparison of validation results before and after parameter estimation is given in Table 9 for link 1 and in Table 10 for link 2. It is clear from the indicators that the estimation of parameters significantly helped to achieve better model

outputs in the process validation. If the indicators for link 1 are compared, the most significant improvement occurred at the initial joint angle 40◦, when the original output of the model reached MAE = 6.1717 and NRMSE = 16.69%, while after estimation, the output of the model was obtained at MAE = 1.1045 and NRMSE = 84.62%. For link 2, the results in the initial validation were even worse than for link 1. As mentioned above, the cause could be a more pronounced oscillation in this axis. Therefore, this fact was expected to play a role in the results even after optimization. However, it is clear from the results that the optimization within the parameter estimation helped the model so much that the outputs not only achieved much better results than before the estimation, but within the NRMSE indicator, they were also similar in value (stable) at different initial rotation angles. However, the most significant change after the determination occurred at an initial angle of 40◦ (taking into account the values of both indicators).


**Table 9.** Comparison of results validation before and after parameter estimation for link 1.



Since the process of estimating the parameters took place simultaneously for both axis 1 and axis 2, from the results of the estimation, it was necessary to select the model that achieved the best results as a whole (not separately for axis 1 and axis 2). As already mentioned at the beginning of this chapter, the most significant progress in the results occurred at the initial angle of rotation of the joint of 40◦, so this model was chosen as the best. Figure 29 shows the time courses of the joint angle for both axes at the initial joint angle of 40◦ before the model optimization process. For both axes, you can see that the output model (blue solid curve) is different compared to the required mean response (red dashed curve). After optimizing the model, however, significant differences (improvements) for both axes can be seen in the model outputs in Figure 30. The output of the model no outrun the required output (this phenomenon was visible especially for axis 1 before estimation). The model also achieves the amplitude of the measured values much better (especially with axis 2 before estimation, this problem is visible).

**Figure 29.** The best model output before parameter estimation for the joint angle 40◦.

**Figure 30.** The best model output after parameter estimation for the joint angle 40◦.

Since the essence of the optimization of the model was the estimation of the model parameters using the chosen algorithm, it is important to analyze the estimated parameters themselves, which ensured the improvement of the results. These were the parameters of the manipulator together with the coefficients found in the friction components. In connection with this, it should be noted that since the goal was set in the validations—the best possible model outputs and results within the indicators. Therefore the arm parameters, such as length, weight, etc., were also included in the estimation. Table 11 shows a comparison of the values of these parameters before estimation and after estimation.


**Table 11.** Comparison of parameters before and after parameter estimation.

The weights of the arm parts *m*<sup>1</sup> and *m*<sup>2</sup> were changed by the estimation—the weight of the upper part of the arm was estimated to 0.9 kg extra, while the weight of the lower part of the arm was estimated to be 0.076 kg less. The parameter, Distance to the center of mass, for both links were increased by a value of 0.053 m for *c*<sup>1</sup> and a value of 0.022 m for *c*2. The length of the arm on axis 1 is less by 0.063 m. Due to these changes, of course, the values of the moments of inertia *I*<sup>1</sup> and *I*<sup>2</sup> have also changed. The Viscous friction coefficient and Coulomb friction magnitude were set experimentally before the validations. However, after optimization, the values of the Viscous friction coefficient and Coulomb friction magnitude are smaller than the initial value of 0.35 (especially the value of *fd*2).

The parameters that were the subject of estimation and used for the subsequent optimization of the model were estimated in the Parameter Estimator. During the estimation process, this tool (among other things) offers the possibility to monitor how the values of the parameters are estimated (changed) depending on the iterations until all the values of the estimated parameters are stabilized. The result of this phase from the given toolbox can be seen in Figures 31 and 32. During the estimation, only one graph could be viewed, where the waveforms of all 11 estimated parameters were. However, since the masses *m*<sup>1</sup> and *m*<sup>2</sup> are represented by much higher values than the remaining parameters, the mass curves were exported to the second graph for better clarity. Based on the waveforms of Figure 32, it can be seen that the smallest changes in the parameter estimation occurred in the parameters *c*<sup>2</sup> (Distance to the center of mass for link 2) and *Iz*<sup>2</sup> (Moment of inertia for link 2). On the contrary, significant changes took place in the values of Coulomb friction magnitude *fd*1, *fd*2. In Figure 31, it is possible to observe how the weight estimate has changed. From the course of estimating the mass *m*2, it can be seen how the weight changed in the first iterations, but the result of the estimation in the last iteration was, in the end, very similar to the original value.

**Figure 31.** Estimation parameters—part 1.

**Figure 32.** Estimation parameters—part 2.

Significant estimates and changes in some parameters can be evaluated based on the following hypotheses:


#### **9. Conclusions**

The research described in this article aimed to compare the model with the parameters that come from CAD and its improvement by estimating the parameters based on gravitational measurements. An improvements model could be subsequently applied to the control phase.

The study described in this article can be summarized in four stages:


Based on the results of parameter validation and estimation, it can be stated that the values of system parameters, which serve as input data for the model, significantly affect the validation results, which was demonstrated after estimating these parameters significantly affected the results. Within the course of the estimated parameters, a certain connection between the parameters was seen, which gives rise to the assumption that the values of the parameters not only influence the results of the validation but also interact with each other. The more accurate the system parameters are obtained in the step, the more accurate the validation results will be. By reducing the oscillations (especially for link 2 and at higher initial rotation angles) which occur in the system due to the individual structural bonds, it would be possible to achieve a more accurate stabilization of the arm movement to zero, and thus to improve the validation results.

In further research, the authors would like to extend this research to the following parts:


**Author Contributions:** Conceptualization, A.H. and M.T.; methodology, A.H.; model design, M.T. and T.C.; measuring data, A.H. and T.K.; validation in Matlab, T. ˇ C.; estimation in Matlab, T. ˇ C.; formal ˇ analysis, M.T.; investigation, M.T.; writing—original draft preparation, M.T.; writing—review and editing, M.T., A.H., T.K. and T.C.; visualization, M.T. and T. ˇ C.; project administration, T.K.; funding ˇ acquisition, A.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The project was supported by the Project VEGA no. 1/0393/18—Research of Methods for Modeling and Compensation of Hysteresis in Pneumatic Artificial Muscles and PAM-actuated Mechanisms to Improve the Control Performance Using Computational Intelligence and by the Project of the Structural Funds of the EU, ITMS code: 26220220103.

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

#### **References**


198

