EMG-Based Biomechanical Cybernetic Control of an Avatar Robotic Fish

1. Introduction

Avatar robotics involves remotely controlling a robot to interact with the physical environment on behalf of a human operator, enabling them to virtually embody the robot and perform actions as if physically present. This transformative technology extends human presence to remote or hazardous locations, with applications spanning space exploration, disaster response, remote inspection, telemedicine, and diverse domains. Leveraging progress in robotics, teleoperation systems, sensory feedback interfaces, and communication networks, avatar robotics enhances human capabilities, ensures safer operations, and broadens human presence and expertise in various fields.

Furthermore, cybernetic control functions as a regulatory system utilizing feedback mechanisms to uphold stability and achieve desired outcomes. The incorporation of feedback loops is central to cybernetic control systems, continuously monitoring a system’s behavior, comparing it to a reference state, and generating corrective actions to address any deviations. This iterative feedback process facilitates self-regulation and goal attainment within the system. The application domains of cybernetic control span engineering, biology, and psychology, with the goal of enabling robots to interact with humans in more intuitive ways. This involves adapting their actions and responses based on human feedback and behavior, cultivating a more seamless and responsive human-robot interaction.

Cybernetic biorobotics, at its core, is an interdisciplinary frontier that harmonizes principles from cybernetics, biology, and robotics. Its primary mission is the exploration and development of robots or robotic systems intricately inspired by the marvels of biological organisms. This field is driven by the ambition to conceive robots capable of mimicking and integrating the sophisticated principles and behaviors observed in living entities. Researchers draw inspiration from the intricate control systems of biological organisms and this creative synthesis results in the creation of robots characterized by adaptive and intelligent behaviors, thus mirroring the intricacies found in the natural world. Bioinspired robotics, a central focus within this discipline, involves distilling the fundamental principles and behaviors intrinsic to biological entities and skillfully incorporating them into the design and control of robotic systems. It has the potential to advance the development of robots endowed with locomotion and manipulation capabilities akin to animals, as well as robots capable of adapting to dynamic environments or interacting with humans in more natural and intuitive ways. Moreover, research in cybernetic biorobotics can offer valuable insights into comprehending biological systems, fostering advancements in disciplines like neuroscience and biomechanics. Furthermore, remote cybernetic robots may rely on haptic systems as essential interfaces. A haptic system, characterized by its ability to provide users with a sense of touch or tactile feedback through force, vibration, or other mechanical means, comprises a haptic interface and a haptic rendering system. Collaboratively, these components simulate touch sensations, enabling users to engage with virtual or remote environments in a tactile manner.

This research introduces a control and sensing architecture that integrates a cybernetic scheme based on the recognition of electromyographic control signals, governing a range of locomotive behaviors in a robotic fish. Conceptually, the human operator receives feedback signals from the sensors of the bio-robotic avatar, conveying information about its remote environment. The proposed approach stands out due to its key features and contributions, which include:

• The exposition of an innovative conceptual cybernetic fish avatar architecture.
• The creation of an EMG data filtering algorithm, coupled with a method for extracting, classifying, and recognizing muscular patterns using a deep ANN, serves as a cybernetic interface for the governance of the fish avatar.
• The development of a fuzzy-based oscillation patterns generator (OPG) designed to generate periodic oscillation patterns around the fish’s caudal fin. These coordinated oscillations are decoupled into right and left step functions, specifically crafted to input into a lateral pair of electromagnetic coils, thereby producing undulating swimming motions of the robot fish.
• The conception of a bioinspired robotic fish mechanism is characterized by the incorporation of underactuated elements propelled by serial links featuring helical springs. This innovative design is empowered by a dual solenoid electromagnetic oscillator and a four-bar linkage, reflecting a novel approach to bioinspired robotics.
• The derivation of closed-form control laws for both the undulation of the underactuated caudal multilink dynamics and the ballasting system.

Section 2 provides a comprehensive discussion on the comparative analysis of the current state of the art. Section 3 provides a detailed description of the proposed architecture cybernetic system model. In Section 4 presents an approach for filtering electromyography (EMG) data and delves into an in-depth discussion of a classifier based on deep ANN for the recognition of hand-motion EMG stimuli patterns. Section 5 presents the development of a fuzzy-based oscillation patterns generator. Section 6 details the robot’s mechanism parts and its dynamic model. Section 7 focuses on the development of a feedback control for the fish’s ballasting system. Finally, Section 8 provides the concluding remarks of the research study.

2. Analysis of the state of the art

This section syntheses the relevant literature and provides insights into the current state of the art. Further, it aims to examine and evaluate existing research and advancements in the field. This brief analysis identifies and compares different aspects, providing a comprehensive overview including relevant research and advancements in the field about methodologies, and outcomes.

Multiple basic concepts of cybernetics [1] in the intersection of physics and control theory and molecular systems were presented in [2] speed-gradient approach to modeling the dynamics of physical systems is discussed. A novel research approach namely Ethorobotics proposes the use and development of advanced bioinspired robotic replicas as a method for investigating animal behaviour [3]. In the domain of telepresence and teleoperation, diverse systems and methodologies have been devised to facilitate remote control of robots [4]. One such system is the multi-robot teleoperation system based on a brain-computer interface, as documented by [6]. This system aims to enable individuals with severe neuromuscular deficiencies to operate multiple robots solely through their brain activity, thus offering telepresence via a thought-based interaction mode. A comprehensive review addressing existing teleoperation methods and techniques for enhancing the control of mobile robots has been presented by [12]. This review critically analyzes, categorizes, and summarizes existing teleoperation methods for mobile robots while highlighting various enhancement techniques that have been employed. It makes clear the relative advantages and disadvantages associated with these methods and techniques. The field of telepresence and teleoperation robotics has witnessed substantial attention and interest over the past decade [14], finding extensive applications in healthcare, education, surveillance, disaster recovery, and corporate/government sectors. In the specific context of underwater robots, gesture recognition-based teleoperation systems have been developed to enable users to control the swimming behavior of these robots. Such systems foster direct interaction between onlookers and the robotic fish, thereby enhancing the intuitive experience of human-robot interaction. Furthermore, efforts have been made to enhance the consistency and quality of robotic fish tails through improved fabrication processes, and target tracking algorithms have been developed to enhance the tracking capabilities of these robots [11]. The work [10] developed teleoperation for remote control of a robotic fish by hand gestures recognition. It allowed direct interaction between onlookers and the biorobot. Another notable system is the assistive telepresence system employing augmented reality in conjunction with a physical robot, as detailed in the work by [7]. This system leverages an optimal non-iterative alignment solver to determine the optimally aligned pose of the 3D human model with the robot, resulting in faster computations compared to baseline solvers and delivering comparable or superior pose alignments. The review presented in [20] analyses the progress of robot skin as multimodal sensing and machine perception for sensory feedback in feeling proximity, pressure, temperature for collaborative robot applications considering immersive teleoperation and affective interaction. The work [24] reported an advanced robotic avatar system designed for immersive teleoperation, having some key functions such as human-like manipulation and communication capabilities, immersive 3D visualization and transparent force-feedback telemanipulation. Suitable human-robot collaboration in medical application has been reported [21], where force perception is augmented to human operator during needle insertion on the soft tissue. Telepresence of mobile robotic systems may incorporate remote video transmission to steer the robot by seeing through its eyes remotely. The work [22] presented an overview including social application domains. Research has been conducted on the utilization of neural circuits to contribute to limb locomotion [15] in the presence of uncertainty. By optimizing data, it showed the combination of circuits necessary for efficient locomotion. A review has also been conducted on central pattern generators (CPGs) employed for locomotion control in robots [16]. This review encompasses neurobiological observations, numerical models, and robotic applications of CPGs. The work [27] describes an extended mathematical model of the CPG supported by two neurophysiological studies: identification of a two-layered CPG neural circuitry and a specific neural model for generating different patterns. The CPG model is used as the low-level controller of a robot to generate walking patterns. Inclusion of ANN as a layer of the CPG, to produce rhythmic and non-rhythmic motion patterns. The work in [18] presented a review of bionic robotic fish, tackling major concepts on kinematics and control, learning, hydrodynamic forces and critical concepts on locomotion coordination. The research presented in [19] reviews the human manual control of devices in cybernetics using mathematical models and advances of theory and applications, from linear time-invariant modeling of stationary conditions to methods and analysis of adaptive and time-varying of cybernetics-human in control tasks. New foundations for cybernetics will emerge and impact numerous domains involving humans in manual and neuromuscular systems modeling control.

Building upon the preceding analysis regarding relevant literature, the subsequent table (Table 1) encapsulates the primary distinctions articulated in this study in relation to the most pertinent literature identified.

As delineated in Table 1, the present study introduces distinctive elements that set it apart from the recognized relevant literature. However, it is noteworthy to acknowledge that various multidisciplinary domains may exhibit commonalities. Across these diverse topics, shared elements encompass robotic avatars, teleoperation, telepresence, immersive human-robot interfaces, as well as haptic or cybernetic systems in different application domains. In this research, the fundamental principle of a robotic avatar entails controlling its swimming response to biological stimuli from the human operator. The human controller is able to gain insight into the surrounding world of the robotic fish avatar through a haptic interface. This interface allows the human operator to yield biological electromyography stimuli as the result of their visual and skin impressions (e.g. pressure, temperature, heading vibrations). The biorobotic fish generates its swimming locomotive behavior that is governed by EMG stimuli yielded in real-time in the human. Through a neuro-fuzzy controller, the neuronal part (cybernetic observer) classifies the type of human EMG reaction, and the fuzzy part determines the swimming behavior.

3. Conceptual System Architecture

This section encompasses a comprehensive framework that highlight the integration of various components to propose a cohesive cybernetic robotic model. In addition, this section outlines the key concepts and elucidates their interactions within the system.

Figure 1 presents an overview of the key components constituting the proposed system architecture. This manuscript thoroughly explores the modeling of four integral elements: i) the cybernetic human controller, employing ANN classification of EMG signals; $ii$) a fuzzy-based locomotion pattern generator; $iii$) an underactuated bioinspired robot fish; and $iv$) the robot’s sensory system, contributing feedback for the haptic system. While we will discuss the relevance and impact of the latter item within the architecture, it is important to note that the detailed exploration of topics related to haptic development and wearable technological devices goes beyond the scope of this paper and will be addressed in future work. Nevertheless, we deduce the observable variables that serve as crucial inputs for the haptic system.

Essentially, there are six haptic feedback sensory inputs of interest for the human, representing the observable state of the avatar robot: Eulerian variables, including angular and linear displacements and their higher-order derivatives; biomechanical caudal motion; hydraulic pressure; scenario temperature; and passive vision. Figure 2-left provides an illustration of the geometric distribution of the sensing devices.

The instrumented robot is an embodiment of the human submerged in water, featuring an undulatory swimming mechanical body imbued with muscles possessing underactuated characteristics. These features empowers the bio-robotic avatar to execute movements and swim its aquatic surroundings.

