TY - JOUR
AU - Palma-Huerta, Andrés Abraham
AB - Introduction Competition among industrial sectors generates new production and control system techniques that maximize processes’ efficiency and flexibility. Developing a control strategy for an automatic system requires physical prototyping, which can be time-consuming and expensive, especially when testing edge cases [1–3]. These tests can push the systems to unsafe conditions and may result in loss of material and equipment, among other expenses and risks. To overcome these difficulties, it is often desirable to rely on the use of simulation technology. Depending on the scope, the simulation environment, and the final purpose, simulation techniques or approaches can be classified into: Model-in-the-loop (MiL) [4]: The control logic is simulated using mathematical models of the system and its environment. Software-in-the-loop (SiL) [5]: The controller is implemented within simulation software that mimics the system’s behavior and its environment. Processor-in-the-loop (PiL) [6]: The control strategy is executed on the actual hardware or processor to govern the behavior of a system within a simulated environment. Hardware-in-the-loop (HiL) [7]: Hardware elements are used to simulate the system, the control strategy, and the operating environment. System-in-the-loop (XiL) [8]: Integrates all components of the system (hardware and software) into the simulation for thorough testing and validation of the entire system. Particularly, MiL offers various advantages: it is low-cost, fast to develop, can detect early problems with control strategies and system models, and is flexible and scalable in terms of the variety of system configurations and simulation conditions that can be imposed on the environment. It is crucial to study, generate, or utilize methods that generate models that closely mimic the real response of a system to develop the MiL approach properly. System identification encompasses various techniques, approaches, and methodologies to obtain models of dynamic systems from information estimated or acquired from their inputs and outputs. The identification process typically includes four main stages [9]: Data acquisition: Where the input and output variables that best represent the behavior of the dynamic system are measured or estimated. This measurement or estimation can be performed in the time or frequency domain for this type of system using sensors or transducers. Selection of a model structure: Here, a mathematical structure is derived that can relate the values of the system inputs to their corresponding outputs. Such structures can be simple algebraic relations, sets of differential equations, or more sophisticated computational models. Model structure tuning: Regardless of their type, model structures have several adjustable parameters on which the accuracy of relating the inputs to the system’s outputs depends on a correct settling. If the parameters are set properly, the model will be able to mimic the behavior to a large extent. On the other hand, a wrong adjustment of these parameters can trigger a behavior of the model very different from that of the real system. Evaluation of the adjusted model: This last stage verifies that the adjusted model obtained is adequate to satisfy the application’s needs. The four stages necessary to identify dynamic systems described above can be performed online or offline. In the offline approach, information on the system’s inputs and outputs is acquired at the beginning of the identification process, ensuring that the data are representative and sufficient and assuming that the system’s behavior will remain unchanged. Subsequently, the remaining three stages are carried out, and the validated model remains invariant and is used for some practical purposes [10], e.g., for the offline tuning of a control system [11]. On the other hand, online identification is a process that is carried out considering that the behavior of the dynamic system will change in the short term. Therefore, the four stages are repeated at certain time intervals during the operation of the dynamic system. One of the most common uses of this identification type is in the indirect adaptive tuning of controllers [12]. At present, different methods implement the above four stages to obtain adequate models of systems. Therefore, identification methods can be classified according to the type of information used to derive an accurate model of a system [13]. In this sense, there are time-domain methods, where information obtained from sensors in the system can be used directly during the identification process [14]. On the other hand, in frequency-based methods, the information from the sensors is not used directly but is converted to the frequency domain [15,16]. Parametric or non-parametric identification methods can be found depending on the model structure. Parametric methods search for a specific model structure’s parameter values that describe the system’s dynamic behavior [17,18]. On the other hand, non-parametric methods do not use a specific model structure but instead occupy generalized mathematical or computational structures with parameters that must be adjusted to replicate a particular system behavior [19,20]. System identification can be performed using state-space methods, which model a system’s behavior through its internal states with a minimal set of variables [21]. Alternatively, non-state-space methods focus only on input-output relationships to derive a model without accounting for internal dynamics [22]. Different models are frequently used among the identification approaches presented up to this point. Representative examples include linear [23] and nonlinear [24] autoregressive models; linear [25] and nonlinear [26] state-space models; transfer function models [27]; Gaussian regression processes [28]; and neural networks [29]. Now, concerning estimation algorithms, some common choices include least squares algorithms [30], classical optimization techniques [31], machine learning training methods [32], and metaheuristics [33]. The choice of a method, model, and estimator in system identification depends on several factors, including the nature of the studied dynamic system (e.g., whether it is complex in terms of its inputs and outputs, whether it can be described using established rules, whether it is subject to change, etc.), the characteristics of the available data (e.g., whether all input and output variables can be measured or estimated, whether the amount of information is sufficient and relevant, etc.), and the specific objectives of the analysis (i.e., the ultimate use of the identified model, e.g., for controller tuning or virtual simulation, etc.). System identification is especially relevant for deriving models of robotic systems. A robotic system can be defined as a collection of mechanical elements interconnected through joints that, when properly planned and controlled, allow very complex movements or tasks. These systems are growing in importance every day due to their ability to accurately perform automated tasks in different application areas ranging from healthcare [34] to space exploration [35]. In the context of robotic manipulators, many factors that can help select methods, models, and estimators for identification are evident. In the first instance, robotic manipulators possess nonlinear (or highly nonlinear) behaviors that can grow in complexity depending on their degrees of freedom. Since these systems are a subset of mechanical systems, it is possible to derive accurate model structures using the axioms of classical mechanics, i.e., Newton’s laws. These structures are naturally found in the state space and time domain and are parameterized by physical variables such as mass, inertia, friction, etc. The states used to describe the configuration of robotic systems at any instant include the positions and velocities of each of their joints, regardless of their type. Likewise, the inputs of this type of system are related to the amount of energy that is included in their motion-generating elements, i.e., in their actuators. In addition, the information on the inputs and outputs, related to the system states and energy inputs, can be acquired relatively simply using transducers or estimated using observers. Due to the inherent complexity of such systems, the identification problems are expected to be equally complex, so the use of metaheuristic techniques may be necessary. Table 1 provides a bibliographic review of recent studies that utilize metaheuristic optimization (Met. Opt.) for estimating model parameters in dynamic systems. This table reveals the frequent use of the parameter estimation in robotic applications, with 50% of the research, followed by the application in the actuators at 30%, and in other mechanical and electrical systems such as phase inverter circuits and mass-spring-damper systems at 20% of the reported works. Fig 1 aids in graphically depicting the algorithms utilized in the studies referenced in Table 1. It can be observed in Fig 1 that particle swarm optimization (PSO) is used in 15% of the studies, followed by the genetic algorithm (GA) at 10%, differential evolution (DE) algorithm and grey wolf optimization (GWO) at 8%. Algorithms such as the artificial bee colony (ABC) and gravitational search algorithm (GSA) accounted for 6% until 4% of the studies. Other algorithms represented 2% collectively in this review. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Percentage of times that metaheuristic algorithms are used in the parameter estimation of dynamic systems. https://doi.org/10.1371/journal.pone.0325168.g001 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Works related to the model’s parameter estimation by Met. Opt. https://doi.org/10.1371/journal.pone.0325168.t001 The limitations identified in the reviewed literature, summarized in Table 1, highlight that system model parameterizations are often tailored to specific applications or individual system components. This approach typically characterizes system dynamics within a restricted movement space or approximates the system using a linear model. As a result, the models tend to overfit the specific task or the assumed linear representation, limiting their generalization ability. It is important to mention that a good generalization is essential for MiL simulations. Another key limitation is that many reviewed studies focus on a single stage of the overall modeling process within a specific domain, aiming to enhance the parameterized model while overlooking other influential factors. For instance, some studies concentrate solely on algorithm development to refine model parameters but neglect other process stages that could further improve model accuracy. Similarly, others apply only one type of optimization algorithm, failing to explore alternative metaheuristic approaches that might yield better parameterization results. Additionally, the reviewed studies often lack flexibility in balancing trade-offs between optimization criteria. In many cases, this trade-off is determined arbitrarily, leading to a parameterization that prioritizes one objective at the expense of others. This limitation restricts the ability to explore more balanced solutions that could enhance overall model performance. Ultimately, these limitations suggest that existing studies do not provide models well-suited for MiL simulation environments. Moreover, a comprehensive, step-by-step methodology for optimizing robotic system parameters using metaheuristic techniques—along with its practical application in Model-in-the-Loop (MiL) testing, is noticeably absent from the literature. Motivated by these limitations and gaps, the main contributions of the paper are listed below: Proposal of the model-based parametric bio-inspired optimization (MbPBO) methodology: This work introduces an offline bio-inspired optimization-based model parameterization methodology for robotic systems, termed MbPBO. The MbPBO is designed for Model-in-the-Loop (MiL) testing and provides a structured approach to achieve accurate parameterization. This methodology significantly improves the alignment between simulation and real-world experimentation. Guidelines for MiL implementation: The MbPBO methodology offers a practical framework for researchers and practitioners looking to implement MiL test systems using bio-inspired optimization. It outlines essential steps and considerations to ensure effective model parameterization. Validation through case studies: The application of MbPBO in two case studies demonstrates the quality, consistency, and adaptability of the obtained models. These case studies highlight the model’s strong correlation with real experimental data and its potential for broader application in MiL test environments. Specifically, the methodology is implemented for a three-degree-of-freedom serial robot and a reaction force-sensing series elastic actuator, validating its flexibility and effectiveness across different electromechanical systems. The rest of this paper details the proposed MbPBO approach for robotic systems based on bio-inspired optimization. It also describes the implementation of the proposed approach on an RRR robot and an electromechanical system, the reaction force-sensing series elastic actuator, to show its scope. The discussions and conclusions of the proposal are drawn at the end of the paper. Model-based parametric bio-inspired optimization methodology The current section presents the proposed model-based parametric bio-inspired optimization (MbPBO) methodology. This methodology allows the dynamic behavior of systems to be identified based on bio-inspired optimization and then validated using the MiL approach. The proposed MbPBO methodology is observed in Fig 2. This figure shows the stages considered in the methodology, landed towards an application for any robotic system. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. General scheme of the MbPBO methodology. https://doi.org/10.1371/journal.pone.0325168.g002 The following subsections describe the stages that conform to the proposed MbPBO methodology. Stage I: Determination of the base model and its key parameters In this stage, the base model of the MbPBO methodology is defined. A white-box model is required, as it incorporates the mathematical relationships governing the system’s physical behavior. Black-box models are excluded because they obscure or eliminate these relationships, which are crucial for validation in the MiL stage. Additionally, black-box models tend to be less interpretable and are typically valid only within a limited workspace region. For mechanical systems, the direct dynamic model [74] is commonly used, as it is suitable for control and numerical simulation. In robotic systems, this model is typically represented by nonlinear differential equations and is preferably expressed in state-space form. Three aspects are considered at this stage: (i) the selection of parameters to be estimated, (ii) the initial values proposed for these parameters, and (iii) the reduction of the model by eliminating dynamics that do not substantially influence the system’s behavior. These aspects are discussed in detail below. (i) Selection of the parameters to be estimated. The model’s parameters define its behavior and must be carefully identified to achieve an accurate representation. These parameters primarily determine the system’s response to external stimuli. Since models are tailored to specific system properties, mathematically similar dynamic models, such as the closed-form dynamics of a robotic manipulator, can exist. However, the specific parameter values differentiate one system from another. In this stage, the n parameters to be identified are defined, forming the design variable vector , which will be used in subsequent stages of the MbPBO methodology. (ii) Initial values proposed for the parameters. As discussed later, the estimation techniques in the MbPBO methodology rely on metaheuristic algorithms, which typically begin their search with randomly initialized values for the design variables in p. However, in some cases, selecting appropriate initial values can enhance the final results. These values should be chosen based on estimates that are close to the real ones. For existing systems, it is recommended to use manufacturer-provided values as a starting point. When such data is unavailable or for entirely theoretical systems, computer-aided design (CAD) tools can be used to approximate key aspects, such as component dimensions and material properties. In the case of non-conservative models, parameters such as friction and hardness are often not provided by manufacturers or specified in CAD programs. Instead, these values can be estimated based on theoretical knowledge, ensuring the model remains solvable and well-defined. For example, the coefficient of dry friction varies depending on material type (e.g., metals, plastics), with accepted ranges documented in specialized literature [75], which can serve as a reference. (iii) Model reduction. In robotic systems, the complexity of the model and the number of required computations generally increase with the number of degrees of freedom (d.o.f.) [74]. Due to this complexity, analytical solutions are often impractical, necessitating the use of numerical methods and computational resources. Since the model is intended for simulation, it is crucial to account for the computational limitations of the simulator. Model complexity is determined by the computational resources required for accurate interpretation [76], often measured by the number of operations involved. To ensure efficient simulation execution, model simplification is necessary. This can be achieved by precomputing frequently recurring terms and applying mathematical properties to reduce redundant expressions. Simplification helps maintain model validity within the experimental framework while optimizing computational performance. Additionally, the design variable vector is updated accordingly. Stage II: Characterization of the system’s behavior In [77], various methods for obtaining system parameters are discussed. One approach involves direct physical measurement, which may require partial disassembly of components, making it impractical in some cases. Another option is to use manufacturer-provided values or computer-aided design (CAD) programs, though inaccuracies in manufacturing and assembly can complicate parameter estimation. A third and often more effective alternative for MiL applications is parametric identification, where parameters are inferred indirectly from the system’s dynamic behavior. At this stage, experimental data at different excitation frequencies are collected for use in model optimization. In MiL testing, system parameters are assumed to remain constant over time, so data collection is performed offline. The experiment involves extracting system response data under external excitation [76]. A key principle in system parameter characterization is the application of a persistent excitation signal [78], which induces a response that exhibits diverse dynamic behaviors. As highlighted in [79], such a signal provides sufficient a priori information about the system, improving parameter identification accuracy and robustness against measurement noise. To generate the excitation signal, [77] suggests using a control signal derived from an optimized trajectory, particularly for robotic manipulators. However, this approach results in parameters specific to the given trajectory, limiting generalization to other tasks. In this study, a sinusoidal or square wave signal with variable periods and bounded amplitude is used as the excitation signal, denoted as ref(t). Experimental system response data are collected based on the applied excitation signal ref(t). The number of recorded measurements depends on available resources, with accuracy closely tied to the quality of the measurement tools. In robotic systems, state variables () evolve over time and represent quantities such as position (), velocity (), acceleration (), among others. In electromechanical systems, excitation signals are typically controllable system variables, such as voltage. If a cascade control architecture is used, these signals correspond to the outer-loop reference signals of the system, indirectly modifying input variables and exciting system dynamics. Once collected, the experimental data will be utilized in the next stage of the methodology. Stage III: Establishment of a mathematical programming problem The parameters p are obtained by solving an optimization problem, which is formulated in this stage of the methodology. Once the experimental data of the system’s state variables (x) have been collected (Stage II) and the theoretical model has been defined (Stage I), the same excitation signal u(t) is applied to the model using a given vector of design variables p. The state error vector is then computed as the difference between the model’s estimated states and the system’s directly measured states over the experiment’s duration. Through a closed-loop, iterative identification process, the objective is to determine a vector that minimizes . In the MbPBO methodology, each objective function typically represents a state error variable. If multiple state errors share the same units, a weighted sum of these errors can also be used as an objective function. To ensure convexity, the methodology minimizes the squared state error, as defined in Eq. (1), where represents the number of state error elements in the i–th objective function (for a single error, ). (1) The optimization problem can be formulated using a multi-objective approach [80], expressed in equations (2) to (6), where represents the vector of objective functions, is the vector of design variables and and are the inequality and equality constraints that limits or provides precise conditions to the problem, respectively. Here, denotes the space of feasible solutions. Additionally, represents the set of states of the system model, where equation (3) describes its dynamics, and establish the minimum and maximum limits of the vector p, respectively. (2) subject to: (3)(4)(5)(6) In the MbPBO methodology, the multi-objective optimization problem is reformulated as a mono-objective optimization problem through the weighted sum approach [81]. In the weighted sum approach, the multi-objective optimization problem is converted into a single-objective problem by assigning a weight to each objective function included in the vector J (2). The objective function is then expressed as a weighted sum of all individual objectives. Adjusting the weights allows for exploring different trade-offs between objectives, allowing for a better understanding of the Pareto-optimal solutions. So, the single-objective function in (7) is the combination of nJ elements in the vector J, each expressed in (1) and weighted by . The i–th weight sets the priority (importance) in terms of the objective function for the optimization process, i.e., higher values indicate more preference in the minimization of the i–th term Ji in the single-objective function . A normalization process in terms of the objective function (Ji) is required to ensure dimensionless values within the range and to set the condition for the weight values as [81]. The normalization is performed by dividing the i–th term Ji by its maximum value . The maximum value is determined by evaluating Ji using the design variable vector obtained when all other terms are minimized individually. (7) Typically, the optimization problem in this stage can be formulated as shown in equations (8) to (12). The objective is to find the vector of design variables p that minimizes the weighted objective function in (7), subject to the system model’s dynamic behavior (9), inequality constraints (10), equality constraints (11), and limits on the design variables (12). (8) subject to: (9)(10)(11)(12) Stage IV: Solving the optimization problem using bio-inspired algorithms In this stage, the optimization problem formulated in Stage III is solved. To align the parameter vector with the experimental system’s behavior, the problem, originally defined in continuous time, is transformed into a discrete-time formulation using the transcription method [82]. This transformation enables the effective application of optimization techniques. The system’s governing differential equations, which serve as inherent constraints in the optimization problem, are discretized into finite states using numerical integration methods such as Euler’s method or Runge-Kutta methods. This approach ensures that dynamic constraints are satisfied simultaneously during the optimization process. Once the transcription is complete, the following steps are carried out: (i) Trade-off selection. Since the multi-objective problem has been transformed into a single-objective optimization problem using the weighted sum approach, it is crucial to establish a desired a priori preference among the terms in the objective function to improve the obtained model in the methodology. To achieve this, different weight configurations () are tested, each yielding a solution that reflects a specific trade-off in the weighted objective function. To determine the best trade-offs, the Pareto front is constructed by evaluating Pareto optimality [83] across all solutions obtained from different weight configurations. The identified Pareto-optimal solutions represent the most balanced trade-offs among the objective function terms. Finally, based on the resulting Pareto front, the designer can select the most suitable trade-off for the specific application. (ii) Algorithm selection. The proposed MbPBO methodology utilizes bio-inspired algorithms for system parameter estimation. These algorithms offer several advantages, including the ability to handle complex problems with reasonable computational costs, support for parallelization, ease of implementation, and adaptability to diverse contexts [84]. It is recommended to begin by selecting algorithms commonly used in the literature (see Table 1). The GA, PSO, and DE were chosen for analysis in this study. These algorithms are widely used due to their adaptability, population-based search strategies, robust performance, and ease of implementation [85–87]. Their effectiveness in both research and real-world applications has contributed to their growing popularity in optimization problems. (iii) Performance analysis of bio-inspired algorithms. Selecting the appropriate optimization technique based on empirical evidence is crucial for obtaining the best model parameters. Since no single algorithm outperforms all others across all problem types, an idea formalized by the “no free lunch theorem" [88], the MbPBO methodology recommends conducting a comparative study to evaluate the performance of representative algorithms. The comparative study involves running each selected algorithm independently with the same number of objective function evaluations. The best solution from each execution is stored, forming a statistical sample for comparison. To ensure robust statistical representation, thirty independent runs are recommended [89]. Additionally, algorithm parameters should be fine-tuned systematically through manual adjustment (“by hand"), informed by state-of-the-art practices “by analogy"), or optimized using automated parameter tuning methods [90]. Once the sample results are collected, statistical analysis is applied to determine which optimization technique achieves the lowest objective function value, indicating minimal error between the model’s estimated states and the plant’s measured states. Comparing algorithm performance is essential for selecting the best parameter set for the system. If significant differences are observed in the solutions obtained by different algorithms, it may be necessary to enhance algorithm operators or incorporate state-of-the-art modifications to improve search efficiency. Upon completing this stage, the optimized parameter vector is integrated into the system model, preparing it for validation in an open-loop control configuration and model-in-the-loop (MiL) simulation. Stage V: Validation of the obtained model and its use in the MiL simulation At this stage, the model undergoes validation by integrating the optimized parameter vector obtained in the previous section. This process consists of two key steps. The first step assesses whether the system model, with , exhibits behavior similar to the plant’s actual states when subjected to the excitation signal from Stage II. This is conducted in an open-loop control configuration. To quantify the model’s accuracy, a correlation measure (e.g., Pearson correlation coefficient) is recommended to evaluate the alignment between the model’s predicted states and the plant’s actual behavior. The second step evaluates the parameterized model within various MiL simulation environments. MiL validation involves a series of tests conducted in the early phases of model-based design to verify model behavior and performance. These tests are performed in a closed-loop control configuration within a simulation environment such as MATLAB/Simulink, assessing predictive accuracy for specific tasks [91]. Unlike software-in-the-loop (SiL), which focuses on control programming in a specific language, or hardware-in-the-loop (HiL), which evaluates component response times and communication delays, MiL validation is primarily concerned with verifying the model’s accuracy and performance before further implementation. Early detection of logical errors at this stage significantly reduces costs and improves model quality. Additionally, the model-based design allows for adaptability to evolving requirements, ensuring alignment with end-user needs [92]. The MbPBO methodology is designed for versatile applications, allowing for the simulation of tasks such as trajectory tracking and regulation, depending on the designer’s objectives. When the physical plant is available, the most effective validation involves deploying the model in a real-world task and comparing its performance against actual system execution [93]. By the end of this stage, the validated model should closely approximate the plant’s operational behavior, facilitating subsequent phases of the model-based design cycle (SiL to HiL). MiL validation is especially critical for applications where direct testing on the physical system may pose risks. Numerical simulations with a validated model provide valuable insights into expected system behavior under various conditions. Testbed platforms are available in the case studies presented below, allowing direct validation using real systems. The Pearson correlation coefficient is used as a metric to quantify the relationship between the model’s predictions and the actual system behavior across different test scenarios. Stage I: Determination of the base model and its key parameters In this stage, the base model of the MbPBO methodology is defined. A white-box model is required, as it incorporates the mathematical relationships governing the system’s physical behavior. Black-box models are excluded because they obscure or eliminate these relationships, which are crucial for validation in the MiL stage. Additionally, black-box models tend to be less interpretable and are typically valid only within a limited workspace region. For mechanical systems, the direct dynamic model [74] is commonly used, as it is suitable for control and numerical simulation. In robotic systems, this model is typically represented by nonlinear differential equations and is preferably expressed in state-space form. Three aspects are considered at this stage: (i) the selection of parameters to be estimated, (ii) the initial values proposed for these parameters, and (iii) the reduction of the model by eliminating dynamics that do not substantially influence the system’s behavior. These aspects are discussed in detail below. (i) Selection of the parameters to be estimated. The model’s parameters define its behavior and must be carefully identified to achieve an accurate representation. These parameters primarily determine the system’s response to external stimuli. Since models are tailored to specific system properties, mathematically similar dynamic models, such as the closed-form dynamics of a robotic manipulator, can exist. However, the specific parameter values differentiate one system from another. In this stage, the n parameters to be identified are defined, forming the design variable vector , which will be used in subsequent stages of the MbPBO methodology. (ii) Initial values proposed for the parameters. As discussed later, the estimation techniques in the MbPBO methodology rely on metaheuristic algorithms, which typically begin their search with randomly initialized values for the design variables in p. However, in some cases, selecting appropriate initial values can enhance the final results. These values should be chosen based on estimates that are close to the real ones. For existing systems, it is recommended to use manufacturer-provided values as a starting point. When such data is unavailable or for entirely theoretical systems, computer-aided design (CAD) tools can be used to approximate key aspects, such as component dimensions and material properties. In the case of non-conservative models, parameters such as friction and hardness are often not provided by manufacturers or specified in CAD programs. Instead, these values can be estimated based on theoretical knowledge, ensuring the model remains solvable and well-defined. For example, the coefficient of dry friction varies depending on material type (e.g., metals, plastics), with accepted ranges documented in specialized literature [75], which can serve as a reference. (iii) Model reduction. In robotic systems, the complexity of the model and the number of required computations generally increase with the number of degrees of freedom (d.o.f.) [74]. Due to this complexity, analytical solutions are often impractical, necessitating the use of numerical methods and computational resources. Since the model is intended for simulation, it is crucial to account for the computational limitations of the simulator. Model complexity is determined by the computational resources required for accurate interpretation [76], often measured by the number of operations involved. To ensure efficient simulation execution, model simplification is necessary. This can be achieved by precomputing frequently recurring terms and applying mathematical properties to reduce redundant expressions. Simplification helps maintain model validity within the experimental framework while optimizing computational performance. Additionally, the design variable vector is updated accordingly. (i) Selection of the parameters to be estimated. The model’s parameters define its behavior and must be carefully identified to achieve an accurate representation. These parameters primarily determine the system’s response to external stimuli. Since models are tailored to specific system properties, mathematically similar dynamic models, such as the closed-form dynamics of a robotic manipulator, can exist. However, the specific parameter values differentiate one system from another. In this stage, the n parameters to be identified are defined, forming the design variable vector , which will be used in subsequent stages of the MbPBO methodology. (ii) Initial values proposed for the parameters. As discussed later, the estimation techniques in the MbPBO methodology rely on metaheuristic algorithms, which typically begin their search with randomly initialized values for the design variables in p. However, in some cases, selecting appropriate initial values can enhance the final results. These values should be chosen based on estimates that are close to the real ones. For existing systems, it is recommended to use manufacturer-provided values as a starting point. When such data is unavailable or for entirely theoretical systems, computer-aided design (CAD) tools can be used to approximate key aspects, such as component dimensions and material properties. In the case of non-conservative models, parameters such as friction and hardness are often not provided by manufacturers or specified in CAD programs. Instead, these values can be estimated based on theoretical knowledge, ensuring the model remains solvable and well-defined. For example, the coefficient of dry friction varies depending on material type (e.g., metals, plastics), with accepted ranges documented in specialized literature [75], which can serve as a reference. (iii) Model reduction. In robotic systems, the complexity of the model and the number of required computations generally increase with the number of degrees of freedom (d.o.f.) [74]. Due to this complexity, analytical solutions are often impractical, necessitating the use of numerical methods and computational resources. Since the model is intended for simulation, it is crucial to account for the computational limitations of the simulator. Model complexity is determined by the computational resources required for accurate interpretation [76], often measured by the number of operations involved. To ensure efficient simulation execution, model simplification is necessary. This can be achieved by precomputing frequently recurring terms and applying mathematical properties to reduce redundant expressions. Simplification helps maintain model validity within the experimental framework while optimizing computational performance. Additionally, the design variable vector is updated accordingly. Stage II: Characterization of the system’s behavior In [77], various methods for obtaining system parameters are discussed. One approach involves direct physical measurement, which may require partial disassembly of components, making it impractical in some cases. Another option is to use manufacturer-provided values or computer-aided design (CAD) programs, though inaccuracies in manufacturing and assembly can complicate parameter estimation. A third and often more effective alternative for MiL applications is parametric identification, where parameters are inferred indirectly from the system’s dynamic behavior. At this stage, experimental data at different excitation frequencies are collected for use in model optimization. In MiL testing, system parameters are assumed to remain constant over time, so data collection is performed offline. The experiment involves extracting system response data under external excitation [76]. A key principle in system parameter characterization is the application of a persistent excitation signal [78], which induces a response that exhibits diverse dynamic behaviors. As highlighted in [79], such a signal provides sufficient a priori information about the system, improving parameter identification accuracy and robustness against measurement noise. To generate the excitation signal, [77] suggests using a control signal derived from an optimized trajectory, particularly for robotic manipulators. However, this approach results in parameters specific to the given trajectory, limiting generalization to other tasks. In this study, a sinusoidal or square wave signal with variable periods and bounded amplitude is used as the excitation signal, denoted as ref(t). Experimental system response data are collected based on the applied excitation signal ref(t). The number of recorded measurements depends on available resources, with accuracy closely tied to the quality of the measurement tools. In robotic systems, state variables () evolve over time and represent quantities such as position (), velocity (), acceleration (), among others. In electromechanical systems, excitation signals are typically controllable system variables, such as voltage. If a cascade control architecture is used, these signals correspond to the outer-loop reference signals of the system, indirectly modifying input variables and exciting system dynamics. Once collected, the experimental data will be utilized in the next stage of the methodology. Stage III: Establishment of a mathematical programming problem The parameters p are obtained by solving an optimization problem, which is formulated in this stage of the methodology. Once the experimental data of the system’s state variables (x) have been collected (Stage II) and the theoretical model has been defined (Stage I), the same excitation signal u(t) is applied to the model using a given vector of design variables p. The state error vector is then computed as the difference between the model’s estimated states and the system’s directly measured states over the experiment’s duration. Through a closed-loop, iterative identification process, the objective is to determine a vector that minimizes . In the MbPBO methodology, each objective function typically represents a state error variable. If multiple state errors share the same units, a weighted sum of these errors can also be used as an objective function. To ensure convexity, the methodology minimizes the squared state error, as defined in Eq. (1), where represents the number of state error elements in the i–th objective function (for a single error, ). (1) The optimization problem can be formulated using a multi-objective approach [80], expressed in equations (2) to (6), where represents the vector of objective functions, is the vector of design variables and and are the inequality and equality constraints that limits or provides precise conditions to the problem, respectively. Here, denotes the space of feasible solutions. Additionally, represents the set of states of the system model, where equation (3) describes its dynamics, and establish the minimum and maximum limits of the vector p, respectively. (2) subject to: (3)(4)(5)(6) In the MbPBO methodology, the multi-objective optimization problem is reformulated as a mono-objective optimization problem through the weighted sum approach [81]. In the weighted sum approach, the multi-objective optimization problem is converted into a single-objective problem by assigning a weight to each objective function included in the vector J (2). The objective function is then expressed as a weighted sum of all individual objectives. Adjusting the weights allows for exploring different trade-offs between objectives, allowing for a better understanding of the Pareto-optimal solutions. So, the single-objective function in (7) is the combination of nJ elements in the vector J, each expressed in (1) and weighted by . The i–th weight sets the priority (importance) in terms of the objective function for the optimization process, i.e., higher values indicate more preference in the minimization of the i–th term Ji in the single-objective function . A normalization process in terms of the objective function (Ji) is required to ensure dimensionless values within the range and to set the condition for the weight values as [81]. The normalization is performed by dividing the i–th term Ji by its maximum value . The maximum value is determined by evaluating Ji using the design variable vector obtained when all other terms are minimized individually. (7) Typically, the optimization problem in this stage can be formulated as shown in equations (8) to (12). The objective is to find the vector of design variables p that minimizes the weighted objective function in (7), subject to the system model’s dynamic behavior (9), inequality constraints (10), equality constraints (11), and limits on the design variables (12). (8) subject to: (9)(10)(11)(12) Stage IV: Solving the optimization problem using bio-inspired algorithms In this stage, the optimization problem formulated in Stage III is solved. To align the parameter vector with the experimental system’s behavior, the problem, originally defined in continuous time, is transformed into a discrete-time formulation using the transcription method [82]. This transformation enables the effective application of optimization techniques. The system’s governing differential equations, which serve as inherent constraints in the optimization problem, are discretized into finite states using numerical integration methods such as Euler’s method or Runge-Kutta methods. This approach ensures that dynamic constraints are satisfied simultaneously during the optimization process. Once the transcription is complete, the following steps are carried out: (i) Trade-off selection. Since the multi-objective problem has been transformed into a single-objective optimization problem using the weighted sum approach, it is crucial to establish a desired a priori preference among the terms in the objective function to improve the obtained model in the methodology. To achieve this, different weight configurations () are tested, each yielding a solution that reflects a specific trade-off in the weighted objective function. To determine the best trade-offs, the Pareto front is constructed by evaluating Pareto optimality [83] across all solutions obtained from different weight configurations. The identified Pareto-optimal solutions represent the most balanced trade-offs among the objective function terms. Finally, based on the resulting Pareto front, the designer can select the most suitable trade-off for the specific application. (ii) Algorithm selection. The proposed MbPBO methodology utilizes bio-inspired algorithms for system parameter estimation. These algorithms offer several advantages, including the ability to handle complex problems with reasonable computational costs, support for parallelization, ease of implementation, and adaptability to diverse contexts [84]. It is recommended to begin by selecting algorithms commonly used in the literature (see Table 1). The GA, PSO, and DE were chosen for analysis in this study. These algorithms are widely used due to their adaptability, population-based search strategies, robust performance, and ease of implementation [85–87]. Their effectiveness in both research and real-world applications has contributed to their growing popularity in optimization problems. (iii) Performance analysis of bio-inspired algorithms. Selecting the appropriate optimization technique based on empirical evidence is crucial for obtaining the best model parameters. Since no single algorithm outperforms all others across all problem types, an idea formalized by the “no free lunch theorem" [88], the MbPBO methodology recommends conducting a comparative study to evaluate the performance of representative algorithms. The comparative study involves running each selected algorithm independently with the same number of objective function evaluations. The best solution from each execution is stored, forming a statistical sample for comparison. To ensure robust statistical representation, thirty independent runs are recommended [89]. Additionally, algorithm parameters should be fine-tuned systematically through manual adjustment (“by hand"), informed by state-of-the-art practices “by analogy"), or optimized using automated parameter tuning methods [90]. Once the sample results are collected, statistical analysis is applied to determine which optimization technique achieves the lowest objective function value, indicating minimal error between the model’s estimated states and the plant’s measured states. Comparing algorithm performance is essential for selecting the best parameter set for the system. If significant differences are observed in the solutions obtained by different algorithms, it may be necessary to enhance algorithm operators or incorporate state-of-the-art modifications to improve search efficiency. Upon completing this stage, the optimized parameter vector is integrated into the system model, preparing it for validation in an open-loop control configuration and model-in-the-loop (MiL) simulation. (i) Trade-off selection. Since the multi-objective problem has been transformed into a single-objective optimization problem using the weighted sum approach, it is crucial to establish a desired a priori preference among the terms in the objective function to improve the obtained model in the methodology. To achieve this, different weight configurations () are tested, each yielding a solution that reflects a specific trade-off in the weighted objective function. To determine the best trade-offs, the Pareto front is constructed by evaluating Pareto optimality [83] across all solutions obtained from different weight configurations. The identified Pareto-optimal solutions represent the most balanced trade-offs among the objective function terms. Finally, based on the resulting Pareto front, the designer can select the most suitable trade-off for the specific application. (ii) Algorithm selection. The proposed MbPBO methodology utilizes bio-inspired algorithms for system parameter estimation. These algorithms offer several advantages, including the ability to handle complex problems with reasonable computational costs, support for parallelization, ease of implementation, and adaptability to diverse contexts [84]. It is recommended to begin by selecting algorithms commonly used in the literature (see Table 1). The GA, PSO, and DE were chosen for analysis in this study. These algorithms are widely used due to their adaptability, population-based search strategies, robust performance, and ease of implementation [85–87]. Their effectiveness in both research and real-world applications has contributed to their growing popularity in optimization problems. (iii) Performance analysis of bio-inspired algorithms. Selecting the appropriate optimization technique based on empirical evidence is crucial for obtaining the best model parameters. Since no single algorithm outperforms all others across all problem types, an idea formalized by the “no free lunch theorem" [88], the MbPBO methodology recommends conducting a comparative study to evaluate the performance of representative algorithms. The comparative study involves running each selected algorithm independently with the same number of objective function evaluations. The best solution from each execution is stored, forming a statistical sample for comparison. To ensure robust statistical representation, thirty independent runs are recommended [89]. Additionally, algorithm parameters should be fine-tuned systematically through manual adjustment (“by hand"), informed by state-of-the-art practices “by analogy"), or optimized using automated parameter tuning methods [90]. Once the sample results are collected, statistical analysis is applied to determine which optimization technique achieves the lowest objective function value, indicating minimal error between the model’s estimated states and the plant’s measured states. Comparing algorithm performance is essential for selecting the best parameter set for the system. If significant differences are observed in the solutions obtained by different algorithms, it may be necessary to enhance algorithm operators or incorporate state-of-the-art modifications to improve search efficiency. Upon completing this stage, the optimized parameter vector is integrated into the system model, preparing it for validation in an open-loop control configuration and model-in-the-loop (MiL) simulation. Stage V: Validation of the obtained model and its use in the MiL simulation At this stage, the model undergoes validation by integrating the optimized parameter vector obtained in the previous section. This process consists of two key steps. The first step assesses whether the system model, with , exhibits behavior similar to the plant’s actual states when subjected to the excitation signal from Stage II. This is conducted in an open-loop control configuration. To quantify the model’s accuracy, a correlation measure (e.g., Pearson correlation coefficient) is recommended to evaluate the alignment between the model’s predicted states and the plant’s actual behavior. The second step evaluates the parameterized model within various MiL simulation environments. MiL validation involves a series of tests conducted in the early phases of model-based design to verify model behavior and performance. These tests are performed in a closed-loop control configuration within a simulation environment such as MATLAB/Simulink, assessing predictive accuracy for specific tasks [91]. Unlike software-in-the-loop (SiL), which focuses on control programming in a specific language, or hardware-in-the-loop (HiL), which evaluates component response times and communication delays, MiL validation is primarily concerned with verifying the model’s accuracy and performance before further implementation. Early detection of logical errors at this stage significantly reduces costs and improves model quality. Additionally, the model-based design allows for adaptability to evolving requirements, ensuring alignment with end-user needs [92]. The MbPBO methodology is designed for versatile applications, allowing for the simulation of tasks such as trajectory tracking and regulation, depending on the designer’s objectives. When the physical plant is available, the most effective validation involves deploying the model in a real-world task and comparing its performance against actual system execution [93]. By the end of this stage, the validated model should closely approximate the plant’s operational behavior, facilitating subsequent phases of the model-based design cycle (SiL to HiL). MiL validation is especially critical for applications where direct testing on the physical system may pose risks. Numerical simulations with a validated model provide valuable insights into expected system behavior under various conditions. Testbed platforms are available in the case studies presented below, allowing direct validation using real systems. The Pearson correlation coefficient is used as a metric to quantify the relationship between the model’s predictions and the actual system behavior across different test scenarios. Case study of the MbPBO methodology: Modeling of an RRR robot A representative case study of robotic system is proposed to verify the effectiveness of the MbPBO methodology. In this case, a planar RRR manipulator robot is selected. This robot consists of three rigid links connected by rotational joints that allow relative movement between them. It is equipped with sensors to obtain the relative position of the links and actuators (motors) to enable their movement. These components allow the robot to position its end-effector, on a two-dimensional plane defined by the coordinates . The manipulator is depicted in the schematic diagram in Fig 3, where the position vector q=[q1,q2,q3] represents the rotational movements in each link, defining its generalized coordinates. Similarly, represents the velocity vector. The dynamic parameters of the i-th link, which are the subject of experimental identification, include the mass mi, the distance from the axis of movement to the center of mass , the moment of inertia , and the length of the link li, with . Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Schematic diagram of the serial robot with three degrees of freedom (d.o.f.). https://doi.org/10.1371/journal.pone.0325168.g003 This system was chosen because it has an open architecture that enables the implementation of different control strategies and tasks that are useful in the validation stage. The experimental platform of the RRR robot, shown in the diagram in Fig 4, consists of a custom-designed and manufactured non-commercial robot with three degrees of freedom. It includes three DC motors (RE-35 from Maxon Motor), their respective gearboxes, and rotary encoders (MR-256 from Maxon Motor). The DC motors are driving by three servo amplifiers LMD18245 configured in the current mode configuration (torque mode configuration). The current mode configuration incorporates a proportional-integral (PI) controller to regulate the armature current of the robot’s motors. So, a cascade control strategy is set up to govern the robot behavior, where the inner loop contains the PI current controller and the outer loop controller provides the current reference to the inner loop one. The current reference, in turn, modifies the system’s input variable (voltage), which excites the system dynamics and the users provide this reference signal. The power supply for the robot is also included, providing a 30V, 6A DC voltage source. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Component diagram of the experimental platform of the RRR robot. https://doi.org/10.1371/journal.pone.0325168.g004 The outer loop controller programming has been carried out on a commercial PC (Intel Core 2 Quad Q8400 at 2.66 GHz with 4 GB of RAM), and data acquisition from the robot is performed through a Sensoray 626 card, which is adapted to one of the PCI ports of the computer system. So, in the open-loop control operation, the outer loop is made by injecting the reference signal (current reference) without feedback into the cascade system. In a closed-loop control operation, the reference signal and feedback are set to the controller. The experimental platform has native integration with MATLAB/Simulink, which is desirable for conducting tests in a MiL environment. Different control laws can be implemented in numerical simulations for the outer loop to validate the obtained model. This setup allows measurable results in the physical system to directly verify the proposed methodology’s effectiveness. The next section describes in detail the application of the MbPBO methodology in the RRR robot case study mentioned. Stage I: Determination of the base model for the RRR robot and its key parameters The considered dynamic system is an autonomous system; as such, its physical parameters are invariant in time. This characteristic makes it suitable for creating a parametric model where the parameters can be identified offline [93]. Additionally, the model is also proposed as a non-conservative system. The frictions between links present a challenge objective for the optimal parameterization proposed by the MbPBO methodology because specifying these parameters is difficult for any manufacturer. A dynamic model of the system is chosen for the primary mathematical modeling effort, where the relationship between the movement and the forces that produce it is explicitly described. For its formulation, it was decided to use the Euler-Lagrange formalism, whose principle is the balance of the kinetic and potential energy of the components of the system. The basis of this model is the Lagrangian equation (13) where and P(q) represent, respectively, the kinetic and potential energy of the links. (13) Likewise, when considering a non-conservative system where energy is dissipated, the Rayleigh dissipation function is included in this formulation. This function models these friction forces as viscous dampers in the term D shown in (14). (14) The term bi in this last equation represents the damper coefficients modeled for the i–th link. Including this dissipation function in the model, the corresponding Lagrange equations for the modeling of non-conservative systems are expressed in equation (15). (15) The model obtained from [94] developed with the Euler-Lagrange energy method was abbreviated in terms to assume the limitations of the simulator and promote shorter calculation times. The latter process considers the third aspect of this first stage related to the model reduction. The effect of viscous friction for each link is represented as , and likewise, the vector τ=[τ1,τ2,τ4] represents the input vector of torques. The closed form of the robot manipulator’s dynamic model is described in (16), where represents the inertia matrix, is the matrix of Coriolis and centrifugal forces, G∈ℝ3 is the gravity vector, and corresponds to the vector of torques produced by friction. (16) where: (17)(18)(19)(20) with: (21)(22)(23) The RRR robot dynamics (16) in the state vector with the control signal is expressed in a simple form as . The following is given concerning the steps (i) and (ii) of the Stage I. (i) Selection of the parameters to be estimated. It is then defined for this system (based on Fig 3, and the previously presented model), the vector p (24) as the set of parameters to be obtained in the identification of the RRR robot. The vector in (24) includes 12 elements to be estimated, where denotes the vector of distances between the axis of movement of each link to its center of mass, is the vector of inertia, is the mass vector, and refers to the viscous friction vector. The parameters in p are essential for identifying the RRR robot’s dynamic model using the MbPBO methodology. (24) (ii) Initial values proposed for the parameters. In the process of parameter identification, it is beneficial to begin with prior knowledge of values close to the parameters to be estimated. This approach allows for a constrained search, enhancing the efficiency of model refinement. The initial values for the parameter vector (24), grouped in , are sourced from [94]. In that study, the SolidWorks mechanical modeling program was utilized to estimate the parameters of the robot, such as mass, inertia, and mass center of links, as detailed in Table 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Proposed initial dynamic and kinematic parameters for the RRR robot. https://doi.org/10.1371/journal.pone.0325168.t002 In this case, neither the manufacturer nor a CAD program can provide the friction parameter . This parameter pertains to the coefficient of friction for dry surfaces, also known as dry friction, which varies depending on the materials in contact. The coefficient of friction between two dry metallic surfaces typically ranges from 0.15 - 0.8 [75]. Therefore, the known range of friction coefficients for metallic surfaces and factors that could result in the highest friction coefficient (like oxidation and high contact pressure, to which the robot joints could be exposed) were considered when proposing the friction coefficient values for the parameter estimation process. In this particular case, one of the worst cases is selected as the initial value for the friction coefficients. So, those values are initially set as 0.9. Stage II: Characterization of the RRR robot’s behavior. To characterize the system behavior, identification experiment on the RRR robot is developed. The experiment used open-loop control signals (open-loop external signal into the outer loop of the cascade system) along with their respective periods, following persistent excitation signals recommended in [95] for this purpose. Each reference signal refi to the inner loop controller for each degree of freedom (d.o.f.) has a square wave signal. These signals varied between the maximum applied current that the actuators could handle. Specific periods T were interleaved within each reference signal, ensuring an experimentation time with a sampling time of . The experiment induced various behaviors in the robot’s states by applying these signals ref in an open-loop control configuration. In this case study, we focused on monitoring the positions and velocities of the robot’s links . A representative diagram of this experimental setup can be seen in Fig 5. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Schematic diagram of the open-loop control experiment conducted on the experimental platform to excite the system states. https://doi.org/10.1371/journal.pone.0325168.g005 The individual periods applied to the three degrees of freedom are: , and . The graphical representation of the reference signals ref1, ref2, and ref3 in the outer loop controller to provide the reference in the inner loop controller of the robot is visualized in Fig 6. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Reference signals used in the open-loop controller on the experimental platform to excite the corresponding system states. https://doi.org/10.1371/journal.pone.0325168.g006 At the end of this stage, a data vector of position states defined as is obtained, which is directly retrieved from the robot’s encoder over the period of (resulting in a state measurement vector of 2600 samples). These measurements are graphically represented with blue lines in the plots of Fig 7a, 7c, 7e. The velocity state, defined as q˙=[q1˙,q2˙,q3˙]T∈R3, is estimated using numerical differentiation in MATLAB. The plots of Fig 7b, 7d, 7f illustrate the estimated velocity values. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Graphs showing the position state vectors q and the estimated velocity state vectors of the real system when the system states are excited by the open-loop control signals. https://doi.org/10.1371/journal.pone.0325168.g007 Stage III: Establishment of a mathematical programming problem for the RRR robot characterization. With the information from the previous stage, the current reference signal ref established in the previous stage is ideally transformed into torque input for only simulation purposes by the equation [96], where km = 1.9839 is the torque constant. So, the control signal is applied in simulation to the degrees of freedom in the robot model defined in Stage I. Its numerical solution provides us with a vector of estimated positions and velocities, defined as and , respectively. The error vector (25) is then defined between the actual position and the estimated position , as well as, the error vector of the velocity states between the experimental data and those associated with the model (26). Thus, the position and velocity errors of the i-th link are represented by the vectors ei and , respectively. (25)(26) So, in the optimization problem, the goal is to minimize the error between the signals obtained from the experiment in Fig 7 and the output of the simulated system in the time domain (error dynamics). Therefore, it is imperative to minimize errors in the states. For this purpose, the criterion of minimizing the integral of the square error (ISE) is chosen, as it is most suitable for this purpose [97]. Once the error vectors have been defined, the objective function terms are described in the space of positions and velocities in (27) and (28), respectively. (27)(28) Following the indication given in Stage III of the MbPBO methodology description, the multi-objective nature of the problem is transformed into a single-objective problem by implementing a weighted sum approach. Here, dimensionless values and represent the weights of the position and velocity errors, respectively. The objective function is minimized based on a normalized weighted expression . For this particular case study, the objective function can be defined as in (29), where only two terms are included to provide the same importance for the position errors and other priority for the velocity errors. (29) In (29), the scalar values and are quantities that represent the maximum possible values of J1 and J2 that will produce a dimensionless and normalized value of (ranging between 0 and 1). Those maximum values are obtained by only minimizing one term in the optimization process, i.e., by only minimizing J1 (27) or J2 (28), and the obtained optimal design variable vector or is evaluated in the other terms. So, in the particular case, the maximum value is calculated by evaluating the term J2 with the obtained vector (when the J1 term is minimized), i.e., . In a similar fashion, the maximum value is given by minimizing J2, obtaining the optimal design variable vector and calculated the maximum value as . The specific values of and are provided in the next Stage IV. The dynamic constraints (35) related to the system behavior are included in the optimization problem to bind the problem solution to the dynamic behavior of the RRR robot (16). A tolerance range is established around the initial (base) value of each variable pi to limit the search space in the optimization process. The increment or reduction factor is denoted as , where represents the desired percentage change over the base value shown in Table 2. The chosen increment/reduction factor values impose practical boundaries in the design variables to ensure that the search space in the optimization process presents valid solutions. The optimization process with these bounds is also more efficient because it does not waste time looking for solutions where the system can not operate. An increase and decrease factor of at most thirty percent is recommended for this purpose. Equations (30) to (33) describe such a limits () and define how the bound constraints are given. The selected upper and lower bound ranges ( ) are obtained from their initial (base) values and the increment/reduction factor given in Table 3. So, the bounds of the design variable vector (30)-(33) are the other constraints taken into account in the optimization problem. In (29), the scalar values and are quantities that represent the maximum possible values of J1 and J2 that will produce a dimensionless and normalized value of (ranging between 0 and 1). Those maximum values are obtained by only minimizing one term in the optimization process, i.e., by only minimizing J1 (27) or J2 (28), and the obtained optimal design variable vector or is evaluated in the other terms. So, in the particular case, the maximum value is calculated by evaluating the obtained vector when the J1 term is minimized, i.e., . In a similar fashion, the maximum value is given by minimizing J2, obtaining the optimal design variable vector and with this information can compute the maximum value as . The specific values of and are provided in the next Stage IV. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. Percentage rate in the design variable changes and the design variable bounds used in the optimization process. The specific quantities are indicated with . https://doi.org/10.1371/journal.pone.0325168.t003 The dynamic constraints (35) related to the system behavior are included in the optimization problem to bind the problem solution to the dynamic behavior of the RRR robot (16). A tolerance range is established around the initial (base) value () of each variable pi to limit the search space in the optimization process. The increment or reduction factor is denoted as , where represents the desired percentage change over the base value shown in Table 2. The chosen increment/reduction factor values impose practical boundaries in the design variables to ensure that the search space in the optimization process presents valid solutions. The optimization process with these bounds is also more efficient because it does not waste time looking for solutions where the system can not operate. An increase and decrease factor of at most thirty percent is recommended for this purpose. Equations (30) to (33) describe such a limits () and define how the bound constraints are given. The selected upper and lower bound ranges () are obtained from their initial (base) values given in Table 2 and the increment/reduction factor given in Table 3. So, the bounds of the design variable vector (30)-(33) are the other constraints taken into account in the optimization problem. (30)(31)(32)(33) Relating the above with the formal approach, the problem consists of finding the vector of design variables p that minimizes the weighted objective function (34), subject to the dynamic behavior of the system model (35), and the constraints on the design variables (36). In this case, the equality constraints are implicitly satisfied by the solution of the model dynamics. (34) subject to: (35)(36) The technique applied to obtain the minimum among these errors will yield the vector of optimal design variables p*, where and n = 12 (24). The optimization techniques used to solve this problem will be addressed in the following Stage IV. Stage IV: Solving the RRR robot’s characterization problem using bio-inspired algorithms It is important to point out that both terms of the objective function (29) present trade-offs for fulfilling both criteria. This means that minimizing both errors is not possible, i.e., the relationship between them turns out to be inverse; they behave like two opposing objectives where only minimizing the first term J1 implies maximizing the second one J2 and vice versa. In order to limit the weight values and to , it is important to normalize the terms (to set the values of such terms in the interval between ) of the objective functions. So, the terms max(J1) and max(J2) of (29) are obtained by individually minimizing the terms J2 and J1, respectively, and using the same constraints in (35)-(36), as described in the previous stage. In this part of the methodology, it is recommended to set several independent runs (for instance, thirty times [89]) for solving both optimization problems due to the stochastic behavior of the algorithm. DE/rand/1/bin is used in this case, but other frequently used bio-inspired optimizers can also be used. The best execution, the one that has the minimum value of the term, is selected from the thirty runs. The result that further minimizes the term J1 provides the objective function vector . Meanwhile, the result that further minimizes the term J2 provides the objective function vector . The values in boldface are the possible extremes sought, i.e., and . For practical simplicity, the values that are considered to carry out the normalization in the equation (29) are and . Once the maximum values of the terms in the optimization problem are obtained, the three steps considered in this stage are described below. (i) Trade-off selection. As mentioned in the description of the MbPBO methodology, the multi-objective problem is transformed into a single-objective problem by a weighted sum approach. A set of non-dominated solutions must be generated by establishing different weights in the optimization problem (34)-(36), forming the approximate Pareto front in the performance function space. The Pareto front helps the designer make decisions on the most suitable trade-off for the specific application. Once the terms of the objective function are normalized, new experiments are carried out to find the best trade-offs in the terms J1 and J2 (quadratic errors of position and velocity). Then, fourteen combinations of weights and in the interval are tested. The first two combinations find the solution of the Pareto front without a trade-off, i.e., by setting and . The next nine combinations consider uniform intervals from to subject to  +  . The three final weights are assigned according to the best promising weight interval found in the previous experiment. The promising weight interval is visualized once the previous experiments with different weights are carried out. Thirty independent executions of the bio-inspired algorithm are performed by establishing different weights in the optimization problem (34)-(36) and using the same DE algorithm with the same parameters given above. This experiment aims to identify the combination of weights that results in a helpful trade-off solution by including the best solution (the solution that further minimizes the performance function from the algorithm executions) found in each combination. Table 4 shows the chosen weight values, the weighted normalized objective function , and its weighted non-normalized terms J1 and J2. The values in boldface indicate the minimum ones in each column. As the uniform distribution of weights does not guarantee uniform distribution of the found Pareto solutions [98], and even the bio-inspired optimizer can stagnate in local solutions [83], the obtained solutions are sorted according to non-dominance verification. The Pareto solutions [83] of results from Table 4 are displayed in Table 5 and the graphical representation of the Pareto front with the normalized terms and is shown in Fig 8. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Normalized Pareto front obtained by using different weights in the optimization process. https://doi.org/10.1371/journal.pone.0325168.g008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Values of the objective function and its non-normalized terms , obtained by using different weights in the optimization process. https://doi.org/10.1371/journal.pone.0325168.t004 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 5. Set of Pareto solutions obtained by the results of Table 4. https://doi.org/10.1371/journal.pone.0325168.t005 Once the Pareto front is obtained, the decision maker’s preference provides one useful trade-off for the application. In this case, the most preferred solution is the Pareto solution with the minimum value in the term J1 and with a suitable value in the term J2. Based on Table 5 and Fig 8, the combination of weights , stands out for its ability to further minimize the trade-off in the performance function and the position quadratic errors, also providing an acceptable velocity quadratic error. The selected weights are employed in the next step to analyze different bio-inspired algorithms to optimize the RRR robot’s characterization problem. (ii) Algorithm selection. The choice of the DE/rand/1/bin algorithm for the experiments performed previously is because it is one of the most widely used alternatives in the reviewed literature [99] (Fig 1). However, the MbPBO methodology description in this work suggests using two more techniques to solve the optimization problem with the selected obtained weights (, ). Therefore, the next frequently used algorithms in the reviewed literature are selected due to their adaptability, population-based strategy, robust performance, and simplicity [85–87]. In particular, the Particle Swarm Optimization (PSO) [39] with the best global topology and the Genetic Algorithm (GA) in its generational variant [100] are also selected to optimize the RRR robot’s characterization problem described before (34)-(36). As mentioned in the methodology description, it is only possible to determine which algorithm presents the best results for a specific problem once it is experimentally proven [88]. Therefore, given the previous algorithms, the present work proposes a comparison based on the results obtained to identify which bio-inspired technique generally performs best when solving the characterization problem. The algorithm parameters must be tuned either “by hand” through systematic manual settings, “by analogy” through the information suggested in the state-of-the-art for similar problems, or using self-tuning mechanisms. The process “by hand” (where algorithm parameters are varied sequentially and ordered to select the best configuration that solves the problem) is also suggested in [101] for algorithms like DE. Therefore, the latter is the process chosen in this work to tune such parameters. The set of parameters presented below was obtained from this trial and error procedure, where the final tuning values for the DE/rand/1/bin algorithm are as follows: Maximum number of generations Gmax = 500, number of individuals NP = 100, scaling factor , and crossover rate CR = 0.2. In the case of PSO, the following values are obtained: Maximum number of generations Gmax = 500, number of individuals NP = 100, weight of the current best C1 = 1.8, weight of the global best C2 = 2.2 and speed factor . The parameters for the GA are set as: Maximum number of generations Gmax = 500, number of individuals NP = 100, mutation probability PM = 1/6, and crossover probability PC = 0.9. It is important to point out that the computational complexity of PSO and GA is similar to DE [102]. This is obtained by the number of calls to the objective function evaluation , the maximum iteration number Gmax, and the population size NP. Then, the computational complexity in a “big O” notation is given by NP Gmax). So, in all algorithms, the total evaluation number of the objective function is 50000 to make a fair comparison in the next step. (iii) Performance analysis of bio-inspired algorithms. Thirty independent runs were performed for each of the above algorithms. Each algorithm was programmed in the M language of MATLABTM, on a commercial PC with an IntelTM Core i7 @ 3.2 GHz processor, 16 GB of RAM. The descriptive statistics are presented below based on the criteria of the normalized weighted error . The average (mean), standard deviation , the best (minimum), and the worst (maximum) are specified based on the data obtained from the best individual of the last generation of each run. Tables 6-8 summarize the results obtained by each algorithm for the different criteria. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 6. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t006 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 7. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t007 The tables show that the PSO algorithm gets the best solution in terms of the objective function , as seen in Table 6. This is highlighted in boldface. Also, PSO presents the same behavior behavior concerning the individual objective function terms J1 and J2 in Tables 7 and 8. The solutions obtained with the DE techniques and GA for turn out to be very close to PSO with a percentage increment of 0.0094% for DE and 0.00000148% for GA. Additionally, the plots in Fig 9a shows a follow-up of the evolution through generations of the mean objective function of the best solutions in thirty executions (convergence graph of the mean objective function) using different optimization techniques. It is observed that while PSO converges the fastest, its overall performance through the algorithm executions is not as good on average. For completeness, Fig 9b and 9c also shows the convergence graph of the mean objective function terms using different optimization techniques, i.e., the evolution of the position error term J1 and the velocity error term J2 through generations are displayed. It is observed that the velocity error term J2 presents more oscillations before the three hundred generations than the term J1, where the J1 behavior is more stable from one hundred generations. This indicates that the achievement of minimizing the position error term J1 is more evident than the velocity error term J2. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Mean convergence graphs of the evolution in the performance function and its terms through generations in thirty executions of each algorithm. https://doi.org/10.1371/journal.pone.0325168.g009 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 8. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 9. Design variable vector obtained with different bioinspired techniques. https://doi.org/10.1371/journal.pone.0325168.t009 Even though in Fig 9a the evolution in PSO of the mean objective function shows the worse performance through generations, PSO obtained a variety of solutions that resulted in an average and standard deviation higher than those of the other two algorithms (see Tables 6–8). PSO also demonstrates a greater search space exploration (see standard deviation), given these sparse solutions, in contrast to the other two algorithms. Despite this, all algorithms can find solutions with comparable performance. In general, the three algorithms, convergence to a value in the case of the weighted target and the term J1 occurs relatively quickly after two hundred generations. In contrast, the term J2 manifests a more irregular behavior until the end of the generations. In order to know about the algorithm performance and to generalize the results, the inferential statistics [89] is applied to confirm the performance among algorithms. In particular, the Friedman test followed by Bonferroni post hoc multiple comparisons with a significance level is used. Fig 10 shows the multiple comparison test considering PSO as the control method, where the winner is the algorithm that ranks closest to the left (x-axis in the illustration). Overlapping gray lines show no discernible difference between the algorithms, meaning that no conclusions about the data distributions can be drawn because the p-values in those comparisons are greater than the selected significance level. The red line that does not overlap indicates that the p-value is below the selected significance level in algorithm outcomes, meaning that the results differ significantly. The red line shows that PSO is similar to GA, indicating that both are the most promising algorithms for this problem. PSO and GA also outperform the DE behavior. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. The Bonferroni multiple comparison post hoc test with a significance level of = 0.05. The control method is PSO. https://doi.org/10.1371/journal.pone.0325168.g010 On the other hand, Table 9 presents the vectors of resulting design variables p* corresponding to the best solution concerning the performance function of the thirty runs reported in Table 6 for all algorithms. Table 10 compares the obtained parameters previously reported in Table 2 with those obtained using different algorithms, evidencing a certain percentage of variation among them. In the case of the center of mass , the variation ranges from to ; the inertia ranges from to ; the mass mi varies between the interval ; and in the case of the friction , the variation interval is . The obtained parameters are remarkably similar across the three algorithms, and those algorithms are suitable for solving the problem. However, the parameters obtained by the PSO are used because it achieves the best performance in minimizing the objective function (albeit by a small margin). Download: PPT PowerPoint slide PNG larger image TIFF original image Table 10. Percentage values of variation of the parameters reported in Table 2 with those obtained using different algorithms in the RRR robot. https://doi.org/10.1371/journal.pone.0325168.t010 On the other hand, the values of the objective function terms with the obtained parameters by the MbPBO methodology using PSO and with the initial parameters in Table 2 are displayed in Table 11. It is observed that the parameters obtained by the proposal, the terms associated with the position error J1 and the velocity error J2 improve their performances in and , respectively, regarding the obtained performance by the initial parameters. So, significant differences in the objective function performance result from the parameters obtained by the proposed MbPBO methodology, which aids in better characterizing the dynamic model in the MiL simulation. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 11. Results in the performance of the objective function terms obtained by the MbPBO methodology using PSO with respect to the initial parameters in Table 2. https://doi.org/10.1371/journal.pone.0325168.t011 Comparative analysis with direct search methods. This section analyzes the performance of direct search algorithms in solving the optimization problem related to the robot’s characterization under the proposed MbPBO methodology. The primary goal is to evaluate the trade-offs involved in using other heuristic-based optimizers with respect to bio-inspired optimization. The selected direct search optimizers include Constrained Optimization BY Linear Approximations (COBYLA), Nelder Mead, and Powell algorithms [103]. Each algorithm is executed independently thirty times using randomly generated initial conditions. The stop criterion is based on the number of objective function evaluations. The objective function evaluation number is 50000, which is related to the same value used in the bio-inspired algorithms. The descriptive statistic is presented in Table 12. The table includes the average (mean), standard deviation (), and the minimum (best) and the maximum (worst) performance function values obtained from executions. Additionally, the number of times each algorithm finds the minimum values is presented in parentheses next to the corresponding result. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 12. Summary of descriptive statistic in for the optimization process of thirty executions with direct search methods. https://doi.org/10.1371/journal.pone.0325168.t012 Based on the results presented in Table 12, the most promising direct search optimizer is Nelder–Mead, as it achieves the lowest performance function value of . Powell’s method takes second place with , followed by COBYLA, which yields a . Notably, each algorithm reaches its minimum value only once out of thirty independent executions, leading to relatively high standard deviation values and indicating a lack of consistency across runs. The comparison between bio-inspired algorithms and direct search methods (Tables 6 and 12) reveals that bio-inspired techniques not only yield superior performance values but also demonstrate higher robustness and repeatability across multiple runs. Comparative analysis with indirect search methods. Similar to the previous section, this section analyzes the performance of indirect search methods (gradient-based algorithms) in solving the optimization problem related to the robot’s characterization. The main objective is to show the importance of bio-inspired optimization in the proposed MbPBO methodology. The selected gradient-based optimizers include Sequential Least Squares Quadratic Programming (SLSQP), the Limited-memory Broyden–Fletcher–Goldfarb–Shanno with Box constraints (L-BFGS-B), the Truncated Newton Conjugate-Gradient (TNC) and Trust-Region (T-R) algorithms [104]. The same experimental setup, namely, thirty independent runs with randomly initialized conditions and a stopping criterion (50000 objective function evaluations), is adopted as in the previous analysis. The average (mean), standard deviation (σ), minimum (best), and maximum (worst) performance function values obtained from the algorithm executions are presented in Table 13. The frequency at which each algorithm achieved the minimum value is shown in parentheses. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 13. Summary of descriptive statistic in for the optimization process of thirty executions with indirect search methods (gradient-based algorithms). https://doi.org/10.1371/journal.pone.0325168.t013 The most promising gradient-based optimizers are SLSQP and L-BFGS-B based on the results presented in Table 13 because both find solutions with the performance function value of . Nevertheless, SLSQP stands out by reaching this optimal value seventeen times, indicating greater consistency. The Trust-Region (T-R) takes second place but only finds two times the performance function value . The worst gradient-based optimizer is TNC because it exhibits the largest performance function value and only finds the sub-optimal solution once out of thirty executions. It is also observed that the most promising gradient-based optimizers (SLSQP and L-BFGS-B) locate the best solution obtained with the most promising bio-inspired algorithms (PSO and GA in Table 6). Nevertheless, the gradient-based algorithms are very sensitive to the initial conditions, providing a diverse set of local solutions, as described in the standard deviation of Table 13. Although certain gradient-based algorithms (such as SLSQP and L-BFGS-B) exhibit moderate consistency achieving the optimal solution in around half of thirty executions, the advantage of bio-inspired algorithms (see Tables 6 and 13) lies in their global search capabilities and robustness to problem characteristics such as non-smoothness noise, and complex constraint landscapes. While gradient-based methods may converge efficiently when the search space is well-behaved (e.g., smooth and convex), they depend highly on initial conditions. They may also become trapped in local minima in more irregular or multimodal functions. On the other hand, bio-inspired algorithms incorporate exploration mechanisms to avoid local optima and search for in a wider design space region. These advantages are observed in the lower standard deviation of Table 6, indicating greater consistency in the obtained results. Thus, even if a subset of gradient-based methods performs similarly on this specific problem, bio-inspired algorithms present a more versatile and reliable option across diverse optimization problems where gradient information, constraint inclusion, multi-objective formulations, or parallel implementation is challenging. This makes bio-inspired algorithms particularly well-suited for the MbPBO methodology, which benefits from robust global search strategies and adaptability to complex problem landscapes. With these findings, the MbPBO methodology can be extended to using hybrid techniques known as memetic algorithms. The memetic algorithm combines the global exploration capabilities of the bio-inspired algorithms with the fine search (exploitation capabilities) of gradient-based techniques. This hybrid extension represents a possible future work in the proposed methodology. Stage V: Validation of the obtained RRR robot model and its use in the MiL simulation Once the comparison described in the step iii) of Stage IV has been carried out, and the design variable vector resulting in the minimum value of the weighted objective function (p* obtained with PSO reported in Table 6) has been obtained, the last stage of the proposed MbPBO methodology is the comparative evaluation concerning the similarity of the system model’s states with those of the plant. So, the first step analyzes the model obtained in open-loop control with the reference signal from Stage II. In addition, this stage also focuses on assessing the parameterized model’s ability to accurately replicate the behavior of the testbed platform in a wide range of applications to test different MiL simulation environments. This evaluation essentially tests the model’s predictive capability with different control strategies and monitoring scenarios by comparing position and velocity states between the model and the actual testbed platform during specific tasks. So, the second step of this stage analyses the model obtained in closed-loop for the regulation and tracking tasks. In this stage, graphic comparisons and the Pearson correlation coefficient analysis between experimental and simulation signals are performed to assess the similarity or variability between the system model’s states and the plant’s numerically. The steps of this stage are described below. (i) Model validation in open-loop control. Figs 11 and 12 show the Cartesian angular position and angular velocity behavior of the RRR robot obtained from the approximated model (signal displayed in solid blue line) and the real experiment (signal represented by the dotted magenta line). The Cartesian position in the real experiment and the simulated plant is represented by and , respectively. Similarly, the angular positions are defined by q1, q2, q3 for the real experiment and , , for the behavior of the obtained model. The angular velocities are set as , , for the real experiment and , , for the simulated plant. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Robot’s end-effector behavior in the Cartesian plane with the solution obtained by PSO in simulation, against the real states of the experiment when system states are excited by the open-loop control signals. https://doi.org/10.1371/journal.pone.0325168.g011 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Angular position and velocity behaviors of the RRR robot when the system states are excited by the open-loop control signals. Real experiment (q1, q2, q3, , , ) vs. approximated model (, , , , , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g012 For a better numerical appreciation of the results observed in Figs 11 and 12, Table 14 presents the correlation coefficient between the real measurements and the results of the approximated model for each state previously examined. The coefficient is a statistical measure that indicates the strength and direction of the linear relationship between two variables (Pearson correlation). So, this measure would indicate the similarity (linear relationship) between the data obtained by the model and the experimental data. A coefficient of -1 would indicate a perfect inverse correlation, i.e., with the same input data, the model provides an output with the same magnitude but in an opposite direction than the output of the real system. The value of 0 indicates no correlation, which means that the model data and the experimental data do not exhibit a linear relationship. Hence, the model could not be helpful to describe the experimental (real) behavior. The value of 1 indicates a perfect positive correlation (strong correlation) between the variables, i.e., the model data matches perfectly to the experimental data with the same input data, meaning that the model is useful for describing the real behavior. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 14. Correlation coefficient between the robot signals and the approximated model behavior when the system states are excited by the open-loop control signals. https://doi.org/10.1371/journal.pone.0325168.t014 The results show suitable relationships between the actual behavior and the approximated model with a average correlation of 0.9015 in the angular space and 0.8969 in the Cartesian space. In particular, the average correlations for the angular position and angular velocity are 0.9656 and 0.8375, respectively. In particular, the velocity states exhibit lower correlation coefficients than the position states concerning the actual experimental behavior, attributed to the noisy measurement in the velocity states. The noise in the identification tasks is a significant challenge, often necessitating well-tuned filtering techniques, advanced measurement instrumentation, or meticulous trajectory planning [77]. (ii) Model validation in closed-loop. In this step, two control strategies are also employed to test two different MiL simulation environments: 1) A PID control strategy for regulation at a single point in the joint space and 2) A PD+G control strategy for following a smooth trajectory in the joint space. The comparisons between the model and the testbed platform, which include graphical representations and corresponding correlation coefficients, form a thorough validation process. This process instills confidence in the model’s accuracy by evaluating how well the model’s predicted states match the actual measured states. The comprehensive analysis aims to validate the model’s effectiveness in capturing the dynamics and behaviors of the testbed platform under different operational conditions. This study’s primary validation measures are the comparisons between measured and estimated position and velocity states. However, additional validation measures, as suggested by [93], could include the prediction of the necessary torque in the actuators to achieve a specific task or how far the original model parameters (perhaps required from the manufacturer) lie against the predicted values (confidence interval). These validation measures typically require more complex and potentially expensive means of data collection and analysis. Therefore, while these measures are valuable for a comprehensive validation of the model’s predictive capability, they are not included in the scope of this particular study. The focus here remains on validating the positional and velocity states as a primary indicator of the model’s accuracy in representing the behavior of the testbed platform. Regulate a point in the joint space with the PID control strategy Table 15 presents the positions to be reached in the joint space for each degree of freedom and the PID controller parameters used for each joint, applicable to both the model and the laboratory testbed platform. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 15. PID controller parameters and the desired position in the joint space used in each d.o.f. for the validation of the obtained RRR robot model in the regulation task. https://doi.org/10.1371/journal.pone.0325168.t015 Fig 13 illustrates the workspace of the manipulator , where the thick black lines represent the robot’s links with the end effector at one end, the solid blue line shows the measured real values , and the thick magenta dotted line represents the predicted values . Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 13. Robot’s end-effector behavior in the workspace of the obtained model against the real experiment for the regulation task using PID controller. Black continuous lines indicate the robot links. https://doi.org/10.1371/journal.pone.0325168.g013 Concerning the controller performances, Fig 14a, 14c, 14e presents the actual and approximated position states for the regulation task. Likewise, Fig 14b, 14d, 14f includes the real and approximate velocity states for the same task. The control signals applied to the model and testbed platform are presented in Fig 15. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 14. Angular position and velocity behaviors of the RRR robot in the regulation task using PID controller. Real experiment (q1, q2, q3, , , ) vs. approximated model (, , , , , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g014 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 15. Control signals () for each d.o.f. obtained from the regulation PID control experiment. Real experiment (u1, u2, u3) vs. approximated model (, , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g015 Table 16 presents the respective correlation coefficients between the states of the testbed platform and those of the model for this experiment. Based on this table and the corresponding Fig 14 in the regulation task with the PID controller, it is observed that the overall correlations in the Cartesian space of the robot are very confident in the interval , meanwhile the angular position space presents a correlation in the interval and the angular velocity space provides a correlation in the interval . With this information, the approximated model successfully behaves as the robot’s end-effector behavior of the testbed platform. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 16. Correlation coefficient between the obtained model signals and the robot signals in the regulation task with the PID controller (closed loop). https://doi.org/10.1371/journal.pone.0325168.t016 Smooth tracking of a path with the PD+G control strategy In this test, the trajectory consists of tracking three points in the joint space, executed in two different quadrants of the Cartesian space. Those points are joined by a third-order Bezier polynomial, resulting in a smooth trajectory. Table 17 displayed the points reached by this trajectory and the controller parameters used for the PD+G control in both the modeling and experimental setups. Fig 16 presents a view of the manipulator’s workspace (Xw,YW), showing the comparison of both the real measured behavior in solid lines and the predicted one in dotted lines, following the convention of previous plots, but this time for the tracking task. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 16. Robot’s behavior in the workspace of the obtained model against the real experiment for the tracking task using PD+G controller. Black continuous lines indicate the robot links. https://doi.org/10.1371/journal.pone.0325168.g016 The plots in Fig 17a, 17c, 17e present the real and approximate position states of for the tracking task. Similarly, the real and approximate velocity states are depicted in Fig 17b, 17d, 17f for the same trajectory tracking task. The control signals of the model and the testbed platform can also be observed in Fig 18. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 17. Angular position and velocity behaviors of the robot in the tracking task using PD+G controller. Real experiment (, , , , , ) vs. approximated model (, , , , , ) with the solution obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g017 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 18. Control signals () for each d.o.f. obtained from the tracking PD+G control experiment. Real experiment (u1, u2, u3) vs. approximated model (, , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g018 Table 18 presents the respective correlation coefficients between the states of the testbed platform and those of the model for this experiment. Based on this table and Fig 17 in the tracking task with the PD+G controller, it is observed that the overall correlations in the Cartesian space of the robot and the angular position space are very confident with a value around one. Meanwhile, the angular velocity space provides a correlation in the interval . With this information, the approximated model successfully behaves as the robot’s end-effector behavior of the testbed platform. Highlights of the validation results. The model validation in the open loop control indicates a high average correlation among all joint states of 0.9015, indicating a high representation of the real system with the obtained parameterized model. In particular, the average correlation in the angular position, angular velocity, and Cartesian position is 0.9656, 0.8375, and 0.8969, respectively. On the other hand, the tests in the MiL simulation using the obtained parameterized model for the regulation and the tracking control problems (closed loop control) shows a high average correlation among all joint states in both control problems of 0.9022, revealing the model’s effectiveness for general-purpose use in MiL tests. In particular, the correlations given for the Cartesian space is in the interval and for the angular position space is in the interval providing an average correlation of 0.9906 and 0.9878, respectively. Meanwhile, the correlation in the angular velocity space is in the interval with an average correlation of 0.8165. The slight decrement in the angular velocity correlation is attributed to the estimation of the angular velocity by using finite difference because it produces noise in the velocity measure and also in the trade-off selected in Stage IV. Application of the MbPBO methodology to the RFSEA The main purpose of this section is to provide useful information about the applicability and the flexibility of the MbPBO methodology in different electro-mechanical systems. In this case, the proposed MbPBO methodology is applied to the model parameterization of the Reaction Force-sensing Series Elastic Actuator (RFSEA), confirming the similarities in the behavior of the obtained model with respect to the real system’s behavior. In what follows, some aspects of the methodology are highlighted. The RFSEA is a series elastic linear actuator [105] where the tip of actuation (the ball screw) involves an active and passive movement through the displacement of the pulley-belt transmission and spring, respectively. The Series Elastic Actuator (SEA) has been widely used in robotic leg [106], biomechanical devices [107–109], mechatronic systems [110], lower and upper limb rehabilitation devices [111,112], among others. The most important aspects of the MbPBO methodology are highlighted next. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 17. PD+G controller parameters and the desired positions in the joint space joined by third-order Bezier polynomials for the validation of the obtained RRR robot model in the tracking task. https://doi.org/10.1371/journal.pone.0325168.t017 Fig 19 represents the RFSEA system divided into the three most important parts of movements in different colors: The motor’s belt-pulley transmission, the spring movement, and the linear motion of the ball screw. The variables associated in the schematic diagram of the RFSEA are: (A) For the belt-pulley transmission: the motor’s angular position , the input torque (motor’s torque), the output torque , the motor’s current ia, the motor’s inertia Jm, the motor’s damping coefficient Bm and the radii of the drive pulleys R1 and R2. (B) For the spring movement, the spring’s linear displacement xs, the spring’s stiffness coefficient Ks, the spring’s mass Ms, the spring’s damping coefficient Bs, the pulley radius re for the linear spring’s displacement and its angular displacement θe. (C) For the ball screw’s linear motion, the ball screw displacement xl, the ball screw’s damping coefficient Bl, the ball screw’s mass Ml, the screw pitch l and the load’s reaction force Dl. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 19. RFSEA system. https://doi.org/10.1371/journal.pone.0325168.g019 The dynamic equations of the RFSEA consist of the mechanical (37)-(39) and electrical (40) dynamics [113–115], where is the transmission ratio of the ball screw with as the belt-pulley transmission ratio obtained when the larger pulley with radius R1 is the output and the smaller pulley with radius R2 is the input. (37)(38)(39)(40) The dynamic coupling of the mechanical and electrical parts of the RFSEA assuming can be expressed in state space with the input signal and the reaction force Dl. The state-space representation of the RFSEA dynamic model results in (41), where , ,  +  ,  +   +  , − , ,  +  , − ,  +   +  , , , , , and . (41) where: Download: PPT PowerPoint slide PNG larger image TIFF original image Table 18. Correlation coefficient between the obtained model signals and the robot signals in the tracking task with the PD+G controller (closed loop). https://doi.org/10.1371/journal.pone.0325168.t018 The Stage I is finished with the obtained model of the RFSEA. The characterization of the system’s behavior in Stage II is carried out in the same manner as described in the MbPBO methodology. The specific set of parameters to be found and the objective function are selected in the optimization process for Stage III. In this case, the selected design parameters vector is related to the damping coefficients with the terms of the objective functions as , where is the estimated state. In addition, the normalization is done with the corresponding maximum values of the objective function terms obtained in the procedure (, , , 10−5, ). The values of the rest model parameters associated to the RFSEA are described through the datasheet of components, or by directly using measuring instruments. Following Stage IV of the proposed model parameterization approach, the selected weights are , giving the same importance to all terms in the objective function. In Table 19, the parameters obtained using the MbPBO optimization process are indicated with an asterisk, along with the remaining parameters employed. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 19. Parameters of the RFSEA. https://doi.org/10.1371/journal.pone.0325168.t019 With the obtained parameters shown in Table 19, the parameterized model is completed, and the validation of the obtained parameterized model using an open-loop controller is carried out as described in Stage V. In this case, the drivers to control the RFSEA motor are configured in voltage mode, so the reference signals in the open-loop controller are related to the voltage and as a consequence the open-loop control signal is u = ref (a cascade controller is not applied in the RFSEA case). The behaviors of the state vectors for the real experiment and the approximated model with the same persistent excitation control signal are given in Fig 20. The bar above the variable indicates the behavior of the obtained model. The Pearson correlation coefficients of the signals related to the screw displacement, the screw velocity, and the current present a suitable correlation of and 0.9226, respectively. Those correlations indicate that the obtained model represents a similar behavior to the real system. The worst correlation coefficients are given by the spring position and spring velocity ( and , respectively). This is attributed because those measures present small magnitudes (see Fig 20b, 20d) due to spring states are not excited in the experiments by applying the external force Dl and are mainly affected by the vibrations of the testbed platform. So, those spring measures are not representative in the study. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 20. State behaviors of the obtained RFSEA model against the real experiment signals when the system states are excited by the open-loop control signals. Real experiment (xl, xs, , , ia) vs. approximated model (, , , , ). https://doi.org/10.1371/journal.pone.0325168.g020 On the other hand, the validation of the obtained parameterized model using a closed-loop controller for testing two different MiL simulation environments is carried out as described in Stage V. Table 20 provides the respective correlation coefficients between the states of the RFSEA testbed platform and those of the model for the regulation and tracking tasks. The representative figures for those signals are presented in Fig 21. Based on this table, it is observed that the average correlation in the regulation and tracking tasks in the ball screw’s linear displacement is 0.9743 and in the ball screw’s linear velocity is 0.8356, indicating high correlation values, i.e., the model data are highly related to the experimental data using the same input data, meaning that the model is useful for describing the real behavior of those signals. In the motor’s current, the average correlation is 0.4263, indicating a moderate correlation value and that the obtained model behavior could present deviations from the real behavior, which is mainly attributed to a low current transducer sensitivity. In the case of the spring position and its velocity, these signals present an average correlation of −0.1749 and 0.0446, respectively. Those correlation values indicate that there is no relationship between the model data and the real behavior. Nevertheless, those weak correlations are attributed to the spring signals not being excited in either the model or the experimentation. So, their movements in the experimental data have small magnitudes (see Fig 21c, 21d, 21g, 21h) mainly given by the vibrations of the testbed platform which are not considered in the model. For instance, the movement of the spring is not perturbed by an external load in both experiments. Then, the spring position and velocity correlations are not representative in the study. So, the results obtained in the MiL simulation environments indicate that the obtained RFSEA model can describe, to a large extent, the real behavior. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 21. Behavior of the RFSEA’s states in the regulation (left) and tracking (right) tasks. Real experiment (xl, xs, , , ia) vs. approximated model (, , , , ). https://doi.org/10.1371/journal.pone.0325168.g021 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 20. Correlation coefficient ρ between the obtained model signals and the system signals in the regulation and tracking tasks with the PID controller in the RFSEA. https://doi.org/10.1371/journal.pone.0325168.t020 General discussion Based on the outcomes achieved at various stages of applying the MbPBO methodology, the following is highlighted: Bio-inspired techniques are widely regarded as effective alternatives for modeling and identification tasks. Based on a comparative analysis involving PSO, GA, and DE algorithms in solving optimization problems, the PSO algorithm consistently converges toward the best solution across multiple executions. A persistent pattern is observed in the percentage variation between the nominal values of the experiment and the solutions obtained using the three algorithms. These findings underscore the importance of the MbPBO methodology in accurately parametrizing models to achieve behavior closer to real systems. Despite unmodeled uncertainties in the investigation, such as noise signal, backlash between gears, and hysteresis, among others, the simulated states exhibit a high average general correlation of 0.9019 considering open-loop control and various tasks and control strategies, including regulation and tracking tasks (MiL tests) in the joint space. Particularly, this model achieves an average correlation of 0.9015 and 0.9022 in open and closed-loop control, respectively. This high positive correlation indicates that the behavior of the obtained model is strongly related to the behavior of the real robotic system. This means that the obtained model can be used to predict the system’s behavior. The MbPBO methodology applied to a different electro-mechanical system (RFSEA) reveals the flexibility of the proposal. In this case, considering the states that were excited in the RFSEA, the correlations in open and closed loop control are 0.9305 and 0.7454, providing a high average general correlation of 0.8379. This correlation indicates a suitable approximated model. Obtaining a representative system model that closely approximates reality offers several advantages. Firstly, it eliminates the need for repetitive tests on the actual system when maximizing task performance, such as in controller tuning. Secondly, the methodology includes model validation within the MiL approach, accelerating the implementation of model-based design (MiL-SiL-HiL). Stage I: Determination of the base model for the RRR robot and its key parameters The considered dynamic system is an autonomous system; as such, its physical parameters are invariant in time. This characteristic makes it suitable for creating a parametric model where the parameters can be identified offline [93]. Additionally, the model is also proposed as a non-conservative system. The frictions between links present a challenge objective for the optimal parameterization proposed by the MbPBO methodology because specifying these parameters is difficult for any manufacturer. A dynamic model of the system is chosen for the primary mathematical modeling effort, where the relationship between the movement and the forces that produce it is explicitly described. For its formulation, it was decided to use the Euler-Lagrange formalism, whose principle is the balance of the kinetic and potential energy of the components of the system. The basis of this model is the Lagrangian equation (13) where and P(q) represent, respectively, the kinetic and potential energy of the links. (13) Likewise, when considering a non-conservative system where energy is dissipated, the Rayleigh dissipation function is included in this formulation. This function models these friction forces as viscous dampers in the term D shown in (14). (14) The term bi in this last equation represents the damper coefficients modeled for the i–th link. Including this dissipation function in the model, the corresponding Lagrange equations for the modeling of non-conservative systems are expressed in equation (15). (15) The model obtained from [94] developed with the Euler-Lagrange energy method was abbreviated in terms to assume the limitations of the simulator and promote shorter calculation times. The latter process considers the third aspect of this first stage related to the model reduction. The effect of viscous friction for each link is represented as , and likewise, the vector τ=[τ1,τ2,τ4] represents the input vector of torques. The closed form of the robot manipulator’s dynamic model is described in (16), where represents the inertia matrix, is the matrix of Coriolis and centrifugal forces, G∈ℝ3 is the gravity vector, and corresponds to the vector of torques produced by friction. (16) where: (17)(18)(19)(20) with: (21)(22)(23) The RRR robot dynamics (16) in the state vector with the control signal is expressed in a simple form as . The following is given concerning the steps (i) and (ii) of the Stage I. (i) Selection of the parameters to be estimated. It is then defined for this system (based on Fig 3, and the previously presented model), the vector p (24) as the set of parameters to be obtained in the identification of the RRR robot. The vector in (24) includes 12 elements to be estimated, where denotes the vector of distances between the axis of movement of each link to its center of mass, is the vector of inertia, is the mass vector, and refers to the viscous friction vector. The parameters in p are essential for identifying the RRR robot’s dynamic model using the MbPBO methodology. (24) (ii) Initial values proposed for the parameters. In the process of parameter identification, it is beneficial to begin with prior knowledge of values close to the parameters to be estimated. This approach allows for a constrained search, enhancing the efficiency of model refinement. The initial values for the parameter vector (24), grouped in , are sourced from [94]. In that study, the SolidWorks mechanical modeling program was utilized to estimate the parameters of the robot, such as mass, inertia, and mass center of links, as detailed in Table 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Proposed initial dynamic and kinematic parameters for the RRR robot. https://doi.org/10.1371/journal.pone.0325168.t002 In this case, neither the manufacturer nor a CAD program can provide the friction parameter . This parameter pertains to the coefficient of friction for dry surfaces, also known as dry friction, which varies depending on the materials in contact. The coefficient of friction between two dry metallic surfaces typically ranges from 0.15 - 0.8 [75]. Therefore, the known range of friction coefficients for metallic surfaces and factors that could result in the highest friction coefficient (like oxidation and high contact pressure, to which the robot joints could be exposed) were considered when proposing the friction coefficient values for the parameter estimation process. In this particular case, one of the worst cases is selected as the initial value for the friction coefficients. So, those values are initially set as 0.9. (i) Selection of the parameters to be estimated. It is then defined for this system (based on Fig 3, and the previously presented model), the vector p (24) as the set of parameters to be obtained in the identification of the RRR robot. The vector in (24) includes 12 elements to be estimated, where denotes the vector of distances between the axis of movement of each link to its center of mass, is the vector of inertia, is the mass vector, and refers to the viscous friction vector. The parameters in p are essential for identifying the RRR robot’s dynamic model using the MbPBO methodology. (24) (ii) Initial values proposed for the parameters. In the process of parameter identification, it is beneficial to begin with prior knowledge of values close to the parameters to be estimated. This approach allows for a constrained search, enhancing the efficiency of model refinement. The initial values for the parameter vector (24), grouped in , are sourced from [94]. In that study, the SolidWorks mechanical modeling program was utilized to estimate the parameters of the robot, such as mass, inertia, and mass center of links, as detailed in Table 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. Proposed initial dynamic and kinematic parameters for the RRR robot. https://doi.org/10.1371/journal.pone.0325168.t002 In this case, neither the manufacturer nor a CAD program can provide the friction parameter . This parameter pertains to the coefficient of friction for dry surfaces, also known as dry friction, which varies depending on the materials in contact. The coefficient of friction between two dry metallic surfaces typically ranges from 0.15 - 0.8 [75]. Therefore, the known range of friction coefficients for metallic surfaces and factors that could result in the highest friction coefficient (like oxidation and high contact pressure, to which the robot joints could be exposed) were considered when proposing the friction coefficient values for the parameter estimation process. In this particular case, one of the worst cases is selected as the initial value for the friction coefficients. So, those values are initially set as 0.9. Stage II: Characterization of the RRR robot’s behavior. To characterize the system behavior, identification experiment on the RRR robot is developed. The experiment used open-loop control signals (open-loop external signal into the outer loop of the cascade system) along with their respective periods, following persistent excitation signals recommended in [95] for this purpose. Each reference signal refi to the inner loop controller for each degree of freedom (d.o.f.) has a square wave signal. These signals varied between the maximum applied current that the actuators could handle. Specific periods T were interleaved within each reference signal, ensuring an experimentation time with a sampling time of . The experiment induced various behaviors in the robot’s states by applying these signals ref in an open-loop control configuration. In this case study, we focused on monitoring the positions and velocities of the robot’s links . A representative diagram of this experimental setup can be seen in Fig 5. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Schematic diagram of the open-loop control experiment conducted on the experimental platform to excite the system states. https://doi.org/10.1371/journal.pone.0325168.g005 The individual periods applied to the three degrees of freedom are: , and . The graphical representation of the reference signals ref1, ref2, and ref3 in the outer loop controller to provide the reference in the inner loop controller of the robot is visualized in Fig 6. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Reference signals used in the open-loop controller on the experimental platform to excite the corresponding system states. https://doi.org/10.1371/journal.pone.0325168.g006 At the end of this stage, a data vector of position states defined as is obtained, which is directly retrieved from the robot’s encoder over the period of (resulting in a state measurement vector of 2600 samples). These measurements are graphically represented with blue lines in the plots of Fig 7a, 7c, 7e. The velocity state, defined as q˙=[q1˙,q2˙,q3˙]T∈R3, is estimated using numerical differentiation in MATLAB. The plots of Fig 7b, 7d, 7f illustrate the estimated velocity values. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Graphs showing the position state vectors q and the estimated velocity state vectors of the real system when the system states are excited by the open-loop control signals. https://doi.org/10.1371/journal.pone.0325168.g007 Stage III: Establishment of a mathematical programming problem for the RRR robot characterization. With the information from the previous stage, the current reference signal ref established in the previous stage is ideally transformed into torque input for only simulation purposes by the equation [96], where km = 1.9839 is the torque constant. So, the control signal is applied in simulation to the degrees of freedom in the robot model defined in Stage I. Its numerical solution provides us with a vector of estimated positions and velocities, defined as and , respectively. The error vector (25) is then defined between the actual position and the estimated position , as well as, the error vector of the velocity states between the experimental data and those associated with the model (26). Thus, the position and velocity errors of the i-th link are represented by the vectors ei and , respectively. (25)(26) So, in the optimization problem, the goal is to minimize the error between the signals obtained from the experiment in Fig 7 and the output of the simulated system in the time domain (error dynamics). Therefore, it is imperative to minimize errors in the states. For this purpose, the criterion of minimizing the integral of the square error (ISE) is chosen, as it is most suitable for this purpose [97]. Once the error vectors have been defined, the objective function terms are described in the space of positions and velocities in (27) and (28), respectively. (27)(28) Following the indication given in Stage III of the MbPBO methodology description, the multi-objective nature of the problem is transformed into a single-objective problem by implementing a weighted sum approach. Here, dimensionless values and represent the weights of the position and velocity errors, respectively. The objective function is minimized based on a normalized weighted expression . For this particular case study, the objective function can be defined as in (29), where only two terms are included to provide the same importance for the position errors and other priority for the velocity errors. (29) In (29), the scalar values and are quantities that represent the maximum possible values of J1 and J2 that will produce a dimensionless and normalized value of (ranging between 0 and 1). Those maximum values are obtained by only minimizing one term in the optimization process, i.e., by only minimizing J1 (27) or J2 (28), and the obtained optimal design variable vector or is evaluated in the other terms. So, in the particular case, the maximum value is calculated by evaluating the term J2 with the obtained vector (when the J1 term is minimized), i.e., . In a similar fashion, the maximum value is given by minimizing J2, obtaining the optimal design variable vector and calculated the maximum value as . The specific values of and are provided in the next Stage IV. The dynamic constraints (35) related to the system behavior are included in the optimization problem to bind the problem solution to the dynamic behavior of the RRR robot (16). A tolerance range is established around the initial (base) value of each variable pi to limit the search space in the optimization process. The increment or reduction factor is denoted as , where represents the desired percentage change over the base value shown in Table 2. The chosen increment/reduction factor values impose practical boundaries in the design variables to ensure that the search space in the optimization process presents valid solutions. The optimization process with these bounds is also more efficient because it does not waste time looking for solutions where the system can not operate. An increase and decrease factor of at most thirty percent is recommended for this purpose. Equations (30) to (33) describe such a limits () and define how the bound constraints are given. The selected upper and lower bound ranges ( ) are obtained from their initial (base) values and the increment/reduction factor given in Table 3. So, the bounds of the design variable vector (30)-(33) are the other constraints taken into account in the optimization problem. In (29), the scalar values and are quantities that represent the maximum possible values of J1 and J2 that will produce a dimensionless and normalized value of (ranging between 0 and 1). Those maximum values are obtained by only minimizing one term in the optimization process, i.e., by only minimizing J1 (27) or J2 (28), and the obtained optimal design variable vector or is evaluated in the other terms. So, in the particular case, the maximum value is calculated by evaluating the obtained vector when the J1 term is minimized, i.e., . In a similar fashion, the maximum value is given by minimizing J2, obtaining the optimal design variable vector and with this information can compute the maximum value as . The specific values of and are provided in the next Stage IV. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. Percentage rate in the design variable changes and the design variable bounds used in the optimization process. The specific quantities are indicated with . https://doi.org/10.1371/journal.pone.0325168.t003 The dynamic constraints (35) related to the system behavior are included in the optimization problem to bind the problem solution to the dynamic behavior of the RRR robot (16). A tolerance range is established around the initial (base) value () of each variable pi to limit the search space in the optimization process. The increment or reduction factor is denoted as , where represents the desired percentage change over the base value shown in Table 2. The chosen increment/reduction factor values impose practical boundaries in the design variables to ensure that the search space in the optimization process presents valid solutions. The optimization process with these bounds is also more efficient because it does not waste time looking for solutions where the system can not operate. An increase and decrease factor of at most thirty percent is recommended for this purpose. Equations (30) to (33) describe such a limits () and define how the bound constraints are given. The selected upper and lower bound ranges () are obtained from their initial (base) values given in Table 2 and the increment/reduction factor given in Table 3. So, the bounds of the design variable vector (30)-(33) are the other constraints taken into account in the optimization problem. (30)(31)(32)(33) Relating the above with the formal approach, the problem consists of finding the vector of design variables p that minimizes the weighted objective function (34), subject to the dynamic behavior of the system model (35), and the constraints on the design variables (36). In this case, the equality constraints are implicitly satisfied by the solution of the model dynamics. (34) subject to: (35)(36) The technique applied to obtain the minimum among these errors will yield the vector of optimal design variables p*, where and n = 12 (24). The optimization techniques used to solve this problem will be addressed in the following Stage IV. Stage IV: Solving the RRR robot’s characterization problem using bio-inspired algorithms It is important to point out that both terms of the objective function (29) present trade-offs for fulfilling both criteria. This means that minimizing both errors is not possible, i.e., the relationship between them turns out to be inverse; they behave like two opposing objectives where only minimizing the first term J1 implies maximizing the second one J2 and vice versa. In order to limit the weight values and to , it is important to normalize the terms (to set the values of such terms in the interval between ) of the objective functions. So, the terms max(J1) and max(J2) of (29) are obtained by individually minimizing the terms J2 and J1, respectively, and using the same constraints in (35)-(36), as described in the previous stage. In this part of the methodology, it is recommended to set several independent runs (for instance, thirty times [89]) for solving both optimization problems due to the stochastic behavior of the algorithm. DE/rand/1/bin is used in this case, but other frequently used bio-inspired optimizers can also be used. The best execution, the one that has the minimum value of the term, is selected from the thirty runs. The result that further minimizes the term J1 provides the objective function vector . Meanwhile, the result that further minimizes the term J2 provides the objective function vector . The values in boldface are the possible extremes sought, i.e., and . For practical simplicity, the values that are considered to carry out the normalization in the equation (29) are and . Once the maximum values of the terms in the optimization problem are obtained, the three steps considered in this stage are described below. (i) Trade-off selection. As mentioned in the description of the MbPBO methodology, the multi-objective problem is transformed into a single-objective problem by a weighted sum approach. A set of non-dominated solutions must be generated by establishing different weights in the optimization problem (34)-(36), forming the approximate Pareto front in the performance function space. The Pareto front helps the designer make decisions on the most suitable trade-off for the specific application. Once the terms of the objective function are normalized, new experiments are carried out to find the best trade-offs in the terms J1 and J2 (quadratic errors of position and velocity). Then, fourteen combinations of weights and in the interval are tested. The first two combinations find the solution of the Pareto front without a trade-off, i.e., by setting and . The next nine combinations consider uniform intervals from to subject to  +  . The three final weights are assigned according to the best promising weight interval found in the previous experiment. The promising weight interval is visualized once the previous experiments with different weights are carried out. Thirty independent executions of the bio-inspired algorithm are performed by establishing different weights in the optimization problem (34)-(36) and using the same DE algorithm with the same parameters given above. This experiment aims to identify the combination of weights that results in a helpful trade-off solution by including the best solution (the solution that further minimizes the performance function from the algorithm executions) found in each combination. Table 4 shows the chosen weight values, the weighted normalized objective function , and its weighted non-normalized terms J1 and J2. The values in boldface indicate the minimum ones in each column. As the uniform distribution of weights does not guarantee uniform distribution of the found Pareto solutions [98], and even the bio-inspired optimizer can stagnate in local solutions [83], the obtained solutions are sorted according to non-dominance verification. The Pareto solutions [83] of results from Table 4 are displayed in Table 5 and the graphical representation of the Pareto front with the normalized terms and is shown in Fig 8. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Normalized Pareto front obtained by using different weights in the optimization process. https://doi.org/10.1371/journal.pone.0325168.g008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Values of the objective function and its non-normalized terms , obtained by using different weights in the optimization process. https://doi.org/10.1371/journal.pone.0325168.t004 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 5. Set of Pareto solutions obtained by the results of Table 4. https://doi.org/10.1371/journal.pone.0325168.t005 Once the Pareto front is obtained, the decision maker’s preference provides one useful trade-off for the application. In this case, the most preferred solution is the Pareto solution with the minimum value in the term J1 and with a suitable value in the term J2. Based on Table 5 and Fig 8, the combination of weights , stands out for its ability to further minimize the trade-off in the performance function and the position quadratic errors, also providing an acceptable velocity quadratic error. The selected weights are employed in the next step to analyze different bio-inspired algorithms to optimize the RRR robot’s characterization problem. (ii) Algorithm selection. The choice of the DE/rand/1/bin algorithm for the experiments performed previously is because it is one of the most widely used alternatives in the reviewed literature [99] (Fig 1). However, the MbPBO methodology description in this work suggests using two more techniques to solve the optimization problem with the selected obtained weights (, ). Therefore, the next frequently used algorithms in the reviewed literature are selected due to their adaptability, population-based strategy, robust performance, and simplicity [85–87]. In particular, the Particle Swarm Optimization (PSO) [39] with the best global topology and the Genetic Algorithm (GA) in its generational variant [100] are also selected to optimize the RRR robot’s characterization problem described before (34)-(36). As mentioned in the methodology description, it is only possible to determine which algorithm presents the best results for a specific problem once it is experimentally proven [88]. Therefore, given the previous algorithms, the present work proposes a comparison based on the results obtained to identify which bio-inspired technique generally performs best when solving the characterization problem. The algorithm parameters must be tuned either “by hand” through systematic manual settings, “by analogy” through the information suggested in the state-of-the-art for similar problems, or using self-tuning mechanisms. The process “by hand” (where algorithm parameters are varied sequentially and ordered to select the best configuration that solves the problem) is also suggested in [101] for algorithms like DE. Therefore, the latter is the process chosen in this work to tune such parameters. The set of parameters presented below was obtained from this trial and error procedure, where the final tuning values for the DE/rand/1/bin algorithm are as follows: Maximum number of generations Gmax = 500, number of individuals NP = 100, scaling factor , and crossover rate CR = 0.2. In the case of PSO, the following values are obtained: Maximum number of generations Gmax = 500, number of individuals NP = 100, weight of the current best C1 = 1.8, weight of the global best C2 = 2.2 and speed factor . The parameters for the GA are set as: Maximum number of generations Gmax = 500, number of individuals NP = 100, mutation probability PM = 1/6, and crossover probability PC = 0.9. It is important to point out that the computational complexity of PSO and GA is similar to DE [102]. This is obtained by the number of calls to the objective function evaluation , the maximum iteration number Gmax, and the population size NP. Then, the computational complexity in a “big O” notation is given by NP Gmax). So, in all algorithms, the total evaluation number of the objective function is 50000 to make a fair comparison in the next step. (iii) Performance analysis of bio-inspired algorithms. Thirty independent runs were performed for each of the above algorithms. Each algorithm was programmed in the M language of MATLABTM, on a commercial PC with an IntelTM Core i7 @ 3.2 GHz processor, 16 GB of RAM. The descriptive statistics are presented below based on the criteria of the normalized weighted error . The average (mean), standard deviation , the best (minimum), and the worst (maximum) are specified based on the data obtained from the best individual of the last generation of each run. Tables 6-8 summarize the results obtained by each algorithm for the different criteria. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 6. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t006 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 7. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t007 The tables show that the PSO algorithm gets the best solution in terms of the objective function , as seen in Table 6. This is highlighted in boldface. Also, PSO presents the same behavior behavior concerning the individual objective function terms J1 and J2 in Tables 7 and 8. The solutions obtained with the DE techniques and GA for turn out to be very close to PSO with a percentage increment of 0.0094% for DE and 0.00000148% for GA. Additionally, the plots in Fig 9a shows a follow-up of the evolution through generations of the mean objective function of the best solutions in thirty executions (convergence graph of the mean objective function) using different optimization techniques. It is observed that while PSO converges the fastest, its overall performance through the algorithm executions is not as good on average. For completeness, Fig 9b and 9c also shows the convergence graph of the mean objective function terms using different optimization techniques, i.e., the evolution of the position error term J1 and the velocity error term J2 through generations are displayed. It is observed that the velocity error term J2 presents more oscillations before the three hundred generations than the term J1, where the J1 behavior is more stable from one hundred generations. This indicates that the achievement of minimizing the position error term J1 is more evident than the velocity error term J2. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Mean convergence graphs of the evolution in the performance function and its terms through generations in thirty executions of each algorithm. https://doi.org/10.1371/journal.pone.0325168.g009 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 8. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 9. Design variable vector obtained with different bioinspired techniques. https://doi.org/10.1371/journal.pone.0325168.t009 Even though in Fig 9a the evolution in PSO of the mean objective function shows the worse performance through generations, PSO obtained a variety of solutions that resulted in an average and standard deviation higher than those of the other two algorithms (see Tables 6–8). PSO also demonstrates a greater search space exploration (see standard deviation), given these sparse solutions, in contrast to the other two algorithms. Despite this, all algorithms can find solutions with comparable performance. In general, the three algorithms, convergence to a value in the case of the weighted target and the term J1 occurs relatively quickly after two hundred generations. In contrast, the term J2 manifests a more irregular behavior until the end of the generations. In order to know about the algorithm performance and to generalize the results, the inferential statistics [89] is applied to confirm the performance among algorithms. In particular, the Friedman test followed by Bonferroni post hoc multiple comparisons with a significance level is used. Fig 10 shows the multiple comparison test considering PSO as the control method, where the winner is the algorithm that ranks closest to the left (x-axis in the illustration). Overlapping gray lines show no discernible difference between the algorithms, meaning that no conclusions about the data distributions can be drawn because the p-values in those comparisons are greater than the selected significance level. The red line that does not overlap indicates that the p-value is below the selected significance level in algorithm outcomes, meaning that the results differ significantly. The red line shows that PSO is similar to GA, indicating that both are the most promising algorithms for this problem. PSO and GA also outperform the DE behavior. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. The Bonferroni multiple comparison post hoc test with a significance level of = 0.05. The control method is PSO. https://doi.org/10.1371/journal.pone.0325168.g010 On the other hand, Table 9 presents the vectors of resulting design variables p* corresponding to the best solution concerning the performance function of the thirty runs reported in Table 6 for all algorithms. Table 10 compares the obtained parameters previously reported in Table 2 with those obtained using different algorithms, evidencing a certain percentage of variation among them. In the case of the center of mass , the variation ranges from to ; the inertia ranges from to ; the mass mi varies between the interval ; and in the case of the friction , the variation interval is . The obtained parameters are remarkably similar across the three algorithms, and those algorithms are suitable for solving the problem. However, the parameters obtained by the PSO are used because it achieves the best performance in minimizing the objective function (albeit by a small margin). Download: PPT PowerPoint slide PNG larger image TIFF original image Table 10. Percentage values of variation of the parameters reported in Table 2 with those obtained using different algorithms in the RRR robot. https://doi.org/10.1371/journal.pone.0325168.t010 On the other hand, the values of the objective function terms with the obtained parameters by the MbPBO methodology using PSO and with the initial parameters in Table 2 are displayed in Table 11. It is observed that the parameters obtained by the proposal, the terms associated with the position error J1 and the velocity error J2 improve their performances in and , respectively, regarding the obtained performance by the initial parameters. So, significant differences in the objective function performance result from the parameters obtained by the proposed MbPBO methodology, which aids in better characterizing the dynamic model in the MiL simulation. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 11. Results in the performance of the objective function terms obtained by the MbPBO methodology using PSO with respect to the initial parameters in Table 2. https://doi.org/10.1371/journal.pone.0325168.t011 Comparative analysis with direct search methods. This section analyzes the performance of direct search algorithms in solving the optimization problem related to the robot’s characterization under the proposed MbPBO methodology. The primary goal is to evaluate the trade-offs involved in using other heuristic-based optimizers with respect to bio-inspired optimization. The selected direct search optimizers include Constrained Optimization BY Linear Approximations (COBYLA), Nelder Mead, and Powell algorithms [103]. Each algorithm is executed independently thirty times using randomly generated initial conditions. The stop criterion is based on the number of objective function evaluations. The objective function evaluation number is 50000, which is related to the same value used in the bio-inspired algorithms. The descriptive statistic is presented in Table 12. The table includes the average (mean), standard deviation (), and the minimum (best) and the maximum (worst) performance function values obtained from executions. Additionally, the number of times each algorithm finds the minimum values is presented in parentheses next to the corresponding result. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 12. Summary of descriptive statistic in for the optimization process of thirty executions with direct search methods. https://doi.org/10.1371/journal.pone.0325168.t012 Based on the results presented in Table 12, the most promising direct search optimizer is Nelder–Mead, as it achieves the lowest performance function value of . Powell’s method takes second place with , followed by COBYLA, which yields a . Notably, each algorithm reaches its minimum value only once out of thirty independent executions, leading to relatively high standard deviation values and indicating a lack of consistency across runs. The comparison between bio-inspired algorithms and direct search methods (Tables 6 and 12) reveals that bio-inspired techniques not only yield superior performance values but also demonstrate higher robustness and repeatability across multiple runs. Comparative analysis with indirect search methods. Similar to the previous section, this section analyzes the performance of indirect search methods (gradient-based algorithms) in solving the optimization problem related to the robot’s characterization. The main objective is to show the importance of bio-inspired optimization in the proposed MbPBO methodology. The selected gradient-based optimizers include Sequential Least Squares Quadratic Programming (SLSQP), the Limited-memory Broyden–Fletcher–Goldfarb–Shanno with Box constraints (L-BFGS-B), the Truncated Newton Conjugate-Gradient (TNC) and Trust-Region (T-R) algorithms [104]. The same experimental setup, namely, thirty independent runs with randomly initialized conditions and a stopping criterion (50000 objective function evaluations), is adopted as in the previous analysis. The average (mean), standard deviation (σ), minimum (best), and maximum (worst) performance function values obtained from the algorithm executions are presented in Table 13. The frequency at which each algorithm achieved the minimum value is shown in parentheses. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 13. Summary of descriptive statistic in for the optimization process of thirty executions with indirect search methods (gradient-based algorithms). https://doi.org/10.1371/journal.pone.0325168.t013 The most promising gradient-based optimizers are SLSQP and L-BFGS-B based on the results presented in Table 13 because both find solutions with the performance function value of . Nevertheless, SLSQP stands out by reaching this optimal value seventeen times, indicating greater consistency. The Trust-Region (T-R) takes second place but only finds two times the performance function value . The worst gradient-based optimizer is TNC because it exhibits the largest performance function value and only finds the sub-optimal solution once out of thirty executions. It is also observed that the most promising gradient-based optimizers (SLSQP and L-BFGS-B) locate the best solution obtained with the most promising bio-inspired algorithms (PSO and GA in Table 6). Nevertheless, the gradient-based algorithms are very sensitive to the initial conditions, providing a diverse set of local solutions, as described in the standard deviation of Table 13. Although certain gradient-based algorithms (such as SLSQP and L-BFGS-B) exhibit moderate consistency achieving the optimal solution in around half of thirty executions, the advantage of bio-inspired algorithms (see Tables 6 and 13) lies in their global search capabilities and robustness to problem characteristics such as non-smoothness noise, and complex constraint landscapes. While gradient-based methods may converge efficiently when the search space is well-behaved (e.g., smooth and convex), they depend highly on initial conditions. They may also become trapped in local minima in more irregular or multimodal functions. On the other hand, bio-inspired algorithms incorporate exploration mechanisms to avoid local optima and search for in a wider design space region. These advantages are observed in the lower standard deviation of Table 6, indicating greater consistency in the obtained results. Thus, even if a subset of gradient-based methods performs similarly on this specific problem, bio-inspired algorithms present a more versatile and reliable option across diverse optimization problems where gradient information, constraint inclusion, multi-objective formulations, or parallel implementation is challenging. This makes bio-inspired algorithms particularly well-suited for the MbPBO methodology, which benefits from robust global search strategies and adaptability to complex problem landscapes. With these findings, the MbPBO methodology can be extended to using hybrid techniques known as memetic algorithms. The memetic algorithm combines the global exploration capabilities of the bio-inspired algorithms with the fine search (exploitation capabilities) of gradient-based techniques. This hybrid extension represents a possible future work in the proposed methodology. (i) Trade-off selection. As mentioned in the description of the MbPBO methodology, the multi-objective problem is transformed into a single-objective problem by a weighted sum approach. A set of non-dominated solutions must be generated by establishing different weights in the optimization problem (34)-(36), forming the approximate Pareto front in the performance function space. The Pareto front helps the designer make decisions on the most suitable trade-off for the specific application. Once the terms of the objective function are normalized, new experiments are carried out to find the best trade-offs in the terms J1 and J2 (quadratic errors of position and velocity). Then, fourteen combinations of weights and in the interval are tested. The first two combinations find the solution of the Pareto front without a trade-off, i.e., by setting and . The next nine combinations consider uniform intervals from to subject to  +  . The three final weights are assigned according to the best promising weight interval found in the previous experiment. The promising weight interval is visualized once the previous experiments with different weights are carried out. Thirty independent executions of the bio-inspired algorithm are performed by establishing different weights in the optimization problem (34)-(36) and using the same DE algorithm with the same parameters given above. This experiment aims to identify the combination of weights that results in a helpful trade-off solution by including the best solution (the solution that further minimizes the performance function from the algorithm executions) found in each combination. Table 4 shows the chosen weight values, the weighted normalized objective function , and its weighted non-normalized terms J1 and J2. The values in boldface indicate the minimum ones in each column. As the uniform distribution of weights does not guarantee uniform distribution of the found Pareto solutions [98], and even the bio-inspired optimizer can stagnate in local solutions [83], the obtained solutions are sorted according to non-dominance verification. The Pareto solutions [83] of results from Table 4 are displayed in Table 5 and the graphical representation of the Pareto front with the normalized terms and is shown in Fig 8. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Normalized Pareto front obtained by using different weights in the optimization process. https://doi.org/10.1371/journal.pone.0325168.g008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Values of the objective function and its non-normalized terms , obtained by using different weights in the optimization process. https://doi.org/10.1371/journal.pone.0325168.t004 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 5. Set of Pareto solutions obtained by the results of Table 4. https://doi.org/10.1371/journal.pone.0325168.t005 Once the Pareto front is obtained, the decision maker’s preference provides one useful trade-off for the application. In this case, the most preferred solution is the Pareto solution with the minimum value in the term J1 and with a suitable value in the term J2. Based on Table 5 and Fig 8, the combination of weights , stands out for its ability to further minimize the trade-off in the performance function and the position quadratic errors, also providing an acceptable velocity quadratic error. The selected weights are employed in the next step to analyze different bio-inspired algorithms to optimize the RRR robot’s characterization problem. (ii) Algorithm selection. The choice of the DE/rand/1/bin algorithm for the experiments performed previously is because it is one of the most widely used alternatives in the reviewed literature [99] (Fig 1). However, the MbPBO methodology description in this work suggests using two more techniques to solve the optimization problem with the selected obtained weights (, ). Therefore, the next frequently used algorithms in the reviewed literature are selected due to their adaptability, population-based strategy, robust performance, and simplicity [85–87]. In particular, the Particle Swarm Optimization (PSO) [39] with the best global topology and the Genetic Algorithm (GA) in its generational variant [100] are also selected to optimize the RRR robot’s characterization problem described before (34)-(36). As mentioned in the methodology description, it is only possible to determine which algorithm presents the best results for a specific problem once it is experimentally proven [88]. Therefore, given the previous algorithms, the present work proposes a comparison based on the results obtained to identify which bio-inspired technique generally performs best when solving the characterization problem. The algorithm parameters must be tuned either “by hand” through systematic manual settings, “by analogy” through the information suggested in the state-of-the-art for similar problems, or using self-tuning mechanisms. The process “by hand” (where algorithm parameters are varied sequentially and ordered to select the best configuration that solves the problem) is also suggested in [101] for algorithms like DE. Therefore, the latter is the process chosen in this work to tune such parameters. The set of parameters presented below was obtained from this trial and error procedure, where the final tuning values for the DE/rand/1/bin algorithm are as follows: Maximum number of generations Gmax = 500, number of individuals NP = 100, scaling factor , and crossover rate CR = 0.2. In the case of PSO, the following values are obtained: Maximum number of generations Gmax = 500, number of individuals NP = 100, weight of the current best C1 = 1.8, weight of the global best C2 = 2.2 and speed factor . The parameters for the GA are set as: Maximum number of generations Gmax = 500, number of individuals NP = 100, mutation probability PM = 1/6, and crossover probability PC = 0.9. It is important to point out that the computational complexity of PSO and GA is similar to DE [102]. This is obtained by the number of calls to the objective function evaluation , the maximum iteration number Gmax, and the population size NP. Then, the computational complexity in a “big O” notation is given by NP Gmax). So, in all algorithms, the total evaluation number of the objective function is 50000 to make a fair comparison in the next step. (iii) Performance analysis of bio-inspired algorithms. Thirty independent runs were performed for each of the above algorithms. Each algorithm was programmed in the M language of MATLABTM, on a commercial PC with an IntelTM Core i7 @ 3.2 GHz processor, 16 GB of RAM. The descriptive statistics are presented below based on the criteria of the normalized weighted error . The average (mean), standard deviation , the best (minimum), and the worst (maximum) are specified based on the data obtained from the best individual of the last generation of each run. Tables 6-8 summarize the results obtained by each algorithm for the different criteria. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 6. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t006 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 7. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t007 The tables show that the PSO algorithm gets the best solution in terms of the objective function , as seen in Table 6. This is highlighted in boldface. Also, PSO presents the same behavior behavior concerning the individual objective function terms J1 and J2 in Tables 7 and 8. The solutions obtained with the DE techniques and GA for turn out to be very close to PSO with a percentage increment of 0.0094% for DE and 0.00000148% for GA. Additionally, the plots in Fig 9a shows a follow-up of the evolution through generations of the mean objective function of the best solutions in thirty executions (convergence graph of the mean objective function) using different optimization techniques. It is observed that while PSO converges the fastest, its overall performance through the algorithm executions is not as good on average. For completeness, Fig 9b and 9c also shows the convergence graph of the mean objective function terms using different optimization techniques, i.e., the evolution of the position error term J1 and the velocity error term J2 through generations are displayed. It is observed that the velocity error term J2 presents more oscillations before the three hundred generations than the term J1, where the J1 behavior is more stable from one hundred generations. This indicates that the achievement of minimizing the position error term J1 is more evident than the velocity error term J2. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Mean convergence graphs of the evolution in the performance function and its terms through generations in thirty executions of each algorithm. https://doi.org/10.1371/journal.pone.0325168.g009 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 8. Summary of descriptive statistics in for the optimization process of thirty executions with different bioinspired algorithms. https://doi.org/10.1371/journal.pone.0325168.t008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 9. Design variable vector obtained with different bioinspired techniques. https://doi.org/10.1371/journal.pone.0325168.t009 Even though in Fig 9a the evolution in PSO of the mean objective function shows the worse performance through generations, PSO obtained a variety of solutions that resulted in an average and standard deviation higher than those of the other two algorithms (see Tables 6–8). PSO also demonstrates a greater search space exploration (see standard deviation), given these sparse solutions, in contrast to the other two algorithms. Despite this, all algorithms can find solutions with comparable performance. In general, the three algorithms, convergence to a value in the case of the weighted target and the term J1 occurs relatively quickly after two hundred generations. In contrast, the term J2 manifests a more irregular behavior until the end of the generations. In order to know about the algorithm performance and to generalize the results, the inferential statistics [89] is applied to confirm the performance among algorithms. In particular, the Friedman test followed by Bonferroni post hoc multiple comparisons with a significance level is used. Fig 10 shows the multiple comparison test considering PSO as the control method, where the winner is the algorithm that ranks closest to the left (x-axis in the illustration). Overlapping gray lines show no discernible difference between the algorithms, meaning that no conclusions about the data distributions can be drawn because the p-values in those comparisons are greater than the selected significance level. The red line that does not overlap indicates that the p-value is below the selected significance level in algorithm outcomes, meaning that the results differ significantly. The red line shows that PSO is similar to GA, indicating that both are the most promising algorithms for this problem. PSO and GA also outperform the DE behavior. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. The Bonferroni multiple comparison post hoc test with a significance level of = 0.05. The control method is PSO. https://doi.org/10.1371/journal.pone.0325168.g010 On the other hand, Table 9 presents the vectors of resulting design variables p* corresponding to the best solution concerning the performance function of the thirty runs reported in Table 6 for all algorithms. Table 10 compares the obtained parameters previously reported in Table 2 with those obtained using different algorithms, evidencing a certain percentage of variation among them. In the case of the center of mass , the variation ranges from to ; the inertia ranges from to ; the mass mi varies between the interval ; and in the case of the friction , the variation interval is . The obtained parameters are remarkably similar across the three algorithms, and those algorithms are suitable for solving the problem. However, the parameters obtained by the PSO are used because it achieves the best performance in minimizing the objective function (albeit by a small margin). Download: PPT PowerPoint slide PNG larger image TIFF original image Table 10. Percentage values of variation of the parameters reported in Table 2 with those obtained using different algorithms in the RRR robot. https://doi.org/10.1371/journal.pone.0325168.t010 On the other hand, the values of the objective function terms with the obtained parameters by the MbPBO methodology using PSO and with the initial parameters in Table 2 are displayed in Table 11. It is observed that the parameters obtained by the proposal, the terms associated with the position error J1 and the velocity error J2 improve their performances in and , respectively, regarding the obtained performance by the initial parameters. So, significant differences in the objective function performance result from the parameters obtained by the proposed MbPBO methodology, which aids in better characterizing the dynamic model in the MiL simulation. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 11. Results in the performance of the objective function terms obtained by the MbPBO methodology using PSO with respect to the initial parameters in Table 2. https://doi.org/10.1371/journal.pone.0325168.t011 Comparative analysis with direct search methods. This section analyzes the performance of direct search algorithms in solving the optimization problem related to the robot’s characterization under the proposed MbPBO methodology. The primary goal is to evaluate the trade-offs involved in using other heuristic-based optimizers with respect to bio-inspired optimization. The selected direct search optimizers include Constrained Optimization BY Linear Approximations (COBYLA), Nelder Mead, and Powell algorithms [103]. Each algorithm is executed independently thirty times using randomly generated initial conditions. The stop criterion is based on the number of objective function evaluations. The objective function evaluation number is 50000, which is related to the same value used in the bio-inspired algorithms. The descriptive statistic is presented in Table 12. The table includes the average (mean), standard deviation (), and the minimum (best) and the maximum (worst) performance function values obtained from executions. Additionally, the number of times each algorithm finds the minimum values is presented in parentheses next to the corresponding result. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 12. Summary of descriptive statistic in for the optimization process of thirty executions with direct search methods. https://doi.org/10.1371/journal.pone.0325168.t012 Based on the results presented in Table 12, the most promising direct search optimizer is Nelder–Mead, as it achieves the lowest performance function value of . Powell’s method takes second place with , followed by COBYLA, which yields a . Notably, each algorithm reaches its minimum value only once out of thirty independent executions, leading to relatively high standard deviation values and indicating a lack of consistency across runs. The comparison between bio-inspired algorithms and direct search methods (Tables 6 and 12) reveals that bio-inspired techniques not only yield superior performance values but also demonstrate higher robustness and repeatability across multiple runs. Comparative analysis with indirect search methods. Similar to the previous section, this section analyzes the performance of indirect search methods (gradient-based algorithms) in solving the optimization problem related to the robot’s characterization. The main objective is to show the importance of bio-inspired optimization in the proposed MbPBO methodology. The selected gradient-based optimizers include Sequential Least Squares Quadratic Programming (SLSQP), the Limited-memory Broyden–Fletcher–Goldfarb–Shanno with Box constraints (L-BFGS-B), the Truncated Newton Conjugate-Gradient (TNC) and Trust-Region (T-R) algorithms [104]. The same experimental setup, namely, thirty independent runs with randomly initialized conditions and a stopping criterion (50000 objective function evaluations), is adopted as in the previous analysis. The average (mean), standard deviation (σ), minimum (best), and maximum (worst) performance function values obtained from the algorithm executions are presented in Table 13. The frequency at which each algorithm achieved the minimum value is shown in parentheses. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 13. Summary of descriptive statistic in for the optimization process of thirty executions with indirect search methods (gradient-based algorithms). https://doi.org/10.1371/journal.pone.0325168.t013 The most promising gradient-based optimizers are SLSQP and L-BFGS-B based on the results presented in Table 13 because both find solutions with the performance function value of . Nevertheless, SLSQP stands out by reaching this optimal value seventeen times, indicating greater consistency. The Trust-Region (T-R) takes second place but only finds two times the performance function value . The worst gradient-based optimizer is TNC because it exhibits the largest performance function value and only finds the sub-optimal solution once out of thirty executions. It is also observed that the most promising gradient-based optimizers (SLSQP and L-BFGS-B) locate the best solution obtained with the most promising bio-inspired algorithms (PSO and GA in Table 6). Nevertheless, the gradient-based algorithms are very sensitive to the initial conditions, providing a diverse set of local solutions, as described in the standard deviation of Table 13. Although certain gradient-based algorithms (such as SLSQP and L-BFGS-B) exhibit moderate consistency achieving the optimal solution in around half of thirty executions, the advantage of bio-inspired algorithms (see Tables 6 and 13) lies in their global search capabilities and robustness to problem characteristics such as non-smoothness noise, and complex constraint landscapes. While gradient-based methods may converge efficiently when the search space is well-behaved (e.g., smooth and convex), they depend highly on initial conditions. They may also become trapped in local minima in more irregular or multimodal functions. On the other hand, bio-inspired algorithms incorporate exploration mechanisms to avoid local optima and search for in a wider design space region. These advantages are observed in the lower standard deviation of Table 6, indicating greater consistency in the obtained results. Thus, even if a subset of gradient-based methods performs similarly on this specific problem, bio-inspired algorithms present a more versatile and reliable option across diverse optimization problems where gradient information, constraint inclusion, multi-objective formulations, or parallel implementation is challenging. This makes bio-inspired algorithms particularly well-suited for the MbPBO methodology, which benefits from robust global search strategies and adaptability to complex problem landscapes. With these findings, the MbPBO methodology can be extended to using hybrid techniques known as memetic algorithms. The memetic algorithm combines the global exploration capabilities of the bio-inspired algorithms with the fine search (exploitation capabilities) of gradient-based techniques. This hybrid extension represents a possible future work in the proposed methodology. Stage V: Validation of the obtained RRR robot model and its use in the MiL simulation Once the comparison described in the step iii) of Stage IV has been carried out, and the design variable vector resulting in the minimum value of the weighted objective function (p* obtained with PSO reported in Table 6) has been obtained, the last stage of the proposed MbPBO methodology is the comparative evaluation concerning the similarity of the system model’s states with those of the plant. So, the first step analyzes the model obtained in open-loop control with the reference signal from Stage II. In addition, this stage also focuses on assessing the parameterized model’s ability to accurately replicate the behavior of the testbed platform in a wide range of applications to test different MiL simulation environments. This evaluation essentially tests the model’s predictive capability with different control strategies and monitoring scenarios by comparing position and velocity states between the model and the actual testbed platform during specific tasks. So, the second step of this stage analyses the model obtained in closed-loop for the regulation and tracking tasks. In this stage, graphic comparisons and the Pearson correlation coefficient analysis between experimental and simulation signals are performed to assess the similarity or variability between the system model’s states and the plant’s numerically. The steps of this stage are described below. (i) Model validation in open-loop control. Figs 11 and 12 show the Cartesian angular position and angular velocity behavior of the RRR robot obtained from the approximated model (signal displayed in solid blue line) and the real experiment (signal represented by the dotted magenta line). The Cartesian position in the real experiment and the simulated plant is represented by and , respectively. Similarly, the angular positions are defined by q1, q2, q3 for the real experiment and , , for the behavior of the obtained model. The angular velocities are set as , , for the real experiment and , , for the simulated plant. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Robot’s end-effector behavior in the Cartesian plane with the solution obtained by PSO in simulation, against the real states of the experiment when system states are excited by the open-loop control signals. https://doi.org/10.1371/journal.pone.0325168.g011 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Angular position and velocity behaviors of the RRR robot when the system states are excited by the open-loop control signals. Real experiment (q1, q2, q3, , , ) vs. approximated model (, , , , , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g012 For a better numerical appreciation of the results observed in Figs 11 and 12, Table 14 presents the correlation coefficient between the real measurements and the results of the approximated model for each state previously examined. The coefficient is a statistical measure that indicates the strength and direction of the linear relationship between two variables (Pearson correlation). So, this measure would indicate the similarity (linear relationship) between the data obtained by the model and the experimental data. A coefficient of -1 would indicate a perfect inverse correlation, i.e., with the same input data, the model provides an output with the same magnitude but in an opposite direction than the output of the real system. The value of 0 indicates no correlation, which means that the model data and the experimental data do not exhibit a linear relationship. Hence, the model could not be helpful to describe the experimental (real) behavior. The value of 1 indicates a perfect positive correlation (strong correlation) between the variables, i.e., the model data matches perfectly to the experimental data with the same input data, meaning that the model is useful for describing the real behavior. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 14. Correlation coefficient between the robot signals and the approximated model behavior when the system states are excited by the open-loop control signals. https://doi.org/10.1371/journal.pone.0325168.t014 The results show suitable relationships between the actual behavior and the approximated model with a average correlation of 0.9015 in the angular space and 0.8969 in the Cartesian space. In particular, the average correlations for the angular position and angular velocity are 0.9656 and 0.8375, respectively. In particular, the velocity states exhibit lower correlation coefficients than the position states concerning the actual experimental behavior, attributed to the noisy measurement in the velocity states. The noise in the identification tasks is a significant challenge, often necessitating well-tuned filtering techniques, advanced measurement instrumentation, or meticulous trajectory planning [77]. (ii) Model validation in closed-loop. In this step, two control strategies are also employed to test two different MiL simulation environments: 1) A PID control strategy for regulation at a single point in the joint space and 2) A PD+G control strategy for following a smooth trajectory in the joint space. The comparisons between the model and the testbed platform, which include graphical representations and corresponding correlation coefficients, form a thorough validation process. This process instills confidence in the model’s accuracy by evaluating how well the model’s predicted states match the actual measured states. The comprehensive analysis aims to validate the model’s effectiveness in capturing the dynamics and behaviors of the testbed platform under different operational conditions. This study’s primary validation measures are the comparisons between measured and estimated position and velocity states. However, additional validation measures, as suggested by [93], could include the prediction of the necessary torque in the actuators to achieve a specific task or how far the original model parameters (perhaps required from the manufacturer) lie against the predicted values (confidence interval). These validation measures typically require more complex and potentially expensive means of data collection and analysis. Therefore, while these measures are valuable for a comprehensive validation of the model’s predictive capability, they are not included in the scope of this particular study. The focus here remains on validating the positional and velocity states as a primary indicator of the model’s accuracy in representing the behavior of the testbed platform. Regulate a point in the joint space with the PID control strategy Table 15 presents the positions to be reached in the joint space for each degree of freedom and the PID controller parameters used for each joint, applicable to both the model and the laboratory testbed platform. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 15. PID controller parameters and the desired position in the joint space used in each d.o.f. for the validation of the obtained RRR robot model in the regulation task. https://doi.org/10.1371/journal.pone.0325168.t015 Fig 13 illustrates the workspace of the manipulator , where the thick black lines represent the robot’s links with the end effector at one end, the solid blue line shows the measured real values , and the thick magenta dotted line represents the predicted values . Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 13. Robot’s end-effector behavior in the workspace of the obtained model against the real experiment for the regulation task using PID controller. Black continuous lines indicate the robot links. https://doi.org/10.1371/journal.pone.0325168.g013 Concerning the controller performances, Fig 14a, 14c, 14e presents the actual and approximated position states for the regulation task. Likewise, Fig 14b, 14d, 14f includes the real and approximate velocity states for the same task. The control signals applied to the model and testbed platform are presented in Fig 15. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 14. Angular position and velocity behaviors of the RRR robot in the regulation task using PID controller. Real experiment (q1, q2, q3, , , ) vs. approximated model (, , , , , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g014 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 15. Control signals () for each d.o.f. obtained from the regulation PID control experiment. Real experiment (u1, u2, u3) vs. approximated model (, , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g015 Table 16 presents the respective correlation coefficients between the states of the testbed platform and those of the model for this experiment. Based on this table and the corresponding Fig 14 in the regulation task with the PID controller, it is observed that the overall correlations in the Cartesian space of the robot are very confident in the interval , meanwhile the angular position space presents a correlation in the interval and the angular velocity space provides a correlation in the interval . With this information, the approximated model successfully behaves as the robot’s end-effector behavior of the testbed platform. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 16. Correlation coefficient between the obtained model signals and the robot signals in the regulation task with the PID controller (closed loop). https://doi.org/10.1371/journal.pone.0325168.t016 Smooth tracking of a path with the PD+G control strategy In this test, the trajectory consists of tracking three points in the joint space, executed in two different quadrants of the Cartesian space. Those points are joined by a third-order Bezier polynomial, resulting in a smooth trajectory. Table 17 displayed the points reached by this trajectory and the controller parameters used for the PD+G control in both the modeling and experimental setups. Fig 16 presents a view of the manipulator’s workspace (Xw,YW), showing the comparison of both the real measured behavior in solid lines and the predicted one in dotted lines, following the convention of previous plots, but this time for the tracking task. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 16. Robot’s behavior in the workspace of the obtained model against the real experiment for the tracking task using PD+G controller. Black continuous lines indicate the robot links. https://doi.org/10.1371/journal.pone.0325168.g016 The plots in Fig 17a, 17c, 17e present the real and approximate position states of for the tracking task. Similarly, the real and approximate velocity states are depicted in Fig 17b, 17d, 17f for the same trajectory tracking task. The control signals of the model and the testbed platform can also be observed in Fig 18. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 17. Angular position and velocity behaviors of the robot in the tracking task using PD+G controller. Real experiment (, , , , , ) vs. approximated model (, , , , , ) with the solution obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g017 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 18. Control signals () for each d.o.f. obtained from the tracking PD+G control experiment. Real experiment (u1, u2, u3) vs. approximated model (, , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g018 Table 18 presents the respective correlation coefficients between the states of the testbed platform and those of the model for this experiment. Based on this table and Fig 17 in the tracking task with the PD+G controller, it is observed that the overall correlations in the Cartesian space of the robot and the angular position space are very confident with a value around one. Meanwhile, the angular velocity space provides a correlation in the interval . With this information, the approximated model successfully behaves as the robot’s end-effector behavior of the testbed platform. Highlights of the validation results. The model validation in the open loop control indicates a high average correlation among all joint states of 0.9015, indicating a high representation of the real system with the obtained parameterized model. In particular, the average correlation in the angular position, angular velocity, and Cartesian position is 0.9656, 0.8375, and 0.8969, respectively. On the other hand, the tests in the MiL simulation using the obtained parameterized model for the regulation and the tracking control problems (closed loop control) shows a high average correlation among all joint states in both control problems of 0.9022, revealing the model’s effectiveness for general-purpose use in MiL tests. In particular, the correlations given for the Cartesian space is in the interval and for the angular position space is in the interval providing an average correlation of 0.9906 and 0.9878, respectively. Meanwhile, the correlation in the angular velocity space is in the interval with an average correlation of 0.8165. The slight decrement in the angular velocity correlation is attributed to the estimation of the angular velocity by using finite difference because it produces noise in the velocity measure and also in the trade-off selected in Stage IV. (i) Model validation in open-loop control. Figs 11 and 12 show the Cartesian angular position and angular velocity behavior of the RRR robot obtained from the approximated model (signal displayed in solid blue line) and the real experiment (signal represented by the dotted magenta line). The Cartesian position in the real experiment and the simulated plant is represented by and , respectively. Similarly, the angular positions are defined by q1, q2, q3 for the real experiment and , , for the behavior of the obtained model. The angular velocities are set as , , for the real experiment and , , for the simulated plant. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Robot’s end-effector behavior in the Cartesian plane with the solution obtained by PSO in simulation, against the real states of the experiment when system states are excited by the open-loop control signals. https://doi.org/10.1371/journal.pone.0325168.g011 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Angular position and velocity behaviors of the RRR robot when the system states are excited by the open-loop control signals. Real experiment (q1, q2, q3, , , ) vs. approximated model (, , , , , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g012 For a better numerical appreciation of the results observed in Figs 11 and 12, Table 14 presents the correlation coefficient between the real measurements and the results of the approximated model for each state previously examined. The coefficient is a statistical measure that indicates the strength and direction of the linear relationship between two variables (Pearson correlation). So, this measure would indicate the similarity (linear relationship) between the data obtained by the model and the experimental data. A coefficient of -1 would indicate a perfect inverse correlation, i.e., with the same input data, the model provides an output with the same magnitude but in an opposite direction than the output of the real system. The value of 0 indicates no correlation, which means that the model data and the experimental data do not exhibit a linear relationship. Hence, the model could not be helpful to describe the experimental (real) behavior. The value of 1 indicates a perfect positive correlation (strong correlation) between the variables, i.e., the model data matches perfectly to the experimental data with the same input data, meaning that the model is useful for describing the real behavior. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 14. Correlation coefficient between the robot signals and the approximated model behavior when the system states are excited by the open-loop control signals. https://doi.org/10.1371/journal.pone.0325168.t014 The results show suitable relationships between the actual behavior and the approximated model with a average correlation of 0.9015 in the angular space and 0.8969 in the Cartesian space. In particular, the average correlations for the angular position and angular velocity are 0.9656 and 0.8375, respectively. In particular, the velocity states exhibit lower correlation coefficients than the position states concerning the actual experimental behavior, attributed to the noisy measurement in the velocity states. The noise in the identification tasks is a significant challenge, often necessitating well-tuned filtering techniques, advanced measurement instrumentation, or meticulous trajectory planning [77]. (ii) Model validation in closed-loop. In this step, two control strategies are also employed to test two different MiL simulation environments: 1) A PID control strategy for regulation at a single point in the joint space and 2) A PD+G control strategy for following a smooth trajectory in the joint space. The comparisons between the model and the testbed platform, which include graphical representations and corresponding correlation coefficients, form a thorough validation process. This process instills confidence in the model’s accuracy by evaluating how well the model’s predicted states match the actual measured states. The comprehensive analysis aims to validate the model’s effectiveness in capturing the dynamics and behaviors of the testbed platform under different operational conditions. This study’s primary validation measures are the comparisons between measured and estimated position and velocity states. However, additional validation measures, as suggested by [93], could include the prediction of the necessary torque in the actuators to achieve a specific task or how far the original model parameters (perhaps required from the manufacturer) lie against the predicted values (confidence interval). These validation measures typically require more complex and potentially expensive means of data collection and analysis. Therefore, while these measures are valuable for a comprehensive validation of the model’s predictive capability, they are not included in the scope of this particular study. The focus here remains on validating the positional and velocity states as a primary indicator of the model’s accuracy in representing the behavior of the testbed platform. Regulate a point in the joint space with the PID control strategy Table 15 presents the positions to be reached in the joint space for each degree of freedom and the PID controller parameters used for each joint, applicable to both the model and the laboratory testbed platform. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 15. PID controller parameters and the desired position in the joint space used in each d.o.f. for the validation of the obtained RRR robot model in the regulation task. https://doi.org/10.1371/journal.pone.0325168.t015 Fig 13 illustrates the workspace of the manipulator , where the thick black lines represent the robot’s links with the end effector at one end, the solid blue line shows the measured real values , and the thick magenta dotted line represents the predicted values . Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 13. Robot’s end-effector behavior in the workspace of the obtained model against the real experiment for the regulation task using PID controller. Black continuous lines indicate the robot links. https://doi.org/10.1371/journal.pone.0325168.g013 Concerning the controller performances, Fig 14a, 14c, 14e presents the actual and approximated position states for the regulation task. Likewise, Fig 14b, 14d, 14f includes the real and approximate velocity states for the same task. The control signals applied to the model and testbed platform are presented in Fig 15. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 14. Angular position and velocity behaviors of the RRR robot in the regulation task using PID controller. Real experiment (q1, q2, q3, , , ) vs. approximated model (, , , , , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g014 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 15. Control signals () for each d.o.f. obtained from the regulation PID control experiment. Real experiment (u1, u2, u3) vs. approximated model (, , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g015 Table 16 presents the respective correlation coefficients between the states of the testbed platform and those of the model for this experiment. Based on this table and the corresponding Fig 14 in the regulation task with the PID controller, it is observed that the overall correlations in the Cartesian space of the robot are very confident in the interval , meanwhile the angular position space presents a correlation in the interval and the angular velocity space provides a correlation in the interval . With this information, the approximated model successfully behaves as the robot’s end-effector behavior of the testbed platform. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 16. Correlation coefficient between the obtained model signals and the robot signals in the regulation task with the PID controller (closed loop). https://doi.org/10.1371/journal.pone.0325168.t016 Smooth tracking of a path with the PD+G control strategy In this test, the trajectory consists of tracking three points in the joint space, executed in two different quadrants of the Cartesian space. Those points are joined by a third-order Bezier polynomial, resulting in a smooth trajectory. Table 17 displayed the points reached by this trajectory and the controller parameters used for the PD+G control in both the modeling and experimental setups. Fig 16 presents a view of the manipulator’s workspace (Xw,YW), showing the comparison of both the real measured behavior in solid lines and the predicted one in dotted lines, following the convention of previous plots, but this time for the tracking task. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 16. Robot’s behavior in the workspace of the obtained model against the real experiment for the tracking task using PD+G controller. Black continuous lines indicate the robot links. https://doi.org/10.1371/journal.pone.0325168.g016 The plots in Fig 17a, 17c, 17e present the real and approximate position states of for the tracking task. Similarly, the real and approximate velocity states are depicted in Fig 17b, 17d, 17f for the same trajectory tracking task. The control signals of the model and the testbed platform can also be observed in Fig 18. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 17. Angular position and velocity behaviors of the robot in the tracking task using PD+G controller. Real experiment (, , , , , ) vs. approximated model (, , , , , ) with the solution obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g017 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 18. Control signals () for each d.o.f. obtained from the tracking PD+G control experiment. Real experiment (u1, u2, u3) vs. approximated model (, , ) with the solution p* obtained by PSO. https://doi.org/10.1371/journal.pone.0325168.g018 Table 18 presents the respective correlation coefficients between the states of the testbed platform and those of the model for this experiment. Based on this table and Fig 17 in the tracking task with the PD+G controller, it is observed that the overall correlations in the Cartesian space of the robot and the angular position space are very confident with a value around one. Meanwhile, the angular velocity space provides a correlation in the interval . With this information, the approximated model successfully behaves as the robot’s end-effector behavior of the testbed platform. Highlights of the validation results. The model validation in the open loop control indicates a high average correlation among all joint states of 0.9015, indicating a high representation of the real system with the obtained parameterized model. In particular, the average correlation in the angular position, angular velocity, and Cartesian position is 0.9656, 0.8375, and 0.8969, respectively. On the other hand, the tests in the MiL simulation using the obtained parameterized model for the regulation and the tracking control problems (closed loop control) shows a high average correlation among all joint states in both control problems of 0.9022, revealing the model’s effectiveness for general-purpose use in MiL tests. In particular, the correlations given for the Cartesian space is in the interval and for the angular position space is in the interval providing an average correlation of 0.9906 and 0.9878, respectively. Meanwhile, the correlation in the angular velocity space is in the interval with an average correlation of 0.8165. The slight decrement in the angular velocity correlation is attributed to the estimation of the angular velocity by using finite difference because it produces noise in the velocity measure and also in the trade-off selected in Stage IV. Application of the MbPBO methodology to the RFSEA The main purpose of this section is to provide useful information about the applicability and the flexibility of the MbPBO methodology in different electro-mechanical systems. In this case, the proposed MbPBO methodology is applied to the model parameterization of the Reaction Force-sensing Series Elastic Actuator (RFSEA), confirming the similarities in the behavior of the obtained model with respect to the real system’s behavior. In what follows, some aspects of the methodology are highlighted. The RFSEA is a series elastic linear actuator [105] where the tip of actuation (the ball screw) involves an active and passive movement through the displacement of the pulley-belt transmission and spring, respectively. The Series Elastic Actuator (SEA) has been widely used in robotic leg [106], biomechanical devices [107–109], mechatronic systems [110], lower and upper limb rehabilitation devices [111,112], among others. The most important aspects of the MbPBO methodology are highlighted next. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 17. PD+G controller parameters and the desired positions in the joint space joined by third-order Bezier polynomials for the validation of the obtained RRR robot model in the tracking task. https://doi.org/10.1371/journal.pone.0325168.t017 Fig 19 represents the RFSEA system divided into the three most important parts of movements in different colors: The motor’s belt-pulley transmission, the spring movement, and the linear motion of the ball screw. The variables associated in the schematic diagram of the RFSEA are: (A) For the belt-pulley transmission: the motor’s angular position , the input torque (motor’s torque), the output torque , the motor’s current ia, the motor’s inertia Jm, the motor’s damping coefficient Bm and the radii of the drive pulleys R1 and R2. (B) For the spring movement, the spring’s linear displacement xs, the spring’s stiffness coefficient Ks, the spring’s mass Ms, the spring’s damping coefficient Bs, the pulley radius re for the linear spring’s displacement and its angular displacement θe. (C) For the ball screw’s linear motion, the ball screw displacement xl, the ball screw’s damping coefficient Bl, the ball screw’s mass Ml, the screw pitch l and the load’s reaction force Dl. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 19. RFSEA system. https://doi.org/10.1371/journal.pone.0325168.g019 The dynamic equations of the RFSEA consist of the mechanical (37)-(39) and electrical (40) dynamics [113–115], where is the transmission ratio of the ball screw with as the belt-pulley transmission ratio obtained when the larger pulley with radius R1 is the output and the smaller pulley with radius R2 is the input. (37)(38)(39)(40) The dynamic coupling of the mechanical and electrical parts of the RFSEA assuming can be expressed in state space with the input signal and the reaction force Dl. The state-space representation of the RFSEA dynamic model results in (41), where , ,  +  ,  +   +  , − , ,  +  , − ,  +   +  , , , , , and . (41) where: Download: PPT PowerPoint slide PNG larger image TIFF original image Table 18. Correlation coefficient between the obtained model signals and the robot signals in the tracking task with the PD+G controller (closed loop). https://doi.org/10.1371/journal.pone.0325168.t018 The Stage I is finished with the obtained model of the RFSEA. The characterization of the system’s behavior in Stage II is carried out in the same manner as described in the MbPBO methodology. The specific set of parameters to be found and the objective function are selected in the optimization process for Stage III. In this case, the selected design parameters vector is related to the damping coefficients with the terms of the objective functions as , where is the estimated state. In addition, the normalization is done with the corresponding maximum values of the objective function terms obtained in the procedure (, , , 10−5, ). The values of the rest model parameters associated to the RFSEA are described through the datasheet of components, or by directly using measuring instruments. Following Stage IV of the proposed model parameterization approach, the selected weights are , giving the same importance to all terms in the objective function. In Table 19, the parameters obtained using the MbPBO optimization process are indicated with an asterisk, along with the remaining parameters employed. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 19. Parameters of the RFSEA. https://doi.org/10.1371/journal.pone.0325168.t019 With the obtained parameters shown in Table 19, the parameterized model is completed, and the validation of the obtained parameterized model using an open-loop controller is carried out as described in Stage V. In this case, the drivers to control the RFSEA motor are configured in voltage mode, so the reference signals in the open-loop controller are related to the voltage and as a consequence the open-loop control signal is u = ref (a cascade controller is not applied in the RFSEA case). The behaviors of the state vectors for the real experiment and the approximated model with the same persistent excitation control signal are given in Fig 20. The bar above the variable indicates the behavior of the obtained model. The Pearson correlation coefficients of the signals related to the screw displacement, the screw velocity, and the current present a suitable correlation of and 0.9226, respectively. Those correlations indicate that the obtained model represents a similar behavior to the real system. The worst correlation coefficients are given by the spring position and spring velocity ( and , respectively). This is attributed because those measures present small magnitudes (see Fig 20b, 20d) due to spring states are not excited in the experiments by applying the external force Dl and are mainly affected by the vibrations of the testbed platform. So, those spring measures are not representative in the study. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 20. State behaviors of the obtained RFSEA model against the real experiment signals when the system states are excited by the open-loop control signals. Real experiment (xl, xs, , , ia) vs. approximated model (, , , , ). https://doi.org/10.1371/journal.pone.0325168.g020 On the other hand, the validation of the obtained parameterized model using a closed-loop controller for testing two different MiL simulation environments is carried out as described in Stage V. Table 20 provides the respective correlation coefficients between the states of the RFSEA testbed platform and those of the model for the regulation and tracking tasks. The representative figures for those signals are presented in Fig 21. Based on this table, it is observed that the average correlation in the regulation and tracking tasks in the ball screw’s linear displacement is 0.9743 and in the ball screw’s linear velocity is 0.8356, indicating high correlation values, i.e., the model data are highly related to the experimental data using the same input data, meaning that the model is useful for describing the real behavior of those signals. In the motor’s current, the average correlation is 0.4263, indicating a moderate correlation value and that the obtained model behavior could present deviations from the real behavior, which is mainly attributed to a low current transducer sensitivity. In the case of the spring position and its velocity, these signals present an average correlation of −0.1749 and 0.0446, respectively. Those correlation values indicate that there is no relationship between the model data and the real behavior. Nevertheless, those weak correlations are attributed to the spring signals not being excited in either the model or the experimentation. So, their movements in the experimental data have small magnitudes (see Fig 21c, 21d, 21g, 21h) mainly given by the vibrations of the testbed platform which are not considered in the model. For instance, the movement of the spring is not perturbed by an external load in both experiments. Then, the spring position and velocity correlations are not representative in the study. So, the results obtained in the MiL simulation environments indicate that the obtained RFSEA model can describe, to a large extent, the real behavior. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 21. Behavior of the RFSEA’s states in the regulation (left) and tracking (right) tasks. Real experiment (xl, xs, , , ia) vs. approximated model (, , , , ). https://doi.org/10.1371/journal.pone.0325168.g021 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 20. Correlation coefficient ρ between the obtained model signals and the system signals in the regulation and tracking tasks with the PID controller in the RFSEA. https://doi.org/10.1371/journal.pone.0325168.t020 General discussion Based on the outcomes achieved at various stages of applying the MbPBO methodology, the following is highlighted: Bio-inspired techniques are widely regarded as effective alternatives for modeling and identification tasks. Based on a comparative analysis involving PSO, GA, and DE algorithms in solving optimization problems, the PSO algorithm consistently converges toward the best solution across multiple executions. A persistent pattern is observed in the percentage variation between the nominal values of the experiment and the solutions obtained using the three algorithms. These findings underscore the importance of the MbPBO methodology in accurately parametrizing models to achieve behavior closer to real systems. Despite unmodeled uncertainties in the investigation, such as noise signal, backlash between gears, and hysteresis, among others, the simulated states exhibit a high average general correlation of 0.9019 considering open-loop control and various tasks and control strategies, including regulation and tracking tasks (MiL tests) in the joint space. Particularly, this model achieves an average correlation of 0.9015 and 0.9022 in open and closed-loop control, respectively. This high positive correlation indicates that the behavior of the obtained model is strongly related to the behavior of the real robotic system. This means that the obtained model can be used to predict the system’s behavior. The MbPBO methodology applied to a different electro-mechanical system (RFSEA) reveals the flexibility of the proposal. In this case, considering the states that were excited in the RFSEA, the correlations in open and closed loop control are 0.9305 and 0.7454, providing a high average general correlation of 0.8379. This correlation indicates a suitable approximated model. Obtaining a representative system model that closely approximates reality offers several advantages. Firstly, it eliminates the need for repetitive tests on the actual system when maximizing task performance, such as in controller tuning. Secondly, the methodology includes model validation within the MiL approach, accelerating the implementation of model-based design (MiL-SiL-HiL). Conclusions This paper proposes the Model-Based Parametric Bioinspired Optimization (MbPBO) methodology. The proposal utilizes a nominal model based on mathematical relationships describing the system and identifies specific parameters through a bio-inspired optimization process framed as an identification experiment. The MbPBO methodology aims to develop an effective model for predicting the behavior of robotic systems integrated within the MiL simulation framework. The MbPBO methodology can also be applied to different systems, revealing a suitable obtained model when applied to the RFSEA electro-mechanical system. The application of MbPBO methodology to a three-degree-of-freedom manipulator robot case study demonstrates that the methodology produces a parametrized model with a high average general correlation of 0.9019 considering open-loop control and closed-loop control strategies (MiL tests). When the MbPBO methodology is applied to the RFSEA, a high average general correlation of 0.8379 is obtained. Those results validate the parametrized model’s effectiveness and methodology, showing the model’s capacity to accurately forecast the behavior of the system and its strong correlation with the real system. The parametrized model obtained through MbPBO methodology enhances the efficiency of the model-based design cycle (MiL-SiL-HiL). Additionally, it facilitates iterative adjustments or learning processes to optimize the performance of robotic systems without requiring direct interaction with the physical system. Future research consists of developing virtual laboratories for robotic systems, where the MiL test will be applied to validate the controller design and the optimal tuning of its parameters. Virtual laboratories provide a controlled, cost-effective environment for testing and optimizing robotic control strategies through high-fidelity simulations. They reduce risks and expenses by allowing researchers to refine controller parameters and detect failures before real-world deployment. Additionally, these labs enable remote collaboration, scalability, and experimentation across different scenarios, fostering advancements in adaptive control, human-robot interaction, and reinforcement learning-based optimizations to enhance robotic system efficiency and reliability. Acknowledgments Authors acknowledge the support from the Comisión de Operación y Fomento de Actividades Académicas of the Instituto Politécnico Nacional (COFAA-IPN), the Secretaría de Investigación y Posgrado of the IPN (SIP-IPN), the Expert’s Network in Robotics and Mechatronics of the IPN, the Colegio de Ciencia y Tecnología (CCyT) of the Universidad Autónoma de la Ciudad de México (UACM) and the Secretaría de Ciencia, Humanidades, Tecnología e Innovación (SECIHTI). The first author acknowledges the support from the CONAHCYT in Mexico through a scholarship to pursue graduate studies at CIDETEC-IPN.
TI - Model parameterization of robotic systems through the bio-inspired optimization
JF - PLoS ONE
DO - 10.1371/journal.pone.0325168
DA - 2025-06-04
UR - https://www.deepdyve.com/lp/public-library-of-science-plos-journal/model-parameterization-of-robotic-systems-through-the-bio-inspired-REv2d8nCp7
SP - e0325168
VL - 20
IS - 6
DP - DeepDyve
ER -