The observation models aim to provide insights into how these sensory perceptions are conveyed to the haptic helmet with haptic devices, including wheel reaction mechanism. A comprehensive schema emerges wherein the haptic nexus, bolstered by pivotal human’s biosensorial components including gravireceptors, Ruffini corpuscles, Paccinian receptors, and retinal photoreceptors, converges to interface with the sensory substrate of the human operator. Such a convergence engenders a cascading sequence wherein biological input stimuli coalesce to yield discernible encephalographic activity, the primary layer of subsequent electromyographic outputs. These consequential EMG outputs, undergo processing a swim oscillatory pattern generator, thereby embodiment control of biomechanical cybernetic governance.

In accordance with Figure 1, it is noteworthy that among the various haptic input variables such as temperature, pressure, and the visual camera images represent direct sensory measurements transmitted from the robotic fish to the components of the haptic interface. Conversely, the robotic avatar takes on the role of a thermosensory adept, in order for the human to assess the ambient thermal landscape. Thus, from this discernment, a surrogate thermal approach is projected onto the tactile realm of the human operator, through modulation of thermally responsive plates enmeshed within the haptic interface. Therefore, a crafted replica of the temperature patterns detected by the robotic aquatic entity is seamlessly integrated into the human sensory experience. This intertwining of thermal emulation is reached by the network of Ruffini corpuscles, intricately nestled within the human skin, thereby enhancing the experiential authenticity of this multisensory convergence. As for interaction through the haptic functions, the Paccinian corpuscles function as discerning receptors, proficiently registering subtle haptic pressures. It finds its origin in the dynamic tactile signals inherent to the aquatic habitat, intricately associated with the underwater depth traversed by the robotic avatar. Integral to the comprehensive sensory scheme, the optical sensors housed within the robotic entity acquires visual data. This visual data is subsequently channeled to the human’s cognition, through the haptic interface’s perceptual canvas. Within this, the human sensory apparatus assumes the role of an engaged receptor, duly transducing these visual envoys through the lattice of retinal photoreceptors.

Embedded within the robotic fish’s body, several Inertial Measurement Units (IMU) plays a pivotal role in quantifying Eulerian inclinations intrinsic to the aquatic environment. These intricate angular displacements are subsequently channeled to the human operator, thereby initiating an engagement of the reaction wheel mechanism. As a consequential outcome of this interplay, a synchronized emulation of tilting motion is induced, mirroring the nuanced cranial adjustments executed by the human operator. Any alignment of movements assumes perceptible form, relayed through the humans’s network of gravireceptors nestled within the internal auditory apparatus. The Euler angular and linear speeds are not directly measured; instead, they must be integrated using various sensor fusion approaches to enhance the avatar’s fault tolerance in reading its environment. For example, the angular observations of the robot fish are obtained through numerical integro-differential equations, which are solved online as measurements are acquired. Let’s introduce the following notation for inclinometers (${\alpha }_i$) and accelerometers ($a_{\mathrm{a}}$), with the singular direct sensor measurement ${\dot{\alpha }}_g$ derived from the gyroscopes. The observation for the fish’s roll velocity combining the three inertial sensors is

$${\omega }_{\alpha }={\displaystyle \frac{\mathrm{d}{\alpha }_{\iota }}{\mathrm{d}t}}+\dot{{\alpha }_g}+\frac{1}{d_{\alpha }}{\int }_t{a}_{\alpha }\mathrm{d}t,$$

(1a)

while the pitch velocity is modelled by

$${\omega }_{\beta }={\displaystyle \frac{\mathrm{d}{\beta }_{\iota }}{\mathrm{d}t}}+\dot{{\beta }_g}+\frac{1}{d_{\beta }}{\int }_t{a}_{\beta }\mathrm{d}t,$$

(1b)

and the yaw velocity is obtained by

$${\omega }_{\gamma }={\displaystyle \frac{\mathrm{d}{\gamma }_{\iota }}{\mathrm{d}t}}+\dot{{\gamma }_g}+\frac{1}{d_{\gamma }}{\int }_t{a}_{\gamma }\mathrm{d}t.$$

(1c)

Within this context, the tangential accelerations experienced by the robot body are denoted as $a_{\alpha ,\beta ,\gamma }$ [m/s${ }^2$]. Additionally, the angular velocities measured by the gyroscopes are represented by $\dot{\alpha },\dot{\beta },\dot{\gamma }$ [rad/s${ }^2$]. Correspondingly, the inclinometers provide angle measurements denoted as $\alpha ,\beta ,\gamma$ [rad]. These measurements collectively contribute to the comprehensive observability and characterization of the robot’s dynamic behavior and spatial orientation.

Furthermore, the oscillations of the caudal tail are reflections of the dynamics of the underactuated spine. These dynamics are captured by quantifying encoder pulses, denoted as ${\eta }_t$, which provide precise angular positions for each vertebra. Given that real-time angular measurements of the vertebrae are desired, higher-order data is prioritized. Consequently, derivatives are computed by initiating from the Taylor series to approximate the angle of each vertebral element with respect to time, denoted as t.

$${\varphi }_i\approx {\displaystyle \frac{{\varphi }_i^{\left(0\right)}}{0!}}{(t_2-t_1)}^0+{\displaystyle \frac{{\varphi }_i^{\left(1\right)}}{1!}}{(t_2-t_1)}^1+{\displaystyle \frac{{\varphi }_i^{\left(2\right)}}{2!}}{(t_2-t_1)}^2+\cdots$$

(2a)

thus, rearranging Math notation and trunking up to the first-derivative,

$${\varphi }_i\approxeq {\varphi }_i+{\varphi }_i^{\left(1\right)}(t_2-t_1),$$

(2b)

dropping ${\varphi }_i^{\left(1\right)}$ off as a state variable, the first-order derivative (${\varphi }^{\left(1\right)}\left(t\right)\equiv \dot{\varphi }\left(t\right)$) is given by

$$\dot{\varphi }\left(t\right)={\displaystyle \frac{{\varphi }_2-{\varphi }_1}{t_2-t_1}},$$

(2c)

and assuming that a vertebra’s angular measurement model and in terms of the encoder’s pulses $\eta$ with resolution R, then it is stated that by substituting the pulses encoder model into the angular speed function for the first vertebra,

$${\dot{\varphi }}_1=\left({\displaystyle \frac{2\pi }{R(t_2-t_1)}}{\eta }_2\right)-\left({\displaystyle \frac{2\pi }{R(t_2-t_1)}}{\eta }_1\right)=\left(\frac{2\pi }{R}\right)\left({\displaystyle \frac{{\eta }_2-{\eta }_1}{t_2-t_1}}\right),$$

(3a)

as for the second vertebra,

$${\dot{\varphi }}_2={\dot{\varphi }}_1+{\displaystyle \frac{2\pi }{R}}\left({\displaystyle \frac{{\eta }_2-{\eta }_1}{t_2-t_1}}\right)$$

(3b)

$${\dot{\varphi }}_3={\dot{\varphi }}_1+{\dot{\varphi }}_2+{\displaystyle \frac{2\pi }{R}}\left({\displaystyle \frac{{\eta }_2-{\eta }_1}{t_2-t_1}}\right)$$

(3c)

and

$${\dot{\varphi }}_4={\dot{\varphi }}_1+{\dot{\varphi }}_2+{\dot{\varphi }}_3+{\displaystyle \frac{2\pi }{R}}\left({\displaystyle \frac{{\eta }_2-{\eta }_1}{t_2-t_1}}\right)$$

(3d)

The preliminary sensing models serve as a comprehensive representation, strategically integrated into the control models as crucial feedback terms. A detailed exploration of this integration is elucidated in Section 6 and Section 7.

4. Deep ANN-based EMG Data Classification

This section details the experimental acquisition of EMG data, its spatial filtering, and pattern extraction achieved through the statistical combination of linear envelopes. Additionally, an adaptive method for class separation and data dispersion reduction is described. The section also covers the structure of a deep neural network, presenting its classification output results from mapping input EMG stimuli.

The related work reported a system for automatic pattern generation for neurosimulation in [25], where a neurointerface was used as a neuro-protocol for outputting fingers deflection and nerves stimulation. In the present research, numerous experiments were carried out to pinpoint the optimal electrode placement and achieve precise electromyographic readings for each predefined movement in the experiment. The positions of the electrodes were systematically adjusted, and the results from each trial were compared. Upon data analysis, it was discerned that the most effective electrode placement is on the Ulnar nerve, situated amidst the muscles Flexor Digitorium Superficialis, Flexor Digitorium Profundus, and Flexor Carpi Ulnaris. A series of more than ten experiments was executed for each planned stimulus or action involving hands, allowing a 2s interval between each action, including the opening and closing of hands, as well as the extension and flexion of the thumb, index, middle, ring and little fingers. The data was measured by a g.MOBIlab+ device with two-channel electrodes and quantified in microvolts per second [$\mu v/s$], as depicted in Figure 3.

The data acquired from the electromyogram often exhibit substantial noise attributed to both the inherent nature of the signal and external vibrational factors. To refine the data quality by mitigating this noise, a filtering process is essential.

In this context, a second-order Notch filter was utilized. This filter is tailored to target specific frequencies linked to noise, proving particularly effective in eliminating electrical interferences and other forms of stationary noise [26]. A Notch filter is a band-rejection filter to greatly reduce interference caused by a specific frequency component or a narrow band signal. Hence, in the Laplace space, the second-order filter is represented by the analog Laplace domain transfer function:

$$H\left(s\right)=\frac{s^2+{\omega }_o^2}{s^2+2s\xi {\omega }_o+{\omega }_o^2},$$

(4)

where ${\omega }_0$ signifies the cut angular frequency targeted for elimination, and $2\xi$ signifies the damping factor or filter quality, determining the bandwidth. Consequently, solving to obtain its solution in the physical variable space. The bilinear transformation relates the variable s from the Laplace domain to the variable $z_k$ in the $\mathbb{Z}$ domain, considering T as the sampling period, and is defined as follows,

$$s\doteq \left({\displaystyle \frac{1}{T}}\right)\left({\displaystyle \frac{z_k-1}{z_k+1}}\right),$$

(5)

upon substituting the previous expression into the transfer function of the analog Notch filter and algebraically simplifying, the following transfer function in the $\mathbb{Z}$ domain is obtained, redefining the notation as ${\upsilon }_t=z_k$ just to meet equivalence with the physical variable,

$$\mathfrak{h}\left({\upsilon }_t\right)=\frac{1-2cos\left({\omega }_{0}t\right){\upsilon }_t^{-1}+{\upsilon }_t^{-2}}{1-2\xi cos\left({\omega }_{0}t\right){\upsilon }_t^{-1}+\xi {\upsilon }_t^{-2}}.$$

(6)

The 2^nd-order Notch filter $h\left(\upsilon \right)$ was employed on raw EMG data to alleviate noise resulting from electrical impedance and vibrational electrode interference, with parameters set at ${\omega }_0=256$Hz and $\xi =0.1$ and results depicted in Figure 4.

Following this, leveraging the outcomes of the Notch filter, the data undergoes processing through three distinct filters or linear envelopes. This serves as a secondary spatial filter and functions as a pattern extraction mechanism. These include a filter for average variability, one for linear variability, and another for average dispersion. Each filter serves a specific purpose, enabling the analysis of different aspects of the signal. Consider n as the number of measurements constituting a single experimental stimulus, and let N represent the entire sampled data space obtained from multiple measurements related to the same stimulus. Furthermore, denote ${\widehat{v}}_i$ as the $i^{th}$ EMG measurement of an upper limb, measured in microvolts ($\mu V$). From such statements, the following Propositions 1, 2 and 3 are introduced as new data patterns.

Proposition 1

(Filter $\gamma$). The γ pattern refers to a statistical linear envelope described by the difference of a local mean ${\widehat{v}}_k$ in a window of samples and the statistical mean ${\widehat{\upsilon }}_i$ of all samples population an experiment.

$$\gamma \left(v_k\right)=\left|{\widehat{v}}_i-\frac{1}{n}\sum _{k=1}^n{v}_k\right|.$$

(7)

Proposition 2

(Filter $\lambda$). The $\lambda \left(v_k\right)$ pattern refers to an statistical linear envelope denoted by the difference of a local mean ${\widehat{v}}_k$ in a window of samples and the statistical mean ${\widehat{\upsilon }}_i$ of the whole population of experiments of same type,

$$\lambda \left(v_k\right)=\left|{\widehat{v}}_i-\frac{1}{N_k}\sum _{k=1}^{N_k}\widehat{v_k}\right|.$$

(8)

Proposition 3

(filter $\Omega$). The Ω pattern refers to an statistical linear envelope denoted by the difference of statistical means between population of one experiment ${\widehat{v}}_i$ and the whole population of numerous experiments of the same type:

$$\Omega \left(v_k\right)=\left|\frac{1}{n}\sum _{k=1}^n{v}_k-\frac{1}{N}\sum _{k=1}^N{v}_k\right|.$$

(9)

Hence, let the vector $\overrightarrow{\mathbf{\delta }}\in {\mathbb{R}}^3$ such that ${\mathbf{\delta }}_k={({\Omega }_k,{\gamma }_k,{\lambda }_k)}^{\top }$ and represent filtered data points in the $\Omega \gamma \lambda$-space. This work considered to include a brief data preprocessing as a method to improve multi-class separability and data scattering reduction in pattern extraction. Three distinctive patterns—$\gamma ,\lambda ,\Omega$—captivatingly converge in Figure 5. This illustration exclusively features patterns associated with sequences of muscular stimuli from both the right and left hands. For supplementary stimulus plots, refer to the appendix in Appendix A at the end of this manuscript.

From numerous laboratory experiments, over $75\%$ of the sampled raw data fall within the range of one standard deviation. Consider the vector $\overrightarrow{\mathbf{\sigma }}\in {\mathbb{R}}^3$ such that the standard deviations vector $\overrightarrow{\mathbf{\sigma }}={({\sigma }_{\Omega },{\sigma }_{\gamma },{\sigma }_{\lambda })}^{\top }$ encompasses the three spatial components by its norm.

$$\mathbf{\sigma }=\sqrt[2]{\frac{1}{N}\sum _{k=1}^N{\left(\left(\begin{array}{c}{\Omega }_k\\ {\gamma }_k\\ {\lambda }_k\end{array}\right)-\left(\begin{array}{c}{\mu }_{\gamma }\\ {\mu }_{\lambda }\\ {\mu }_{\Omega }\end{array}\right)\right)}^2}.$$

(10)

Building upon the preceding statement, we can formulate an adaptive discrimination criterion as elucidated in Definition 1.

Definition 1

(Discrimination condition). Consider the scalar value ${\mathbf{\delta }}_j$ as preprocessed EMG data located within a radius of magnitude ${\kappa }_d$ times the standard deviation $\parallel \mathbf{\sigma }\parallel$:

$${\mathbf{\delta }}_j=\left\{\begin{array}{cc}{\mathbf{\delta }}_k,\hfill & {\kappa }_d\parallel \mathbf{\sigma }\parallel \le \parallel \overrightarrow{\mathbf{\delta }}\parallel \hfill \\ \mathbf{0},\hfill & {\kappa }_d\parallel \mathbf{\sigma }\parallel >\parallel \overrightarrow{\mathbf{\delta }}\parallel \hfill \end{array}\right.$$

(11)

where $\mathbf{0}={(0,0,0)}^{\top }$ represents discriminated data.

Hence, considering the recent Definition 1, and in the current scenario with ${\kappa }_d=1.0$, which serves as a tuning discrimination factor. Therefore, the norm $l_h$ represents the distance between the frame origin and any class in the $\Omega \gamma \lambda$-space. This distance is adaptively calculated based on the statistics of each EMG class.

$$l_h=\sqrt[2]{{\left({\kappa }_{\Omega }{\sigma }_{\Omega }\right)}^2+{\left({\kappa }_{\gamma }{\sigma }_{\gamma }\right)}^2+{\left({\kappa }_{\lambda }{\sigma }_{\lambda }\right)}^2}$$

(12)

where, the coefficients ${\kappa }_{\Omega }$, ${\kappa }_{\gamma }$, ${\kappa }_{\lambda }$ are smooth adjustment parameters to set separability along axes. Hence, relocating each class center to a new position is stated by Proposition 4.

Proposition 4

(Class separability factor). New class position ${\mu }_{\Omega ,\gamma ,\lambda }^{+}$ in the $\Omega \gamma \lambda$-space, is established by the statistically adaptive linear relationship:

$${\mu }_{\Omega }^{+}={\mu }_{\Omega }+{\zeta }_{\Omega }{l}_h,$$

(13a)

$${\mu }_{\sigma }^{+}={\mu }_{\gamma }+{\zeta }_{\gamma }{l}_h$$

(13b)

and

$${\mu }_{\lambda }^{+}={\mu }_{\lambda }+{\zeta }_{\lambda }{l}_h.$$

(13c)

where ${\zeta }_{\Omega \gamma \lambda }$ are coarse in-space separability factors. The mean values ${\mu }_{\Omega \gamma \lambda }$ are the actual class positions obtained from the linear envelopes $\Omega \left({\upsilon }_k\right)$, $\gamma \left({\upsilon }_k\right)$ and $\lambda \left({\upsilon }_k\right)$.

Thus, by following the step-by-step method outlined earlier, Figure 6 showcases the extracted features of the EMG data, representing diverse experimental muscular stimuli. These results hold notable significance in the research, as they successfully achieve the desired class separability and data scattering, serving as crucial inputs for the multilayer ANN.

Henceforth, the focus lies in identifying and interpreting the EMG patterns projected in the $\gamma \lambda \Omega$-space, as illustrated in Figure 6. The subsequent part of the section delve into the architecture and structure of the deep ANN employed as a classifier, providing a detailed account of the training process. Additionally, the section highlights the performance metrics and results achieved by the classifier, offering insights into its effectiveness. Despite the challenges posed by nonlinearity, multidimensionality, and extensive datasets, various neural network structures were configured and experimented with. These configurations involved exploring different combinations of hidden layers, neurons, and the number of outputs in the ANN. To achieve the highest success rate in accurate data classification through experimentation, the final ANN was designed with perceptron units, as depited in Figure 7. It featured three inputs corresponding to the three EMG patterns $\gamma ,\lambda ,\Omega$ and included 12 hidden layers, each with 20 neurons. The supervised training process, conducted on a standard-capability computer, took approximately 20–30 minutes, resulting in nearly 1% error in pattern classification.

However, in the initial stages of the classification computations, with some mistuned adaptive parameters, the classification error was notably higher, even with much deeper ANN structures, such as 100 hidden layers with 99 neurons per layer. To facilitate the implementation, this work utilized the C/C++ library Fast Artificial Neural Networks (FANN), generating extremely fast binary code once the ANN was trained. In the training process of this research, about 50 datasets from separate experiments for each type of muscular stimulus were collectively stored, each comprising several thousand muscular repetitions. A distinct classification label was assigned a priori for each class type within the patterns space. To demonstrate the reliability of the approach, 16 different stimuli per ANN were established for classification and recognition, resulting in the ANN having 4 combinatory outputs, each with two possible states. Figure 8 depicts mixed sequences encompassing all types of EMG stimuli, with the ANN achieving a 100% correct classification rate.

Moreover, Table 2 delineates the mapping relationship between the ANN’s input, represented by the EMG stimuli, and the ANN’s output linked to a swimming behavior for controlling the robotic avatar.

5. Fuzzy-based Oscillation patterns generator

This section delineates the methodology utilized to produce electric oscillatory signals, essential for stimulating the inputs of electromagnetic devices (solenoids) embedded within the mechanized oscillator of the bio-robotic fish. The outlined approach for generating periodic electric signals encompasses three key components: a) the implementation of a fuzzy controller; b) the incorporation of a set of periodic functions dictating angular oscillations to achieve desired behaviors in the caudal undulation of the fish; and c) the integration of a transformation model capable of adapting caudal oscillation patterns into step signals, facilitating the operation of the dual-coil electromagnetic oscillator.

Another study [28] reported a different approach, a neuro-fuzzy-topological biodynamical controller for muscular-like joint actuators. In the present research, an innovative strategy suggested for the fuzzy controller involves the combination of three distinct input fuzzy sets: the artificial neural network outputs transformed into Crisp sets, and the linear and angular velocities of the robot derived from sensor measurements. Simultaneously, the fuzzy outputs correspond to magnitudes representing the periods of time utilized to regulate the frequency and periodicity of the caudal oscillations. This comprehensive integration enables the fuzzy controller to effectively process both neural network-derived information and real-time sensor data, dynamically adjusting the temporal parameters that govern the fish’s undulatory motions.

The depiction of the outputs from the EMG pattern recognition neural network is outlined in Table 2. Each binary output in the table is linked to its respective crisp-type input fuzzy sets when represented in the decimal numerical base $y_C$, as illustrated in Figure 9a. Moreover, Definition 2 details the parametric nature of the input sets.

Definition 2

(Input fuzzy sets). The output of the artificial neural network (ANN) corresponds to the fuzzy input, denoted as $y_C$, and only when it falls within the crisp set $\mathcal{C}=\left[0,1,\cdots ,15\right]$, the membership in the crisp set is referred to as ${\mu }_y\left(y_C\right)$.

$${\mu }_y\left(y_C\right)=\left\{\begin{array}{cc}0,\hfill & y_C\notin \mathcal{C}\hfill \\ 1,\hfill & y_C\in \mathcal{C}.\hfill \end{array}\right.$$

(14a)

In relation to the sensor observations of the bio-robot, the thrusting velocity v [cm/s] and angular velocity ω [rad/s] exhibit S-shaped sets modeled by Sigmoid membership functions. This modeling approach is applied consistently to both types of input variables, capturing their extreme-sided characteristics. Let ${\mu }_{s,f}\left(v\right)$ define the sets labeled as ‘stop’and ‘fast’in relation to the thrusting velocity v are elucidated by,

$${\mu }_{s,f}\left(v\right)=\frac{1}{1+{\mathrm{e}}^{\pm v\mp a}}$$

(14b)

Likewise, for the sets designated as ‘left-fast’($lr$) and ‘right-fast’($rf$) concerning the angular velocity ω, the S-shaped sets are modeled as,

$${\mu }_{lf,rf}\left(\omega \right)=\frac{1}{1+{\mathrm{e}}^{\pm v\mp b}}$$

(14c)

In addition, the rest of the sets in-between are any of the $k^{th}$ Gauss membership functions (‘slow’, ‘normal’and ‘agile’), for the robot’s thrusting velocity with parametric mean-valued $\overline{v}$ and standard deviation ${\sigma }_{v_k}$

$${\mu }_k\left(v\right)={\mathrm{e}}^{-\left(\frac{{(\overline{v}-v)}^2}{2{\sigma }_{v_k}^2}\right)},$$

(14d)

and for its angular velocity (‘left-slow’, ‘no turn’and ‘right-slow’), with parametric mean-valued $\overline{\omega }$ and standard deviation ${\sigma }_{{\omega }_k}$,

$${\mu }_k\left(\omega \right)={\mathrm{e}}^{-\left(\frac{{(\overline{\omega }-v)}^2}{2{\sigma }_{{\omega }_k}^2}\right)},$$

(14e)

Therefore, the reasoning rules articulated in the inference engine have been devised by applying inputs derived from Table 2, specifically tailored to generate desired outputs that align with the oscillation periods T [s] of the fish’s undulation frequency (see Figure 10).

Definition 3

($v,\omega =any$). For any linguistic value v, representing sensor observations of the fish’s thrusting velocity,

v = any ≐ stop or slow or normal or agile or fast.

Likewise, for any linguistic value ω, representing sensor observations of the fish’s angular velocity,

ω = any ≐ left-fast or left-slow or no turn or right-slow or right-fast.

Therefore, the following inference rules describe the essential robot fish swimming behavior.

• if $y_C$=sink and $v=$any or $\omega =$any then $T_{f,r,l}=$too_slow
• if $y_C$=buoyant and $v=$any or $\omega =$any then $T_{f,r,l}=$too_slow
• if $y_C$=gliding and $v=$any and $\omega =$any then $T_f=$slow, $T_{r,l}=$too_slow
• if $y_C$=slow_thrust and $v=$any and $\omega =$any then $T_f=$slow, $T_{r,l}=$too_slow
• if $y_C$=medium_thrust and $v=$any and $\omega =$any then $T_f=$normal, $T_{r,l}=$too_slow
• if $y_C$= fast_thrust and $v=$any and $\omega =$any then $T_f=$agile, $T_{r,l}=$too_slow
• if $y_C$=slow-right_maneuvering and $v=$any and $\omega =$any then $T_{f,l}=$too_slow, $T_r=$normal
• if $y_C$= medium-right_maneuvering and $v=$any and $\omega =$any then $T_f=$slow, $T_r=$agile, $T_l=$too_slow
• if $y_C$=fast-right_maneuvering and $v=$any and $\omega =$any then $T_f=$normal, $T_r=$fast, $T_l=$too_slow
• if $y_C$=slow-left_maneuvering and $v=$any and $\omega =$any then $T_{f,r}=$too_slow, $T_l=normal$
• if $y_C$=medium-left_maneuvering and $v=$any and $\omega =$any then $T_f=$slow, $T_r=$too_slow, $T_l=$agile
• if $y_C$= fast-left_maneuvering and $v=$any and $\omega =$any then $T_f=$normal, $T_r=$too_slow, $T_l=$fast
• if $y_C$=speed-up_right-turn and $v=$any and $\omega =$any then $T_{f,l}=$too_slow, $T_r=$fast
• if $y_C$=speed-up_left-turn and $v=$any and $\omega =$any then $T_{f,r}=$too_slow, $T_l=$fast
• if $y_C$=slow-down_right-turn and $v=$any and $\omega =$any then $T_{f,l}=$too_slow, $T_r=slow$
• if y=slow-down_left-turn and $v=$any and $\omega =$any then $T_{f,r}=$too_slow, $T_l=$slow

Hence, the output fuzzy sets delineate the tail undulation speeds of the robot fish across the caudal oscillation period T [s], as illustrated in Figure 10. Notably, three identical output sets with distinct concurrent values correspond to the periods of three distinct periodic functions forward: (f), right (r) and left (l); as subsequently defined by equations (22a), (22b), and (22c).

The inference engine’s rules dictate the application of fuzzy operators to assess the terms involved in fuzzy decision-making. As for the thrusting velocity fuzzy sets, where $v=any$ was previously established, and by applying Definition 3 the fuzzy operator is described by

$${\mu }_v^{max}=\underset{{\mu }_k\in v}{max}\left({\mu }_{stop},{\mu }_{slow},{\mu }_{normal},{\mu }_{agile},{\mu }_{fast}\right).$$

Likewise, the angular velocity fuzzy sets, where previously $\omega =any$ was stated, by applying the second part of Definition 3, the following fuzzy operator is described by

$${\mu }_{\omega }^{max}=\underset{{\mu }_k\in \omega }{max}\left({\mu }_{left-fast},{\mu }_{left-slow},{\mu }_{noturn},{\mu }_{right-slow},{\mu }_{right-fast}\right).$$

Therefore, according to premise ($v=any$ or $\omega =any$), the following fuzzy operator applies,

$${\mu }_{v,\omega }^{max}=\underset{{\mu }_k\in v\cup \omega }{max}\left({\mu }_v^{max},{\mu }_{\omega }^{max}\right)$$

Essentially, the fuzzification process applies strictly similar for the rest of the inference rules, according to the following Proposition,

Proposition 5

(Combined rules ${\mu }_i^{*}$). The general fuzzy membership expression for the $i^{th}$ inference rule that combines multiple inference propositions is

$${\mu }_i^{*}=\underset{\mu \in y_C\cap v\cup \omega }{min}\left({\mu }_{y_i}\left(y_{C_i}\right)=1,{\mu }_k^{max}(v,\omega )\right).$$

(15)

In any Crisp set, each ${\mu }_i$ attains a distinct value of 1, irrespective of the corresponding inference rule indexed by i. This value aligns with the $i^{th}$ entry in the neural network outputs outlined in Table 2. Additionally, k represents a specific fuzzy set associated with the same input.

Executing the previously mentioned proposition, Remark Section 5 provides a demonstration of its application.

Remark 

(Proposition 5 example). Let us consider rule $i=1$, where $y_{C_1}=\text{’}\sin \mathrm{k}\text{’}$ and either ($v=10.0$ cm/s or $\omega =0.0$ rad/s). Thus, articulated in the context of the resulting fuzzy operator,

$${\mu }_1^{*}=\underset{{\mu }_k\in y_C\cap v\cup \omega }{min}\left({\mu }_{sink}\left(y_C\right),\underset{v\cup \omega }{max}{\mu }_{v,\omega }(1,1)\right)=\underset{{\mu }_k\in y_C\cap v\cup \omega }{min}\left(1,1\right)=1,$$

Here, ${\mu }_v\left(10.0\right)=1$ and ${\mu }_{\omega }\left(0.0\right)=1$. Based on the earlier proposition, the resulting ${\mu }_1^{*}=1$, and given that rule 1 indicates an output period $T=$’too-slow’, its inverse outcome $T\left({\mu }_1^{*}\right)=6.0$ seconds. This outcome is entirely accurate, because the fish’s swim undulation will slow down up to 6s which is the slowest period oscillation.

Moving forward, during the defuzzification process, the primary objective is to attain an inverse solution. The three output categories for periods T include ‘forward’, ‘right’, and ‘left’, all sharing identical output fuzzy sets of T (Figure 10). Nevertheless, the output fuzzy sets consist of two categories of distributions: Gauss and sigmoid distribution sets, as outlined in Definition 4. Regarding the Gauss distributions, their functional form is specified by:

Definition 4

(Output fuzzy sets). The membership functions for both extreme-sided output sets are defined as ‘fast’with ${\mu }_f$ and ‘too slow’with ${\mu }_{ts}$, such that

$${\mu }_{f,ts}\left(T\right)={\displaystyle \frac{1}{1+{\mathrm{e}}^{\pm T_k\mp c_k}}}.$$

(16)

Here, T [s] denotes the period of time for oscillatory functions, with the slope direction determined by its sign. The parameter c represents an offset, and k is the numerical index of a specific set.

Furthermore, the membership functions for three intermediate output sets are defined as ‘agile’with ${\mu }_a$, ‘normal’with ${\mu }_n$, and ‘slow’with ${\mu }_s$, such that:

$${\mu }_{a,n,s}\left(T\right)={\mathrm{e}}^{-\left(\frac{{({\overline{T}}_k-T)}^2}{2{\sigma }_{T_k}^2}\right)},$$

(17)

Here, ${\overline{T}}_k$ represents the mean value, and ${\sigma }_{T_k}$ denotes the standard deviation of set k, with k serving as the numeric index of a specific set.

In accordance with Definition 4, any ${\mu }_k$ possesses a normalized membership outcome within the interval ${\mu }_k\in [0,\cdots ,1]$. The inverse sigmoid membership function, denoted as $T_k\in \mathbb{R}$∀${\mu }_k$, is determined by the general inverse expression:

$$T\left({\mu }_k\right)=\mp ln\left({\displaystyle \frac{1}{{\mu }_k}}\pm c_k\right).$$

(18)

Similarly, the inverse Gaussian membership function, where $T_k\in \mathbb{R}$∀${\mu }_k$, is defined by the inverse function:

$$T\left({\mu }_k\right)={\mu }_k-\sqrt[2]{-2{\sigma }_k^{2}ln\left({\mu }_k\right)}.$$

(19)

Hence, exclusively for the $j^{th}$ category among the output fuzzy sets affected by the fuzzy inference rule essential for estimating the value of T, the Centroid method is employed for deffuzification through the following expression:

$$T_{f,r,l}={\displaystyle \frac{{\sum }_j{\mu }_j^{*}{T}_j\left({\mu }_j^{*}\right)}{{\sum }_j{\mu }_j^{*}}},$$

(20)

or more specifically,

$$T_{f,r,l}={\displaystyle \frac{{\mu }_f^{*}T\left({\mu }_f^{*}\right)+{\mu }_a^{*}T\left({\mu }_a^{*}\right)+{\mu }_n^{*}T\left({\mu }_n^{*}\right)+{\mu }_s^{*}T\left({\mu }_s^{*}\right)+{\mu }_{ts}^{*}\left({\mu }_{ts}^{*}\right)}{{\mu }_f^{*}+{\mu }_a^{*}+{\mu }_n^{*}+{\mu }_s^{*}+{\mu }_{ts}^{*}}},$$

(21)

For terms $T\left({\mu }_f\right)$ and $T\left({\mu }_{ts}\right)$, the inverse membership function (18) is applicable, whereas for the remaining sets in the $j^{th}$ category, the inverse membership (19) is applied.

The work [29] reported a CPG model to control a robot fish’s motion in swimming and crawling, and let perform different motions influenced by sensory input from light, water, and touch sensors. Oscillators and Central Pattern Generators (CPGs) are closely related concepts. Oscillators are mathematical or physical systems exhibiting periodic behavior and are characterized by the oscillation around a stable equilibrium point (limit cycle). In the context of CPGs, these are neural networks that utilize oscillators that create positive and negative feedback loops, allowing for self-sustaining oscillations and the generation of rhythmic patterns, particularly implemented in numerous robotic systems [30]. CPGs are neural networks found in the central nervous system of animals (e.g. fish swimming [31]), that generate rhythmic patterns of motor activity and are responsible for generating and coordinating optimized [32] repetitive movements.

The present research proposes a different approach from the basic CPG model, and as a difference from other wire-driven robot fish’s motion approaches [33], this study introduces three fundamental undulation functions: forward, right-turning, and left-turning. These functions are derived from empirical measurements of the robot’s caudal fin oscillation angles. However, a distinctive behavioral undulation swim is achieved by blending these three oscillation functions, each incorporating corresponding estimation magnitudes derived from the fuzzy controller outputs. The formulation of each function involves fitting empirical data through Fourier series. As a difference from other approaches on CPG parameters adjustment[34], the preceding fuzzy outputs obtained from (21) to estimate the time periods $T_f,T_r,T_l$ play a pivotal role in parameterizing the time periods for the periodic oscillation functions, as outlined in Proposition 6.

Proposition 6

(Oscillation patterns function). Three fundamental caudal oscillation patterns, designed to generate swimming undulations, are introduced, each characterized by 11 pre-defined numerical coefficients. These patterns are described by amplitude functions, denoted as $\psi (\varphi ,T)$, where ϕ represents oscillation angles, and the time period T is an adjustable parameter.

The undulation pattern for forward motion is provided by the following function

$$\begin{array}{cc}\hfill & {\psi }_f(\varphi ,T_f)=0.0997+0.3327cos\left(\varphi \frac{2\pi }{T_f}\right)-0.1297sin\left(\varphi \frac{2\pi }{T_f}\right)-0.5760cos\left(2\varphi \frac{2\pi }{T_f}\right)+\hfill \\ \hfill & 0.3701sin\left(2\varphi \frac{2\pi }{T_f}\right)-0.1431cos\left(3\varphi \frac{2\pi }{T_f}\right)+0.1055sin\left(3\varphi \frac{2\pi }{T_f}\right)-0.0870cos\left(4\varphi \frac{2\pi }{T_f}\right)+\hfill \\ \hfill & 0.06323sin\left(4\varphi \frac{2\pi }{T_f}\right)-0.0664cos\left(5\varphi \frac{2\pi }{T_f}\right)-0.0664sin\left(5\varphi \frac{2\pi }{T_f}\right).\hfill \end{array}$$

(22a)

Likewise, the undulation pattern for right-turn motion is given by the function,

$$\begin{array}{cc}\hfill & {\psi }_r(\varphi ,T_r)=0.3324+0.1915cos\left(\varphi \frac{2\pi }{T_f}\right)+0.0622sin\left(\varphi \frac{2\pi }{T}\right)+0.4019cos\left(2\varphi \frac{2\pi }{T}\right)+\hfill \\ \hfill & 0.2920sin\left(2\varphi \frac{2\pi }{T}\right)-0.264cos\left(3\varphi \frac{2\pi }{T}\right)-0.3634sin\left(3\varphi \frac{2\pi }{T}\right)-0.0459cos\left(4\varphi \frac{2\pi }{T}\right)\hfill \\ \hfill & -0.1413sin\left(4\varphi \frac{2\pi }{T}\right)+0.0665sin\left(5\varphi \frac{2\pi }{T}\right).\hfill \end{array}$$

(22b)

Finally, the undulation pattern for left-sided turning motion is established by expression:

$$\begin{array}{cc}\hfill & {\psi }_l(\varphi ,T_l)=-0.1994+0.125cos\left(\varphi \frac{2\pi }{T_l}\right)-0.0622sin\left(\varphi \frac{2\pi }{T_l}\right)+0.3354cos\left(2\varphi \frac{2\pi }{T_l}\right)-\hfill \\ \hfill & 0.292sin\left(2\varphi \frac{2\pi }{T_l}\right)-0.3305cos\left(3\varphi \frac{2\pi }{T_l}\right)+0.3634sin\left(3\varphi \frac{2\pi }{T_l}\right)-0.1124cos\left(4\varphi \frac{2\pi }{T_l}\right)+\hfill \\ \hfill & 0.1413sin\left(4\varphi \frac{2\pi }{T_l}\right)-0.0664cos\left(5\varphi \frac{2\pi }{T_l}\right)+0.2659sin\left(5\varphi \frac{2\pi }{T_l}\right).\hfill \end{array}$$

(22c)

The approaches to forward, right-turn, and left-turn based on the findings of Proposition 6 are illustrated in Figure 11. Additionally, a novel combined oscillation pattern emerges by blending these three patterns (23), each assigned distinct numerical weights through the neuro-fuzzy controller.

$$\psi (\varphi ,T_f,T_r,T_l)={\psi }_f(\varphi ,T_f)+{\psi }_r(\varphi ,T_r)+{\psi }_l(\varphi ,T_l).$$

(23)

The proposed robotic mechanism features a multilink-based propulsive spine, driven by an electromagnetic oscillator composed of a pair of antagonistic solenoids that necessitate a synchronized sequence of electric pulses (see Figure 12b). The amplitudes generated by $\psi (\varphi ,T_f,T_r,T_l)$ in Equation 23 essentially represent the desired undulation pattern for the robotic fish’s caudal fin. However, these oscillations are not directly suitable for the inputs of the coils. To address this, our work introduces a decomposition of $\psi$ into two step signals centered around a stable equilibrium point (limit cycle), one for the right coil (positive with respect to the limit cycle) and another for the left coil (negative with respect to the limit cycle). The coil’s step function, either for the right-sided or left-sided coil is given by $s_{r,l}$, taken the equilibrium point as their limit value $\xi$,

$$s_{r,l}=\left\{\begin{array}{cc}0,\hfill & \psi \le \xi \hfill \\ 1,\hfill & \psi >\xi \hfill \end{array}\right.$$

(24)

In contrast to the work presented in [35] focusing on the swimming modes and gait transition of a robotic fish, the current study, as depicted in Figure 11, introduces a distinctive context. The three oscillatory functions, ${\psi }_{f,r,l}$, are displayed both overlapped and separated, highlighting their unique decoupled step signals. Assuming $\xi =0$ for all $\varphi$ in each case, in Figure 11a, ${\psi }_r={\psi }_l\approx 0$, with $T_{r,l}\ge 6$; in Figure 11b, ${\psi }_f={\psi }_l\approx 0$, with $T_{f,l}\ge 6$; and in Figure 11c, ${\psi }_f={\psi }_r\approx 0$, with $T_{f,r}\ge 6$. For a more comprehensive understanding, Figure 12b presents the electromechanical components of the caudal motion oscillator.

6. Robot Fish Biomechanical Model

This section introduces the design of the robotic fish mechanism and explores the model of the underactuated physical system to illustrate the fish undulation motions. The conceptualization of the proposed system is inspired by an underactuated structure featuring a links-based caudal spine with passive joints, utilizing helical springs to facilitate undulatory locomotion (see Figure 12a). The robotic fish structure introduces a mechanical oscillator comprising a pair of solenoids activated through coordinated sequences of step signals, as described by (24). Essentially, the electromagnetic coils generate antagonistic attraction/repulsion linear motions, translating into rhythmic oscillations within a mechanized four-bar linkage (depicted in Figure 12b). This linkage takes on the form of a trapezoid, composed of two parallel rigid links and two lateral linear springs functioning as antagonistic artificial muscles. Moreover, beneath the electromagnetic oscillator, there is a ballasting device for either submersion or buoyancy (Figure 12c). The robot’s fixed reference system consists of the X axis, which intersects the lateral sides, and the Y axis aligned with the robot’s longitudinal axis.

In Figure 12a, the electromagnetic oscillator of the robotic avatar responds to opposing coordinated sequences of step signals. The right-sided (R) and left-sided (L) solenoids counteract each other’s oscillations, generating angular moments in the trapezoid linkage (first vertebra). Both solenoids are identical, each comprising a coil and a cylindrical neodymium magnet nucleus. The trapezoid linkage, depicted in Figure 12b, experiences magnetic forces $\pm f_{o{s}_{RL}}$ at the two neodymium magnet attachments situated at a radius of $r_{os}$, resulting in two torques, ${\tau }_{os}$ and ${\tau }_s$, with respect to their respective rotation centers. As input forces $\pm f_{o{s}_{RL}}$ come into play, the linear muscle in its elongated state stores energy. Upon restitution contraction, this stored energy propels the rotation of the link $r_s$, which constitutes the first vertebra of the fish.

Furthermore, the caudal musculoskeletal structure, comprising four links (${\ell }_1,{\ell }_2,{\ell }_3,{\ell }_4$) and three passive joints (${\theta }_1,{\theta }_2,{\theta }_3$), facilitates a sequential rotary motion transmitted from link 1 to link 4. This transmission is accompanied by an incremental storage of energy in each helical spring that is serially connected. Consequently, the last link (link 4) undulates with significantly greater mechanical advantage. In summary, a single electrical pulse in any coil is sufficient to induce a pronounced undulation in the swimming motion of the robot’s skeleton.

As for the ballasting control device situated beneath the floor of the electromechanical oscillator, activation occurs only when either of two possible outputs from the artificial neural network (ANN) is detected: when $y_C$ equals ’sink’ or ’buoyancy.’ However, the fuzzy nature of these inputs results in a gradual slowing down of the fish’s undulation to its minimum speed. Additionally, both actions are independently regulated by a dedicated feedback controller overseeing the ballasting device.

Now, assuming knowledge of the input train of electrical step signals $s_{r,l}$ applied to the coils, let us derive the dynamic model of the biorobot, starting from the electromagnetic oscillator and extending to the motion transmitted to the last caudal link. Thus, as illustrated in Figure 12ab, the force f [N] of the solenoid’s magnetic field oscillator is established on either side (right, denoted as R, or left, denoted as L),

$$f=\frac{B^{2}A}{2{\mu }_o},$$

(25)

in this context, A [m${ }^2$] represents the area of the solenoid’s pole. The symbol ${\mu }_o$ denotes the magnetic permeability of air, expressed as ${\mu }_o=4\pi \times {10}^{-7}$ H/m (henries per meter). Hence, the magnetic field B (measured in Teslas) at one extreme of a solenoid is approximated by:

$$B=\frac{{\mu }_{o}iN}{l},$$

(26)

where i represents the coil current [A], N is the number of wire turns in a coil, and l denotes the coil length [m]. Furthermore, a coil’s current is described by the following linear differential equation as a function of time t (in seconds), taking into account a potential difference $\mathrm{v}$ (in volts):

$$i=\frac{1}{L}{\int }_0^T\mathrm{v}\mathrm{d}t+i_0,$$

(27)

here, L represents the coil’s inductance (measured in henries, H) with an initial current condition denoted as $i_0$. Additionally, the coil’s induction model is formulated by:

$$L=\frac{{\mu }_o{N}^{2}A}{l}.$$

(28)

In essence, this study states that both lateral solenoids exhibit linear motion characterized by an oscillator force $f_{os}$. This force is expressed as:

$$f_{os}={\displaystyle \frac{i^{2}Al}{2{\mu }_o{N}^{2}A}}$$

(29)

and due to the linear impacts of solenoids at both sides R and L (refer to Figure 12a), the first bar of the oscillator mechanism generates a torque, expressed as:

$${\tau }_{os}=\left(f_{os}\right)\left(r_{os}\right),$$

(30)

it is theorized that the restitution/elongation force along the muscle is denoted as $f_m$ (R or L), with this force being transmitted from the electromechanical oscillator to the antagonistic muscle in the opposite direction (refer to Figure 12b). This implies that the force generated from the linear motion solenoid in the oscillator’s right-sided coil, denoted as $f_{o{s}_R}$, is applied at point R and subsequently reflected towards point L with an opposite direction, represented as $-f_{o{s}_L}$. Similarly, conversely from the oscillator’s left-sided solenoid, the force $f_{o{s}_L}$ is applied at point L and transmitted to point R as $-f_{o{s}_L}$. For $f_{o{s}_L}$ applied in L,

$$f_{m_R}={\displaystyle \frac{-f_{o{s}_L}}{cos\left({\alpha }_R\right)}};f_{m_L}={\displaystyle \frac{-f_{o{s}_R}}{sin\left({\alpha }_L\right)}},$$

(31)

Hence, the angles ${\alpha }_{R,L}$ assume significance as the forces acting along the muscles $f_m$ differ, resulting in distinct instant elongations $x_m\left(t\right)$. Consequently, the four-bar trapezoid-shaped oscillator mechanism manifests diverse inner angles, namely ${\theta }_{1,2}$, $\beta$, and ${\gamma }_{1,2}$ as illustrated in Figure 12b.

Thus, prior to deriving an analytical solution for $\alpha$, it is imperative to formulate a theoretical model for the muscle. In this study, a Hill’s model is adopted, as depicted in Figure 12a (on the right side). The model incorporates a serial element SE (overdamped), a contractile element CE (critically damped), and a parallel element PE (critically damped), each representing distinct spring-mass-damper systems.

The generalized model for the antagonistic muscle is conceptualized in terms of the restitution force, and it is expressed as:

$$f_m=f_{SE}-(f_{CE}+f_{PE}),$$

(32)

Therefore, by postulating an equivalent restitution/elongation mass $m_w$ associated with instantaneous weight-force loads w (such as due to hydrodynamic flows), the preceding model is replaced with Newton’s second law of motion,

$$f_m=m_w{\ddot{x}}_{SE}-m_w{\ddot{x}}_{CE}-m_w{\ddot{x}}_{PE}.$$

(33)

Furthermore, through the independent solution of each element within the system in terms of elongations, the SE model can be expressed as:

$$x_{SE}\left(t\right)=s_1{\mathrm{e}}^{{\lambda }_1}t+s_2{\mathrm{e}}^{{\lambda }_2}t,$$

(34)

Here, $s_{1,2}$ represent arbitrary constants representing damping amplitude. The terms ${\lambda }_{1,2}$ denote the root factors,

$${\lambda }_{1,2}=-\frac{c_{SE}}{m_w}\pm {\displaystyle \frac{\sqrt[2]{{\left(\frac{c_{SE}}{m_w}\right)}^2-4\frac{k_{SE}}{m_w}}}{2}}$$

(35)

Here, the factors ${\lambda }_{1,2}$ are expressed in relation to the damping coefficient $c_{SE}$ (in kg/s) and the elasticity coefficient $k_{SE}$ (in kg/s²).

Similarly, for the contractile element CE, its elongation is determined by:

$$x_{CE}\left(t\right)=(c_1+c_{2}t){\mathrm{e}}^{-\frac{c_{CE}}{m_w}t},$$

(36)

With amplitude factors $c_{1,2}$ and damping coefficient $c_{CE},$ a similar expression is obtained for the parallel element PE:

$$x_{PE}\left(t\right)=(p_1+p_{2}t){\mathrm{e}}^{-\frac{c_{PE}}{m_w}t},$$

(37)

with amplitude factors $p_{1,2}$ and damping coefficient $c_{PE},$ the next step involves substituting these functional forms into the general muscle model,

$$f_m\left(t\right)=m_w\left[\frac{d^2}{d{t}^2}{x}_{SE}\left(t\right)+\frac{d^2}{d{t}^2}{x}_{CE}\left(t\right)+\frac{d^2}{d{t}^2}{x}_{PE}\left(t\right)\right]$$

(38)

such that the complete muscle’s force model $f_m$ formulated by

$$\begin{array}{cc}\hfill f_m\left(t\right)=& \frac{s_1{\lambda }_1^2}{m_w}{\mathrm{e}}^{{\lambda }_{1}t}+\frac{s_2{\lambda }_2^2}{m_w}{\mathrm{e}}^{{\lambda }_{2}t}+c_1\frac{c_{CE}^2}{m_w}{\mathrm{e}}^{\frac{-c}{m_w}t}-c_2{c}_{CE}{\mathrm{e}}^{\frac{-c}{m_w}t}-c_2{c}_{CE}{\mathrm{e}}^{\frac{-c}{m_w}t}+\hfill \\ \hfill & c_2\frac{c_{CE}^2}{m_w}t{\mathrm{e}}^{\frac{-c}{m}t}+p_1\frac{c_{PE}^2}{m_w}{\mathrm{e}}^{\frac{-c}{m_w}t}-p_2{c}_{PE}{\mathrm{e}}^{\frac{-c}{m_w}t}-p_2{c}_{PE}{\mathrm{e}}^{\frac{-c}{m_w}t}+\hfill \\ \hfill & p_2\left(\frac{c_{PE}^2}{m_w}\right)t{\mathrm{e}}^{\frac{-c}{m_w}t}.\hfill \end{array}$$

(39)

Subsequently, simplifying the preceding expression leads to the formulation presented in Proposition 7.

Proposition 7

(Muscles force model). The solution to the muscle force model, based on a Hill’s approach, is derived as a time-dependent function $f_m\left(t\right)$ encompassing its three constituent elements (serial, contractile, and parallel). This formulation is expressed as:

$$\begin{array}{cc}\hfill f_m\left(t\right)=& \frac{s_1{\lambda }_1^2{\mathrm{e}}^{{\lambda }_{1}t}+s_2{\lambda }_2^2{\mathrm{e}}^{{\lambda }_{2}t}}{m_w}+\left(\frac{c_1+c_{2}t}{m_w}{c}_{CE}-2c_2\right)c_{CE}{\mathrm{e}}^{\frac{-c_{CE}}{m_w}t}+\hfill \\ \hfill & \left(\frac{p_1+p_{2}t}{m_w}{c}_{PE}-2p_2\right)c_{PE}{\mathrm{e}}^{\frac{-c_{PE}}{m_w}t}.\hfill \end{array}$$

(40)

Thus, without loss of generality, considering a muscle model characterized by elongation $x_m$ and a force-based model $f_m$, we proceed to derive the passive angles of the oscillator and the output forces $f_x$ and $f_y$ for ${\ell }_1$.

Under initial conditions, the trapezoid oscillator bars are assumed to have ${\theta }_0=0^o$, aligning the four-bar mechanism with the X axis. As the bars rotate by an angle ${\theta }_1$ due to solenoid impacts at points R or L, the input bar of the oscillator with a radius of $r_{os}$ undergoes an arc displacement $s_1$. Simultaneously, the output bar of shorter radius $r_s$ experiences a displacement rate of $s_2$, such that:

$$s_2=\left({\displaystyle \frac{r_S}{r_{os}}}\right)s_1.$$

(41)

the arc displacement at points R or L is given by $s_1=r_{os}{\theta }_{os}$. Consequently, the rotation angle of the input oscillator is expressed as:

$${\theta }_{os}={\displaystyle \frac{s_1}{r_{os}}}.$$

(42)

Therefore, by formulating this relationship in the context of forces and subsequently substituting the newly introduced functions, the resulting expression is

$${\theta }_{os}={\displaystyle \frac{1}{r_{os}}}\int {\int }_t\ddot{y}d{t}^2$$

(43)

Here, $\ddot{y}$ denotes the linear acceleration of either point R or L along the robot’s Y axis. By replacing the solenoid’s mass-force formulation,

$${\theta }_{os}={\displaystyle \frac{1}{m{r}_{os}}}\int {\int }_t{f}_{os}d{t}^2={\displaystyle \frac{f{t}^2}{2m{r}_{os}}}.$$

(44)

Hence, the functional expression for $s_1$ takes the form

$$s_2={\displaystyle \frac{r_s{f}_{os}{t}^2}{2m{r}_{os}^2}}.$$

(45)

Without loss of generality, the inner angle ${\theta }_1$ of the oscillator mechanism (refer to Figure 12b) is derived as:

$${\theta }_1={\displaystyle \frac{\pi }{2}}\pm {\theta }_{os}.$$

(46)

Initially, when the oscillator bars are aligned with respect to the X axis, an angular displacement denoted by ${\theta }_1$ occurs as a result of the transfer of motion from the solenoid’s tangential linear motion to the input bar. Similarly, in the output bar, the corresponding angular displacement is represented by ${\theta }_2$,

$${\theta }_2=\pi \pm (\theta +{\Delta }_{\theta }),$$

(47)

here, ${\Delta }_{\theta }$ signifies a minute variation resulting from motion perturbation along the various links of the caudal spine. The selection of the ± operator depends on the robot’s side, whether it is denoted as R or L. As part of the analysis strategy, the four-bar oscillator was geometrically simplified to half a trapezoid for the purpose of streamlining deductions (refer to Figure 12b). Within this reduced mechanism, two triangles emerge. One triangle is defined by the parameters $r_{os}$, ℓ, d, while the other is characterized by $x_m\left(t\right),\ell ,r_s$, where ℓ serves as the hypotenuse and the sides d, $r_{os}$, and $r_s$ remain constant. Consequently, the instantaneous length of the hypotenuse is deduced as follows:

$${\ell }^2=r_{os}^2+d^2-2d{r}_{os}cos\left({\theta }_1\right),$$

(48)

Upon determining the value of ℓ, the inner angle ${\gamma }_1$ can be derived as follows:

$$r_{o{s}^2=d^2+{\ell }^2-2\ell dcos\left({\gamma }_1\right)}$$

(49)

Therefore, by isolating ${\gamma }_1$,

$${\gamma }_1=arccos\left({\displaystyle \frac{r_{os}-{\ell }^2-d^2}{-2d\ell d}}\right).$$

(50)

Until this point, given the knowledge of ${\gamma }_1$ and ${\theta }_2$, it is feasible to determine the inner complementary angle ${\gamma }_2$ through the following process:

$${\gamma }_2={\theta }_2-{\gamma }_1,$$

(51)

Subsequently, the angle formed by the artificial muscle and the output bar can be established according to the following principle:

$${\displaystyle \frac{sin\left({\gamma }_2\right)}{x_m}}={\displaystyle \frac{sin\left(\beta \right)}{\ell }},$$

(52)

thus, the inner angle $\beta$ is

$$\beta =arcsin\left(\frac{\ell }{x_m}sin\left({\gamma }_2\right)\right),$$

(53)

or alternatively, an approximation of the muscle length is

$$x_m=\ell {\displaystyle \frac{sin\left({\gamma }_2\right)}{sin\beta }},$$

(54)

this is the mechanism through which the input bar transmits a force $f_{m_1}$, as defined in expression (31), from the tangent $f_{os}$ to the output bar, achieving a mechanical advantage denoted as $f_{m_2}$,

$$f_{m_2}=\left({\displaystyle \frac{r_{os}}{r_s}}\right)f_{m_1}.$$

(55)

Hence, in accordance with the earlier stipulation in expression (31), Definition 5 delineates the instantaneous angles ${\alpha }_{R,L}$.

Definition 5

(Angles ${\alpha }_{R,L}$). The instantaneous angle α, expressed as a function of the inner angles of the oscillator, is introduced by:

$${\alpha }_{R,L}={\beta }_{R,L}-{\theta }_{1_{R,L}}.$$

(56)

It is noteworthy that, owing to the inertial system of the robot, the longitudinal force output component $f_y$ aligns with the input force $f_{os}$ in direction. Consequently, for a right-sided force, we have ${\alpha }_R\doteq {\beta }_R-{\theta }_{1_R}$, where:

$$f_{x_R}=\left(\frac{r_{os}}{r_s}\right)f_{o{s}_L}\frac{sin\left({\alpha }_R\right)}{cos\left({\alpha }_R\right)}$$

(57a)

and

$$f_{y_R}=\left(\frac{r_{os}}{r_s}\right)f_{o{s}_R}.$$

(57b)

Likewise, for the left-sided ${\alpha }_L\doteq {\beta }_L-{\theta }_{1_L}$,

$$f_{x_L}=\left(\frac{r_{os}}{r_s}\right)f_{o{s}_R}\frac{sin\left({\alpha }_L\right)}{cos\left({\alpha }_L\right)}$$

(58a)

as well as

$$f_{y_L}=\left(\frac{r_{os}}{r_s}\right)f_{o{s}_L}.$$

(58b)

In this scenario, an inverse solution is only applicable for $f_{x_{R,L}}$, with no necessity for determining $f_{y_{R,L}}$. Consequently, the mechanical advantage transferred between the input and output bars can be expressed by a simplified coefficient.

$$\kappa \doteq \frac{r_{os}}{r_s}.$$

(59)

Furthermore, through the utilization of the following trigonometric identity,

$$\frac{sin(\beta -{\theta }_1)}{cos(\beta -{\theta }_1)}\equiv tan(\beta -{\theta }_1)$$

(60)

can substitute and streamline the ensuing system of nonlinear equations by solving them simultaneously. Additionally, let ${\theta }_1$ be defined as:

$${\theta }_{1_{R,L}}=\frac{\pi }{2}\pm \frac{r_s{f}_{o{s}_{R,L}}{t}^2}{2m{r}_{os}^2}.$$

(61)

Hence, the simultaneous nonlinear system is explicitly presented solely for the force components along the X axis:

$$f_{x_R}=\kappa f_{o{s}_L}tan({\beta }_R-{\theta }_{1_R})$$

(62a)

and

$$f_{x_L}=\kappa f_{o{s}_R}tan({\beta }_L-{\theta }_{1_L}).$$

(62b)

Therefore, for the numerical solution of the system, a multidimensional Newton-Raphson approach is employed as outlined in the provided solution:

$${\beta }_{R_{t+1}}={\beta }_{R_t}-\frac{f_{x_R}\frac{\partial f_{x_L}}{\partial {\beta }_L}-f_{x_L}\frac{\partial f_{x_L}}{\partial {\beta }_L}}{\frac{\partial f_{x_R}}{\partial {\beta }_R}\frac{\partial f_{x_L}}{\partial {\beta }_L}-\frac{\partial f_{x_R}}{\partial {\beta }_L}\frac{\partial f_{x_L}}{\partial {\beta }_R}}$$

(63a)

and

$${\beta }_{L_{t+1}}={\beta }_{L_t}-\frac{f_{x_L}\frac{\partial f_{x_R}}{\partial {\beta }_R}-f_{x_R}\frac{\partial f_{x_L}}{\partial {\beta }_R}}{\frac{\partial f_{x_R}}{\partial {\beta }_R}\frac{\partial f_{x_L}}{\partial {\beta }_L}-\frac{\partial f_{x_R}}{\partial {\beta }_L}\frac{\partial f_{x_L}}{\partial {\beta }_R}}.$$

(63b)

Thus, by defining all derivative terms to finalize the system,

$$\frac{\partial f_{x_R}}{\partial {\beta }_R}=\kappa f_{o{s}_R}\left(\frac{-{\theta }_{1_R}}{{cos}^2({\beta }_R-{\theta }_{1_R})}\right),$$

(64a)

$$\frac{\partial f_{x_R}}{\partial {\beta }_L}=0,$$

(64b)

$$\frac{\partial f_{x_L}}{\partial {\beta }_R}=0,$$

(64c)

$$\frac{\partial f_{x_L}}{\partial {\beta }_L}=\kappa f_{o{s}_L}\left(\frac{-{\theta }_{1_L}}{{cos}^2({\beta }_L-{\theta }_{1_L})}\right).$$

(64d)

Therefore, by subsequently organizing and algebraically simplifying,

$${\beta }_{R_{t+1}}={\beta }_{R_t}-\frac{f_{x_R}\frac{\partial f_{x_L}}{\partial {\beta }_L}}{\frac{\partial f_{x_R}}{\partial {\beta }_R}\frac{\partial f_{x_L}}{\partial {\beta }_L}}={\beta }_{R_t}+\frac{tan({\beta }_R-{\theta }_{1_R}){cos}^2({\beta }_R-{\theta }_{1_R})}{{\theta }_{1_R}}$$

(65a)

and

$${\beta }_{L_{t+1}}={\beta }_{L_t}-\frac{f_{x_L}\frac{\partial f_{x_R}}{\partial {\beta }_R}}{\frac{\partial f_{x_R}}{\partial {\beta }_R}\frac{\partial f_{x_L}}{\partial {\beta }_L}}={\beta }_{L_t}+\frac{tan({\beta }_L-{\theta }_{1_L}){cos}^2({\beta }_L-{\theta }_{1_L})}{{\theta }_{1_L}}.$$

(65b)

The objective is to achieve numerical proximity, aiming for ${\beta }_{t+1}\approx {\beta }_t$. Consequently, through this inverse solution, the lateral force components of the first spinal link, denoted by ${\mathbf{f}}_{\mathbf{x}}$, are intended to be estimated because they are perpendicular to the links and produce the angular moments at each passive joint.

$${\mathbf{f}}_{\mathbf{x}}=\left(\begin{array}{c}{f}_{x_R}\\ f_{x_L}\end{array}\right)=\kappa ·\left(\begin{array}{c}{f}_{o{s}_L}tan\left({\beta }_{R_{t+1}}-\frac{\pi }{2}+\frac{r_s{f}_{o{s}_L}{t}^2}{2m{r}_{os}^2}\right)\\ f_{o{s}_R}tan\left({\beta }_{L_{t+1}}-\frac{\pi }{2}-\frac{r_s{f}_{o{s}_R}{t}^2}{2m{r}_{os}^2}\right)\end{array}\right).$$

(66)

Thus, given that the torque of the trapezoid’s second bar is ${\tau }_s=f_s{r}_s$ (see Figure 12b), we establish a torque-angular moment equivalence, denoted as ${\tau }_s\equiv M_1$. Leveraging this equivalence and the prior understanding of the torque ${\tau }_s$ acting on the second bar of the trapezoid, mechanically connected to the first link ${\ell }_1$, we affirm their shared angular moment. Consequently, the general expression for the tangential force $f_k$ applied at the end of each link ${\ell }_k$ is:

$$f_k=\frac{M_k}{{\ell }_k}.$$

(67)

Yet, considering the angular moment $M_k$ for each helical-spring joint, supporting the mass of the successive links, let’s introduce equivalent inertial moments, starting with $I_{{\epsilon }_1}=I_1+I_2+I_3+I_4$. Subsequently, we define $I_{{\epsilon }_2}=I_2+I_3+I_4$, $I_{{\epsilon }_3}=I_3+I_4$, and finally, $I_4$. Thus, in the continuum of the caudal spine, the transmission of energy to each link is contingent upon the preceding joints, as established by:

$$M_1=I_{{\epsilon }_1}\ddot{{\theta }_1},M_2=I_{{\epsilon }_2}\ddot{{\theta }_2},M_3=I_{{\epsilon }_3}\ddot{{\theta }_3},M_4=I_4\ddot{{\theta }_4}.$$

(68)

Each helical spring, connecting pairs of vertebrae, undergoes an input force $f=-kx$, directly proportional to the angular spring deformation indicated by elongation x. Here, k [kg m${ }^2$/s${ }^2$] represents the stiffness coefficient. External forces result in an angular moment, given by $\tau =-k\theta$, where torque serves as an equivalent variable to angular momentum, such that $I\alpha =-k\theta$. Consequently, when expressing the formula as a linear second-order differential equation, we have:

$$\ddot{\theta }+\frac{k}{I}\theta ={\ddot{\theta }}_{L_k}.$$

(69)

Here, ${\ddot{\theta }}_{L_k}$ represents undulatory accelerations arising from external loads or residual motions along the successive caudal links, which are detectable through encoders and IMUs. Assuming an angular frequency ${\omega }^2=k/I$, a period $p=2\pi \sqrt{I/k}$, and moments of inertia expressed as $I_k=r_k^2{m}_k$, the general equation is formulated as follows:

$$M_k=I_{{\epsilon }_k}{\ddot{\theta }}_k,$$

(70)

where ${\ddot{\theta }}_k$ is replaced by the helical spring expression (69) to derive

$$M_k=I_{{\epsilon }_k}\left({\ddot{\theta }}_{L_k}-{\displaystyle \frac{k_k}{I_{{\epsilon }_k}}}\right).$$

(71)

By algebraically extending, omitting terms, and rearranging for all links in the caudal spine, we arrive at the following matrix-form equation:

$$\left(\begin{array}{c}{M}_1\\ M_2\\ M_3\\ M_4\end{array}\right)=\left(\begin{array}{cccc}{I}_{{\epsilon }_1}& 0& 0& 0\\ 0& I_{{\epsilon }_2}& 0& 0\\ 0& 0& I_{{\epsilon }_3}& 0\\ 0& 0& 0& I_{{\epsilon }_4}\end{array}\right)·\left(\begin{array}{c}{\ddot{\theta }}_{L_1}\\ {\ddot{\theta }}_{L_2}\\ {\ddot{\theta }}_{L_3}\\ {\ddot{\theta }}_{L_4}\end{array}\right)-\left(\begin{array}{c}{k}_1{\theta }_1\\ k_2{\theta }_2\\ k_3{\theta }_3\\ k_4{\theta }_4\end{array}\right).$$

(72)

Hence, in accordance to expression (67), the tangential forces exerted on all the caudal links of the robotic fish are delineated by the following expression:

$$\left(\begin{array}{c}{f}_{l1}\\ f_{l2}\\ f_{l3}\\ f_{l4}\end{array}\right)=\left(\begin{array}{cccc}\frac{I_{{\epsilon }_1}}{{\ell }_1}& 0& 0& 0\\ 0& \frac{I_{{\epsilon }_2}}{{\ell }_2}& 0& 0\\ 0& 0& \frac{I_{{\epsilon }_3}}{{\ell }_3}& 0\\ 0& 0& 0& \frac{I_{{\epsilon }_4}}{{\ell }_4}\end{array}\right)·\left(\begin{array}{c}{\ddot{\theta }}_{L_1}\\ {\ddot{\theta }}_{L_2}\\ {\ddot{\theta }}_{L_3}\\ {\ddot{\theta }}_{L_4}\end{array}\right)-\left(\begin{array}{cccc}\frac{k_1}{{\ell }_1}& 0& 0& 0\\ 0& \frac{k_2}{{\ell }_2}& 0& 0\\ 0& 0& \frac{k_3}{{\ell }_3}& 0\\ 0& 0& 0& \frac{k_4}{{\ell }_4}\end{array}\right)·\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\\ {\theta }_3\\ {\theta }_4\end{array}\right)$$

(73)

Alternatively, the last expression can be denoted as the following control law:

$$\mathbf{f}=\mathbf{M}{\ddot{\mathbf{\theta }}}_L-\mathbf{Q}{\mathbf{\theta }}_t\equiv {\displaystyle \frac{\mathbf{M}({\dot{\mathbf{\theta }}}_{t_{t1}}-{\dot{\mathbf{\theta }}}_{t_{t1}})}{t_2-t_1}}-\mathbf{Q}{\mathbf{\theta }}_t,$$

(74)

Where $\mathbf{f}={(f_{l1},f_{l2},f_{l3},f_{l4})}^{\top }$, $\mathbf{M}$ represents masses dispersion, and $\ddot{\mathbf{\theta }}={({\ddot{\theta }}_{L_1},{\ddot{\theta }}_{L_2},{\ddot{\theta }}_{L_3},{\ddot{\theta }}_{L_4})}^{\top }$ denotes the vector of angular accelerations for the caudal vertebrae, including external loads. Additionally, $\mathbf{Q}$ stands for the matrix of stiffness coefficients. Therefore, the inverse dynamics control law, presented in a recursive form, is:

$${\dot{\mathbf{\theta }}}_{L_{t+1}}={\dot{\mathbf{\theta }}}_{L_t}+{\displaystyle \frac{{\mathbf{M}}^{-1}}{t_2-t_1}}\left(\mathbf{f}+\mathbf{Q}{\mathbf{\theta }}_t\right).$$

(75)

Finally, for feedback control, both equations are simultaneously employed within a computational recursive scheme, and angular observations are frequently derived from sensors on both joints: encoders and IMUs.

7. Ballasting Control System

This section delineates the integration of the ballasting control system, crafted to complement the primary structure of the biorobot. It introduces the ballasting model-based control system, selectively activated in response to the artificial neural network’s (ANN) output, particularly triggered when the ANN signals ’sink’ or ’buoyancy.’ Figure 13a visually depicts the biorobot’s ballasting system, while Figure 13b presents a diagram illustrating the fundamental components of the hydraulic piston, crucial for control modeling.

The core operational functions of the ballasting device involve either filling its container chamber with water to achieve submergence or expelling water from the container to attain buoyancy. Both actions entail the application of a linear force for manipulating a plunger or hydraulic cylindrical piston, thereby controlling water flow through either suction or exertion. Consequently, the volume of the liquid mass fluctuates over time, contingent upon a control reference or desired level marked as H, along with quantifying a filling rate $u\left(t\right)$ and measuring the actual liquid level $h\left(t\right)$.

Hence, we can characterize the filling rate $u\left(t\right)$ as the change in volume V with respect to time, expressed as

$$\frac{dV\left(t\right)}{dt}=u\left(t\right),$$

(76)

and assuming a cylindrical plunger-chamber with radius r and area $A=\pi r^2$, the volume is expressed as

$$V\left(t\right)=Ah\left(t\right),$$

(77)

Here, $h\left(t\right)$ represents the actual position of the plunger due to the incoming hydraulic mass volume. Consequently, the filling rate can also be expressed as

$$u\left(t\right)=k(H-h(t\left)\right).$$

(78)

Consider k as an adjustment coefficient, and let H be the reference or desired filling level. The instantaneous longitudinal filling level is denoted as $h\left(t\right)$. By substituting the previous expressions into the initial equation (76), we derive the following first-order linear differential equation:

$$A\frac{dh\left(t\right)}{dt}=k(H-h\left(t\right)).$$

(79)

To solve the aforementioned equation, we employ the integrating factor method, such that

$$\dot{h}\left(t\right)+\frac{k}{A}h\left(t\right)=\frac{k}{A}H.$$

(80)

In this instance, the integrating factor is determined as follows:

$$\mu \left(t\right)=e^{\int \frac{k}{A}dt}=e^{\frac{kt}{A}}.$$

(81)

Thus, by applying the integrating factor, we effectively reduce the order of derivatives in the subsequent steps,

$$e^{\frac{kt}{A}}\dot{h}\left(t\right)+e^{\frac{kt}{A}}\frac{k}{A}h\left(t\right)=\frac{k}{A}H{e}^{\frac{kt}{A}},$$

(82)

Through algebraic simplification of the left side of the aforementioned expression, the following result is determined

$${\left(h\left(t\right)e^{\frac{kt}{A}}\right)}^{\prime }=\frac{k}{A}H{e}^{\frac{kt}{A}}.$$

(83)

Following this, by integrating both sides of the equation with respect to time,

$${\int }_t{\left(h\left(t\right)e^{\frac{kt}{A}}\right)}^{\prime }dt={\int }_t\frac{k}{A}H{e}^{\frac{kt}{A}}dt,$$

(84)

where the expression on the left side undergoes a transformation into

$$h\left(t\right)e^{\frac{kt}{A}}=\frac{k}{A}H{\int }_t{e}^{\frac{kt}{A}}dt,$$

(85)

and the right side of the equation, once solved, transforms into

$$h\left(t\right)e^{\frac{kt}{A}}=H{e}^{\frac{kt}{A}}+c.$$

(86)

Now, to obtain the solution for $h\left(t\right)$, it is isolated by rearranging the term $e^{\frac{kt}{A}}$:

$$h\left(t\right)=H+c{e}^{-\frac{kt}{A}}$$

(87)

For initial conditions where $h\left(t_0\right)=0$ indicates the plunger is completely inside the contained chamber at the initial time $t_0=0$ s, the integration constant c is determined as

$$0=H+c{e}^0.$$

(88)

Therefore, the value of c takes on $c=-H$, and substituting it into the previous obtained solution,

$$h\left(t\right)=H(1-e^{-\frac{kt}{A}}).$$

(89)

In addition, considering that the required force of the piston $f_e$ is hence given by:

$$f_e=m\left(t\right)\frac{dv}{dt}+f_k+{\rho }_{a}A,$$

(90)

where $f_k$ is the friction force of the piston in the cylindrical piston, and ${\rho }_{a}A$ refers to the water pressure at that depth over the piston’s entry area. The instantaneous mass considers the piston’s mass $m_e$ and the liquid mass of the incoming water $m_a$:

$$m\left(t\right)=m_e+m_a,$$

(91)

where the water density is ${\delta }_a=\frac{m_a}{V_a}$ and $V_a=\pi r^{2}h\left(t\right)$, thus completing the mass model:

$$m\left(t\right)=m_e+{\delta }_a\pi r^{2}H(1-e^{-\frac{kt}{A}}).$$

(92)

Therefore, the force required to pull/push the plunge device is stated by the control law given as

$$f_e=\left(m_e+{\delta }_a\pi r^{2}H(1-e^{-\frac{kt}{A}})\right)\frac{dv}{dt}+f_k+{\rho }_{a}A.$$

(93)

8. Conclusion and future work

In summary, this study introduces a cybernetic control approach integrating electromyography, haptic feedback, and an underactuated bio-robotic avatar fish. Human operators control the avatar fish using their muscular stimuli, eliminating the need for handheld apparatus. The incorporation of fuzzy control, combining EMG stimuli with motion sensor observations, has proven highly versatile in influencing the decision-making process governing the fish’s swimming behavior.

The implementation of a deep neural network achieved remarkable accuracy, surpassing 98%, in recognizing sixteen distinct electromyographic gestures. This underscores the system’s robustness, effectively translating human intentions into precise control commands for the underactuated robotic fish.

This manuscript reports results from experimental EMG data classification and recognition using a multilayered artificial neural network. The oscillation pattern generator provided real signals to an experimental prototype of the underactuated robotic fish with its electromagnetic oscillator. Additionally, the fuzzy controller and the fish’s dynamical control model were validated through computer simulations.

While the introduction of haptic feedback and interface is conceptual in the proposed architecture, it represents a promising avenue for future research, aiming to enhance remote operation with immersive experiences. The advancements showcased in this work hold significant potential for future applications in underwater exploration through immersive cybernetic control.