A double variable-length pendulum with counterweight mass, kinematic excitation and electromagnetic forcing

This study introduces a novel double variable-length cable pendulum model and experimental setup featuring elastic suspension and counterweight mass. Our main goal is to investigate the complex dynamics resulting from variable length’s impact on vibration frequency and amplitude. Through numerical simulations and experiments, we explore the system’s response to different external forces. Utilizing methods like phase plots, bifurcation diagrams, and Lyapunov exponents, we delve into nonlinear dynamics. We also use vision-based techniques to assess friction damping-related vibrations and magnetic field interactions. The results reveal diverse behaviors, including chaotic and periodic oscillations, shedding light on control functions and parameter relationships. The developed cable system captures intricate nonlinear dynamics and attains stable vibration modes, as confirmed by vision-based measurements. This platform can analyze and control irregular dynamics in systems with elastically suspended weights driven by motors or mobile cranes. Its nature, encompassing kinematic excitation, electromagnetic interactions, and sliding friction, allows for exploring complex nonlinear dynamics. The system’s capacity to modulate vibration frequencies contributes to mitigating persistent vibrations.


Introduction
This study presents a novel research setup, so far unprecedented.The objective of this investigations is to transition the concept of a variable-length pendulum with accompanying effects described in distinct sections to a laboratory level, thereby establishing the foundation for investigating nonlinear phenomena and intricate dynamic effects.In light of this, an initial, approximate physical description and the associated mathematical models are provided.In the subsequent part of this introduction, we elaborate on key research concepts associated with the described topic, namely, the applications, mathematical modeling of such dynamic systems, the significance of frictional effects, the excitation originating from the applied magnetic field, and finally, the mathematical modeling of this intricate system.
The variable-length pendulum is a physical concept governed by certain differential equations and functional principles and is often associated with parametric oscillations [17].A parametric oscillator can be thought of as a harmonic oscillator with time-dependent phys-ical features [34].These time-dependent variables are related to the resonance frequency or damping of the oscillator.
After conducting an extensive analysis of the existing studies on the variable-length pendulum across different models and methodologies, we have identified specific applications of variable-length pendulum systems.These pendulums have many applications, including energy harvesters, load-lifting equipment, and robotics [1,3,25,27,36,38,41,43,47], such as pendulum-like robots in [18,49], to mention but a few.Especially, as stated by Li et al. many control strategies are designed based on single pendulum crane models, which ignore the hooks mass and treat the payload as a mass point [23].However, the payload will swing around the hook when the payload is too large, or the hook mass cannot be directly ignored, which is called double pendulum effects, making the dynamic more complex.The underactuated dynamical system of a double pendulum gantry crane with the load hoisting or lowering can be controlled with the use a nonlinear coupled tracking anti-swing controller developed in [39].
The experimental study described in [27] involved the development of a testing apparatus coupled with a mathematical model for energy harvesting.This study investigated the dependency of the dynamical response on the pendulum length through both computational simulations and experimentation.The findings revealed that as the pendulum length decreases, the energy gain increases.Furthermore, the study suggested that employing multiple pendulums could enhance the overall dynamic response of the system.In contrary to the attempts of energy gain effects, Mahmoudkhani et al. demonstrated a different application by implementing a flexural pendulum absorber connected to a linear system [25].The system dynamics were simulated using a sequential mass-spring-damper model, and the velocity of a specific point on the balanced beam was described using differential equations.The study unveiled that the flexural behavior of the beam displayed potential local bifurcation modes, such as Neimark-Sacker and pitchfork bifurcations, resulting in an increased response rate for the system.In the same field, Krasilnikov et al. showcased that the auto-parametric absorber rigid pendulum surpassed the performance of the standard rigid pendulum in terms of durability [19].Consequently, auto-parametric pendulum absorption systems demon-strate their effectiveness, particularly at low excitation amplitudes.
A comprehensive set of numerical experiments investigating resonances in a spring pendulum system with variable length is presented in [32].In this study, the equations of motion establish a three-degree-offreedom model system.The equations governing the pendulum angle, elongation, and slider displacement were formulated using the second derivatives of the state variables.The Euler-Lagrange equation and the Rayleigh dissipation function of the pendulum were employed for the pendulum angle and slider displacement equations, respectively.Under resonance conditions, the variable-length spring pendulum, suspended from a periodically forced slider, exhibits phenomena like quasi-periodicity and chaotic motions.Moreover, the investigation revealed that in close proximity to resonance, the influence of body coupling on system dynamics could lead to unpredictable dynamical behaviors.Extensions of this research can be found in [30], where a discrete wave modulated step function of length is used to stimulate a parametric pendulum.This pendulum undergoes mathematical analysis akin to the approach outlined in [32], and the Lyapunov exponents of the system were computed in [40].
Considering that a variable-length pendulum can exhibit both faster and longer oscillations as shown in the research [19], this modification could lead to an associated enhancement in fluid pumping efficiency.Such an improvement aligns with the reduction in human effort resulting from faster and longer oscillations compared to a pendulum with a constant length, as highlighted in another study [45].As a result, it is suggested that, in certain engineering applications, a variable-length pendulum could be preferable over the conventional constant-length pendulum, offering a substantial reduction in the necessary human effort.
Let us summarize the advantages of diverse applications utilizing variable-length pendulum concepts.They can be found in various fields, including energy harvesting, load-lifting equipment, robotics, and pendulum-like robots; research has shown that as the pendulum length decreases, energy gain increases, making variable-length pendulums suitable for energy harvesting applications; the introduction of variablelength pendulums introduces complex dynamics, including double pendulum effects, enhancing the system's dynamic behavior; adaptive performance and durability compared to standard rigid pendulums, especially at low excitation amplitudes is shown on the basis of auto-parametric absorber pendulums; variable-length pendulums have demonstrated potential for enhancing fluid pumping efficiency due to the ability to exhibit both faster and longer oscillations; numerical experiments investigating resonance in spring pendulum systems with variable length reveal phenomena such as quasi-periodicity and chaos.
Various methods are employed to derive mathematical equations for engineering systems.For dynamical systems, the most commonly used approaches include the Euler-Lagrange method, the Hamiltonian principle, and Newton's second law of motion.In this paper, we derive a mathematical model for the double variablelength pendulum with counterweight mass, which has been previously studied on a smaller scale [34].Specifically, we develop a more comprehensive model that captures the dynamics of the double variable-length pendulum with counterweight mass in greater detail.
Lagrange equations are commonly used to derive the equations of motion for double pendulum systems [10,35].For example, in [22], the Euler-Lagrange approach is used to derive the governing equations of a base-excited double pendulum with a magnet attached to the lower end of the second pendulum and equally spaced coils placed near the curvature of the magnet loop.This system can be used to generate electricity, and the harvested energy can be increased by increasing the amplitude and frequency of excitation, which can also be achieved by adding more coils to the harvester [22].Additionally, a double pendulum with one-sided rigid suppression under symmetric excitation leads to an asymmetric and chaotic system, as presented in [14].
When testing dynamic systems, it is important to consider friction since it significantly impacts the system's dynamic properties and operating conditions.The dynamic model presented in this work requires the consideration of friction phenomena between the disks and the cables on which the oscillating bodies are suspended.Therefore, we provide below a few significant publications related to this phenomenon: [2,5,20].The contribution of these publications to the description of the dry contact friction can be summarized as follows.Beginning from the first one, Kraus et al. focus on modeling and compensating for friction, with particular attention to its impact on the rotational motion of pulleys placed in series and on the control of parallel robots with cable drive [20].Incorporating friction into our model to account for dissipative forces, the dynamics of an extensible cable wound around a pulley, with attention to the Coulomb friction between the pendulum cable and its support, is comprehensively addressed in [5].In addition, Bazaei et al. analyzed the friction of limited arms moved at low speeds and developed a new model to compensate for this phenomenon [2].This model was verified through simulation and experimental studies that used a feedback control system.Some significant events associated with the friction modeling addressed in our study encompass pure dry sliding friction, the stick-slip effect, viscous friction, and potentially the Stribeck effect or frictional delay.These phenomena have been investigated in a survey of several friction force models for dynamic analysis of multibody mechanical systems by Marques et al. [26]; a continuous friction model for dynamics and control with physically meaningful parameters' by Brown and McPhee [7]; modeling friction in the dynamics analysis of selected one-degree-of-freedom spatial linkage mechanisms by Harlecki and Urbaś [16]; and also modeling, analysis, and control of dynamical systems with friction by Olejnik et al. [29].
With respect to the above, we have included pure dry sliding friction, stick-slip effect, viscous friction to capture phenomenological model of the nonlinear dynamical system to gain a better understanding and determine the precise reactions of the double pendulum of variable length.For instance, as presented by [26] and others, static models describe the steady-state behavior of relative frictional force or velocity, while dynamic models record more physical responses and properties by incorporating additional state variables [31].For this reason, it may be necessary to incorporate dynamic friction models with additional control parameters in the double variable-length pendulum with counterweight mass to capture the exact physical response of this complex mechanical system with unpredictable dynamics.
Electromagnets commonly use solenoids, which are wire coils arranged in a spring-like configuration.When an electric current flows through the solenoid, it generates a magnetic field.The strength of this field is determined by the properties of the solenoid and produces a force on charged particles that is directly proportional to their magnitude.A plunger and solenoid or a pair of solenoids with uniform winding make a typical tractive electromagnet pair in the investigated system.The solenoids experience forces from the current flow-ing through them, which causes the plunger to move.When the forces acting on the plunger are equalized, the plunger ceases to move.For instance, the plunger is centered in the solenoid and balances the forces.An electromagnet can quickly alter the magnetic field by adjusting the electric current flowing through the winding.A solenoid-generated external magnetic field interacts with a cylindrical permanent magnet plunger to drive a magneto-mechanical oscillator.For instance, this mechanism was used by Polczyński et al. to investigate various pendulum systems [33].
Hao et al. introduce a new magneto-mechanical oscillator model that emphasizes significant parallel leaf spring deflection and the electromagnetic actuator [15].A model of magnetic forces from a single solenoid actuator was developed using the charge model and infinitesimal element analysis.Findings showed that the electromagnetic force depends on the permanent magnet's position, hindering pure harmonic excitation for large oscillation amplitudes.An electromagnetic actuator design with paired solenoids to control double pendulum dynamics was introduced.This design employs quasi-constant force to maintain near-constant excitation.This actuator provides a broad excitation range.Motion of a magnetically forced spring pendulum was studied also using a uniform magnetic field [28].
Due to the increasing need to consider nonlinearities in dynamical systems, bifurcation diagrams have become a highly useful tool in engineering applications for the nonlinear analysis of dynamical systems.Therefore, the analysis of bifurcation phenomena observed in the dynamic response of analyzed systems cannot be avoided, especially in the case of nonlinear multiparameter dynamical systems.Bifurcation diagrams can be constructed simply using appropriate orthogonal decomposition modes that are computed on the fly when the bifurcation parameter is changed [42].
The theory of dynamical systems, in which differential equations play a fundamental role [17], is used to describe complex dynamic behavior.In such systems, period-doubling bifurcation occurs when a small change in a system parameter causes a new periodic trajectory to emerge from the existing periodic trajectory, with a doubled original oscillation period.It is important to distinguish between the halving bifurcation and the period-doubling cascade.This is the standard route to collapsing into chaos.This matter refers to homo-clinic, subharmonic, and superharmonic bifurcations for a periodically varying length pendulum [4].
Understanding the predictability of a dynamical system requires a definition of the concept itself, as well as the application of mathematical tools such as Lyapunov exponents.In this work, we focus on the spectrum of Lyapunov exponents, which enables us to characterize the rate of separation of the double pendulum mass from its initial position in space over time.The system's complete trajectory is determined by its position and momentum, and to analyze its chaotic orbits, we calculate the greatest Lyapunov exponent as described in [6].In general, a positive Lyapunov exponent indicates that the double variable-length pendulum with counterweight mass under study in this work is chaotic, provided that other conditions such as the compactness of the phase space are met [29].
In [9], Lyapunov exponents were investigated as a function of excitation amplitude during modulation of chaotic motion of a pendulum with excitation.In [24], the authors applied the Poincaré-Birkhoff fixedpoint theorem and showed that there is a fixed, nonpermanent periodic solution for a forced pendulum of variable length.In [37], it is indicated that by appropriately selecting the threshold speed during modulation of the rotational pendulum force, a controlled pendulum can achieve constant rotations for both constraining regularities, despite different initial conditions and forcing parameters.The concept of a variable-length pendulum has also been applied in [44,48] to control a crane system and in [12] to improve performance of nonlinear absorber using surrogate optimization.
Based on a literature analysis, it is clear that modeling and analyzing parametric dynamical models of variable-length pendulums with various forms of external forcing can be very challenging.However, the potential applications of these systems in mechanical and mechatronic systems make them an important subject of research, with significant relevance in both theory and practical engineering applications.Additionally, investigating these systems allows for the use of modern tools and methods for solving and analyzing parametric dynamical systems.
In Sect.2, we introduce our original experimental test stand.In Sect.3, we describe the physical and mathematical models used to capture the most important mechanical phenomena and general principles of the dynamical modelling process.These models include Lagrangian mechanics, friction modelling in pulley sliding bearings, electromagnetic actuation mechanisms, and kinematic forcing.Section 4 presents the results of our numerical simulations, nonlinear dynamical analyses, and vision-based measurements on the real experimental stand.Finally, in Sect.5, we provide our observations and conclusions.

The experimental laboratory stand
In this part of the paper, the main principles of construction and purpose of the original research setup are presented, aimed at achieving the desired dynamic functionalities.
In the investigated double variable-length pendulum with counterweight mass, one of the weights, called a pendulum, is set in pendulous motion by a kinematic forcing, while the other weight, a counterweight on the second end of the elastic cable joining both weights, can move freely in the direction of the Earth's gravitational field.This model assumes that both weights move in the same plane and do not interfere with each other.As already mentioned in previous chapters, the practical development undertaken in this work involves the attachment of a second pendulum, which is springily connected to the former.This concept was analyzed with the aim of implementing a real system with the above distribution of vibrating masses, subjected to kinematic forcing applied to the suspension point of the variable-length double pendulum.Given that the described pendulum is suspended by a rotating pulley and its length is variable, the design of the proposed stand takes into account a rigid frame structure equipped with a servo drive that sets the movement of the end-effector-a slider in a linear ball bearing with an attached pulley.
An original laboratory stand, shown in Fig. 1, allows for observation of the dynamical behavior of a real double variable-length pendulum with counterweight mass.Control of the horizontal movement of one of the rotating pulleys, which is placed in the slider, is achieved using a servomotor (2).This allows for limiting the movement of the counterweight of mass M, or for performing kinematic forcing of the suspension point of the flexible pendulum.These actions provide a limited dynamic range of vibration amplitude for both pendulum bodies, m 1 and m 2 , while keeping the system in a sensitive state of dynamic equilibrium.The pendulum's suspension was forced by replacing one pulley with two smaller pulleys.The axis of rotation of the pulleys is located above the pendulum and is not permanently connected to the frame of the station.Instead, it has the ability to move horizontally and relatively to the frame according to a given trajectory of movement, such as the function f 0 sin(ωt) of time t, where f 0 is the amplitude and ω is the angular frequency.Using a position-or velocity-controlled servo, the designed propulsion system is able to accurately reproduce the complex trajectories of the position and speed of the suspension point of the pendulum under test.
To block or partially restrict the movement of the counterweight, two current-fed coils are located near it.With this solution, the displacement of the counterweight with mass M can be periodically limited with different forces and alternating electromagnetic fields.
The stand is constructed in such a way that a double pendulum with masses m 1 and m 2 can move freely in a complex motion, with a pair of electromagnets periodically switched on or off to act on a counterweight of mass M. As a result, the length of the first pendulum does not increase beyond the permissible value determined by the length of the cable.The upper limit of the length of the first pendulum is reached when the sum of the masses m 1 + m 2 is much greater than the mass M of the counterweight or when the sum of the instantaneous contributions of the inertia force in the pendulum motion exceeds the inertia of the counterweight.In both situations, the pendulum system drops downwards, and the cable hits the ground.
Aluminum modular profiles provide the necessary construction and strength to meet the assumptions of the mathematical model and the physical problem defined based on the concept of the variable-length pendulum.The remaining components are made of machined aluminum or 3D printed plastic.Specific components of the experimental stand, designed to operate with a double pendulum slung through pulleys, are shown in the pictures below (see Fig. 2).
A connecting rod (9) with Igus spherical heads attached to the ends (8) was used as a link between the crank (7) and the end-effector, which consists of a linear plain bearing with an attached movable pulley (14).This solution allows us to adjust the length l 2 of the second arm of the mechanism.The ball joints of the Igus heads compensate for inaccuracies that may arise during installation of the drive system, from the motor Fig. 1 A prototype of the experimental realization of the laboratory stand is shown placed on a Stewart platform table.The pendulum masses, shown in white and blue, are fixed to the cable, and the slider is marked in green Fig. 2 A particular view of the experimental station's components is shown, highlighting the crank mechanism (a) and the overhanging pendulum cable.A fast camera uses a green marker to estimate the position of the slider axis to the end-effector.The linear bearing of the slider is achieved using a rail ( 12) and a ball bearing (11).
The stand is equipped with beam mounts (17) for the electromagnet coils (16), which are enclosed in plastic coil mounts (15) that allow for adjustment of their spacing and modification of the shape of the magnetic field acting on the counterweight M.These two coils, located above and below the counterweight, serve as controllable repellers that stop its position.Using the capabilities of the IndraDrive servo amplifier, we can control these coils from an external LabVIEW program.Thanks to this, this pair of coils acts as a generator of an external electromagnetic field, which forces or sometimes initiates the beginning of an experiment.

Mathematical and dynamical modeling
In this section, we outline the process of deriving the equations of motion for a four-degree-of-freedom mechanical system, which serves as a model for a swinging machine [46].We will employ Newton's second law of motion and the Lagrange-Euler method to establish the governing dynamic equations of motion for the double variable-length pendulum with a counterweight mass.Both approaches are anticipated to yield a precise mathematical representation of the analyzed dynamic systems.Prior to deriving the model, several assumptions were made: A1 A damped spring model between the two pendulum masses m 1 and m 2 is incorporated.A2 The mass of the cable that couples all the system masses is negligible.A3 The system is non-symmetric.This means that the kinematic excitation applied to the point of suspension of the two pendulum masses m 1 and m 2 puts more inertia on the right-hand side of the pendulum of the two coupled masses m 1 and m 2 , as shown in Fig. 3. A4 Air resistance-caused drag forces acting on the masses of the system are neglected.

The dynamical modelling
Our derivations in this section are based on the use of Newton's second law and the associated free body diagram of the investigated physical model shown in Fig. 3.All forces acting in the plane (x, y) and dimensions are indicated with respect to the origin O (0,0) .After reviewing the potential fields of application discussed in [46], we added a spring pendulum on the Fig. 3 Free body diagram of the investigated physical model of the double variable-length pendulum with counterweight mass and kinematic excitation opposite side of the counter mass M and a suspension system based on a Maxwell element with a stiffness k and a damper c placed between the two pendulums with masses m 1 and m 2 .Point O 1 is fixed, while O 2 is movable and can oscillate on the line (O, X ), allowing for the variation of the length l 1 (t)-the first state variable-and the double pendulum couplings.The length l 20 between the two point-focused masses of the pendulum is constant due to elongation under a static gravitational force, while l 2 (t) is the state measured after the dynamic extension caused by the spring connection with stiffness k between the masses.Using Newton's second law of dynamics, we established the equations of motion, assuming that the links are massless and the masses are concentrated at the points where the subsequent cords connect.The vector to the counter mass M in the assumed system of coordinates is defined as follows: where X 0 (t) = f 0 sin (ωt + θ) is a time function that represents the periodic kinematic excitation, measured from the origin of the global coordinate system O (0,0) in the X direction.Here, ω and θ are the angular frequency and phase shift of the excitation, s represents the distance in the X direction from O (0,0) to the fixed support point O 1 , while L n represents the length of the entire cable between the centers of the counterweight mass M and the pendulum m 1 .
The sum of forces acting on mass M in the Ydirection is given by: where T R represents the total friction force in both pulley bearings, including Coulomb and viscous friction; see Sect.3.3.By defining the second pendulum's length L 2 (t) = l 20s + l 2 (t), we introduce vectors for both pendulous bodies using the coordinates of vectors r i : where l 20s = l 20 + l 2st , with l 20 being the free length of the spring with stiffness k, and l 2st = m 2 g/k representing the static elongation of the spring.We now define the forces acting on the masses of the first and second pendulums, denoted m 1 and m 2 , respectively: Here, T t,3 = kl 2 (t) + c l2 (t) represents the force in the spring-damper Maxwell coupling-based model.The equations of motion are given in general form: where F i X and F iY are the force components acting on mass m i in the (X, Y )-plane.The Newtonian model is governed by the ODEs (5), which express the dynamic equilibrium of forces acting in the considered mechanical system.Based on the deductions above, the dynamics of the system's response with generalized coordinates l 1 (t), l 2 (t), ϕ 1 (t), and ϕ 2 (t) are described by four secondorder ordinary differential equations (6), which are as follows: In Eq. ( 6), ϕ 21 = ϕ 2 − ϕ 1 , Ẍ0 = −ω 2 f 0 sin ωt, f 0 is the maximum amplitude of the kinematic excitation, and ω is the excitation frequency.Furthermore, ϕ 1 (t) and ϕ 2 (t) represent the angles of deflection of the pendulum arms of length l 1 (t) and l 2 (t) (relative to their predecessors) from the vertical axis crossing the center of the movable pulley at O 2 , as shown in Fig. 3.

Lagrangian of the mechanical system
The system under investigation is structurally complex, so it should be examined from an energetic perspective.Ideally, one would find its Lagrangian and use the Euler-Lagrange method based on energy to derive the mathematical model.This method combines the kinetic and potential energies, denoted by T and U , respectively.In this section, the same symbols are used to derive the equations of motion.
When proposing the formula for potential energy, we assume that the masses in the system are concentrated in points and that the links connecting them are massless.A coupling between masses m 1 and m 2 is modeled using a Maxwell element with spring stiffness k and dashpot damping c.The time-dependent kinematic excitation is denoted by X 0 (t) = f 0 sin (ωt), which represents the distance from O 2 to the right in the X 2 -direction.The masses M and m 1 are hanging at Y 0 each and are separated by a distance of s + X 0 , so the length of the cable is L n = 2Y 0 + X 0 .To calculate the potential energy only, we take Y 0 as the reference coordinate and set it equal to −L n − f 0 + s.
The coordinates for the potential energies of all the masses are measured in the global coordinate system.A change in the length l 1 (t) of the first pendulum will affect only Y M , regardless of the angle ϕ 1 (t).
To calculate the first term of U for m 1 , we use the projection of l 1 (t) onto the Y -axis, i.e., l 1 (t) cos ϕ 1 (t).An increase in Y 1 from Y 0 also appears on the Y -axis.If Y 0 represents the equilibrium position for masses M and m 1 at ϕ 1 (t) = 0, then we introduce the distance from O (0,0) of the global system to m 2 , which is in equilibrium at ϕ 2 (t) = 0.If mass m 2 moves from Y 20 to Y 2 due to the influence of Y 1 , then virtual work is done in the gravitational field at some distance from Y 20 (at the free-hanging lengths of the second and first pendulum, subsequently at Y 0 , when M is also at Y 0 ) to Y 2 .An increase Y 1 affects the potential energy of m 2 .We have all the components needed to write the total potential energy for the spring of stiffness k between m 1 and m 2 .After the analysis above, we can propose the potential energy in the following form: where l 2 (t) represents an incremental change in the length of the second pendulum, measured from l 20 .
The total kinetic energy T as the sum of the kinetic energies T M , T 1 , T 2 possessed by the system masses working in a serial chain connection M-m 1 -m 2 is given in the form: Subsequently, the non-conservative force that captures damping of the system along the displacement l 2 is derived from the Rayleigh dissipation function: The Euler-Lagrange equations, with the Lagrangian L defined in terms of four degrees of freedom, namely the vector of generalized coordinates q i = [l 1 , l 2 , ϕ 1 , ϕ 2 ] and the vector of generalized forces Finally, by solving the Euler-Lagrange equations (10) with the assumptions and definitions given in ( 7), (8), and ( 9), we obtain an identical system of secondorder differential equations (6), which proves the correctness of our derivations.
In summary, the obtained model describes the dynamic behavior of the mechanical system under examination, which varies over time and throughout the available space of the angles of rotation of both mathematical pendulums and the entire length of the pendulum arms, including negative lengths.However, the negative length l 1 of the first arm is beyond the scope of this model, as it would cause the cable to lose contact with the second movable pulley.Thus, the numerical experiment, as well as the initial conditions and later dynamic states, must ensure that l 1 (t) > 0. This poses problems in finding an appropriate set of constant parameters and initial conditions imposed on the tested dynamical system.Additionally, dynamic analysis, including the determination of bifurcation plots, becomes much more challenging, as this condition must be met for all values of the control parameter and in each iteration of the discrete solution.

Modeling friction in pulley bearings
Now, we will present a model of the friction force in the bearings of both pulleys placed horizontally at a distance of X 0 (t) + s, as shown in Fig. 3.We consider Coulomb and viscous friction.In the first and most intuitive assumption, the rotation of the pulleys takes place under friction in the contact surface with the cable.It should be noted that the diameter of the pulley is disregarded, but it cannot be too small formally, because of the difficulty in describing the change in the tension force when bending the cable.
The Coulomb friction model is used to model the tangential forces on the contact surface, which depend on the relative speed of movement of the two dry surfaces of solids forming the frictional contact [29].Coulomb friction models are widely used, and depending on the application, a number of functional forms are distinguished, with a discontinuity expressed mathematically by the sgn function.However, we begin the current derivation by defining a simple frictional contact model as the starting point, which linearly connects the Coulomb and viscous friction models.
We derived the equations of motion for the considered triple pendulum (reduced to a double pendulum in the context of the dependent movement of the countermass M), which is suspended by an inextensible cable through two pulleys.We assumed dry friction occurring on the contact surfaces of both rotary slide bearings, while ignoring the friction on the contact surfaces of the cable and the two roller sleeves mounted on these bearings.To bring our dynamic model closer to the real system used in the experimental part of this work, we introduced the bearing resistance T R,i in the form of the sum of Coulomb and viscous friction forces.
We assumed that the cable connecting bodies M and m 1 experiences no slip on the outer surface of these sleeves.In the ideal case without friction, T t,i = T t,i−1 , so the tensile forces in front of and behind the rolling sleeves are equal.However, after taking into account the effects of the resistance on the rotation of the sliding bearings, i.e., Coulomb friction force T c,i and the force of viscous friction T v,i at the i-th pulley at points O i , we obtain: The above formula describes the balance of forces in the cable guided through rollers around both sliding bearings, taking into account the friction force T R,i at the cable speed l1 .
Viscous friction depends on the speed of movement of the cable, the associated sliding bearing, and a certain factor.Assuming that the cable is constantly stretched over the entire length of L n , which is true for the mechanical system under consideration, its speed around each bearing will be the first derivative of the change in length l 1 (t) of the pendulum of m 1 .In this case, it is sufficient to assume a certain constant value of the viscous friction coefficient, denoted by d i for each of the pulleys.Thus, we arrive at the formula The first component of the drag force, namely Coulomb friction, seems to be more difficult to model.The normal force exerted by pull forces on both sides of the sliding bearings must be determined.One of these forces acts at a constant wrapping angle α 1 of the roller, while the second one depends on the α 2 angle and changes with time due to the dynamic change in the angle ϕ 1 (t) of the m 1 pendulum.The wrapping angles for the first (fixed) and the second pulley moving horizontally with amplitude X 0 (t) are as follows: Equation 12 defines the Coulomb friction force in the bearing of the i-th pulley.This force is calculated based on the basic definition, where T c,i is the friction force, μ i is the coefficient of friction for the i-th pulley, T N ,i is the normal force on the pulley bearing, and l1 is the velocity of the cable.The sign function, sgn, ensures that the direction of the friction force opposes the direction of the cable's motion.
Following [20], T N ,i is the i-th force of normal reaction of a particular pulley bearing: The normal force in Eq. ( 13) acts at the bisector of the wrapping angle formed by the cable and the pulley.This allows for the following approximation: The relationship for the tension forces acting between successive pulleys is unknown.To determine this relationship, we first eliminate T N ,i by substituting Eq. ( 14) into (12) and then, substituting T c,i in (11).We then solve the resulting algebraic equation for T t,i , giving: for Tc,i = sgn The additional term that reflects the total frictional resistance T R , which is an external force to the mechanical system and exists in both pulley bearings O i , can be found as a sum of differences T R,i = T t,i − T t,i−1 (see Eq.( 11)).Specifically, T R = T R,1 + T R,2 , where the first tension force is initially T t,0 = Mg − m 1 l1 (t) (see Fig. 3).

The electromagnetic actuator design and output
Figure 4 shows the details of a double-identical solenoid actuator.This type of actuator uses two solenoids that are mounted in parallel with each other [11,13].The use of two identical solenoids offers several benefits.First, it ensures that the actuator produces a balanced force because the two solenoids are identical in size and strength.Second, it provides redundancy in case one of the solenoids fails, as the other solenoid can continue to provide motion [8,11,13,50].Finally, it allows for greater control over the motion produced by the actuator since the two solenoids can be independently controlled, as discussed in [11] or [33].
As stated by [11], the charge model states that the force exerted by a permanent magnet in an external magnetic field can be expressed as follows: where ρ = −∇ • M and σ = M • n are equivalent volume and surface charge density, respectively; B ext is the external magnetic field, and M is the magnetic moment of the permanent magnet per unit volume.Additionally, the formula M = B r /μ 0 can be used to determine the magnetic moment, where B r is the residual flux density of the permanent magnet.The cylindrical magnet plunger is polarized with fixed and uniform magnetization along its axis, denoted as M = M s i, resulting in ρ = −∇ • M = 0.The magnetic force acting along the x-axis is given by: where B x is a function of the variable x.In our case, the magnetic plunger is a hollow cylinder, then, to evaluate Fig. 4 The magnetic plunger rod is a counterweight placed between two identical solenoids in the twin-solenoid actuator σ , the unit surface normal is determined as follows: where e r is the unit vector in the direction of r .If the permanent magnet's north and south poles are on its left and right sides, respectively, then its end surface's charge density is described as: where l, d 1 , and d 2 are the lengths of the inner and outer radii of the magnet cylinder, x l and x r are the coordinates of the permanent magnet's left and right regions on the x-axis.Equations ( 16) and ( 19) are used to calculate the force acting on the permanent magnet plunger, as described below: The amplification of the driving force is a function of the displacement variance x, as shown by the calculations above for a single solenoid actuator.A more effective electromagnetic actuator is described in [15], which involves two similar and coaxial solenoids with coils wound in opposite directions and fed by the same currents.This results in a constant force that is independent of displacement, making it suitable for mechanical and mechatronic applications.
We assume that the distance between the centers O and O of the wound air coils that create the pair of solenoids passing the same, but mutually reversed direction current is d.The magnetic moving body, with a coordinate system attached at O, operates approximately along the imaginary line created by O − O .
Using Eq.( 20), we can derive the electromagnetic force as shown below: To analyze the nature of the resulting electromagnetic force acting on the counterweight mass, or plunger (as seen in Fig. 4b), in the experimental stand, force curves were constructed and are shown in Fig. 5a.The characteristics of the electromagnetic force generated by a constant direct current, such as those in Fig. 5a, can be used to determine which of the four types of quasi-constant force (QCF), oblique linear force (OLF), single peak force (SPF), and double peak force (DPF) is being produced.
The various types of forces reflect the connections between displacement and force, which can be altered by changing the properties of the electromagnetic actuator, such as the length, thickness, and radius of the solenoid, the number of turns in the coil, and the distance between a pair of solenoids.As d decreases and increases, the force curves transform into SPFs and DPFs, respectively, for instance, at d = 75 and 92 mm .For example, according to Fig. 5a, the QCF zone is present at x = 0 when d = 83 mm .The OLF region is visible in the plot shown in Fig. 5a, where the force is 123 Fig. 5 The actuator's force characteristic, including the QCF, SPF and DPF (a) and its dependence on varying current in the coils (I = 3 A for experiments in Sect.4.3) Fig. 6 The segments of OLF and SPF for the single solenoid actuator.The single solenoid actuator's forces at various parameters d 2 and L between −10 and 15 N .The segments of OLF and SPF for the single solenoid actuator can be seen in Fig. 6a  and b.It should be noted that a single solenoid actuator can only produce SPFs and OLFs.On the other hand, an actuator powered by a pair of identical solenoids can produce all four types of forces.

The kinematic forcing
To maintain the suspension point of the double variablelength pendulum, we used a servo connected to a movable pulley through a crank mechanism (as demonstrated in a similar way in [21]).As shown in Fig. 7a, the crank with a length of l 1 on the segment (x 0 , y 0 )-(x 1 , y 1 ) and the connecting rod with a length of l 2 on the segment (x 1 , y 1 )-(x 2 , y 2 ) form a simple kinematic chain.The end of this chain, located at the point (x 2 , y 2 ), determines the position of the end-effector, which is the slide that supports the suspended pulley.
Therefore, the relationship between the angular position α of the servo shaft and the position x 2 of the pulley in space is generally nonlinear.This means, for example, that the horizontal sinusoidal motion of the end-effector may not correspond exactly to the sine of the servo's α angle.
In the kinematic model of the driving system, we assume that the servo's axis of rotation is located at point (x 0 , y 0 ) in the coordinate system (X, Y ).The position of the servo shaft in this stationary reference frame is determined by the angle α.The point (x 1 , y 1 ) represents the joint connecting the crank (fixed-length arm l 1 ) and the connecting rod (fixed-length arm l 2 ) of the mechanism.The position of the end-effector of this mechanism on the plane is indicated by point (x 2 , y 2 ), which can only move along the X axis.Therefore, its coordinate in the Y direction is a constant value that depends on the dimensions of the frame structure.
When determining the desired functional relationship between the time dependent α angle and the desired position (x 2 (t), y 2 ), i.e., when solving the problem of inverse kinematics, we use Fig. 7a to find: The coordinates of the connecting rod position should be independent of the β angle and the coordinates of the intermediate joint at the point (x 1 , y 1 ).Another relation for y 1 is given by the formula: After inserting y 1 and calculating sin β, we obtain a relationship between angles α and β based on two rectangular triangles with hypotenuses l 1 and l 2 , respectively, as well as the second dependence of β angle on the observation x 1 = x 0 − l 1 sin α, i.e.: We use the Pythagorean identity sin 2 β +cos 2 β = 1 to incorporate the terms in Eq. ( 22) into an equation that includes the desired function of the α angle for an input function of displacement X 0 (t) of the endeffector with amplitude f 0 and angular frequency ω.We also take into account the approximate relationship between x 2 (t) and X 0 (t).
Continuing the transformation, we obtain the final nonlinear algebraic equation for the sought t-dependent α angle that the servo will produce: The problem of inverse kinematics, which consists of solving a nonlinear algebraic equation (23) in which time t is a parameter, is solved numerically in an approximate manner.The values of the α(t) function are obtained in the form of a series [t i , α i ], where i = 1 . . .n, subject to the constraint x 2 (t) = f 0 sin(ωt).
In Fig. 7, the green circle corresponds to the circle located at O 2 in Fig. 3a.The kinematic forcing trajectory α(t) is applied by the servo and shown on a background of a π -phase shifted sine.The function β(t) represents the angle of the connecting rod and is shown in green with an approximate analytical interpolation (b).
The trajectories of the numerical solution are illustrated in Fig. 7b and exhibit three periods of the α(t) function, measured in radians and obtained as the solution of Eq. ( 23) for n = 900 iterations at t ≈ 0.042 s.The parameters of x 2 (t) are ω = 0.5 rad/s , f 0 = 0.1 m , and the crank mechanism's lengths are l 1 = 0.2 m and l 2 = 0.4 m , with an initial position of (x 0 , y 0 ) = (−0.1,0) m.
Based on the analysis of the obtained solution, it can be seen that, with the construction parameters used on the test stand (l 1 = 0.16, l 2 = 0.18, y 0 = 0.08 m , and α(0) = 0 rad), the function describing the angular position of the servo shaft over time almost coincides with an ideal sine function.This analysis of the kinematics of the end-effector's excitation mechanism led to the assumption that the quasi-sinusoidal profile set on the servo will generate, with good approximation, a displacement of the form x 2 (t) = f 0 sin(ωt).Since the profile has a similar shape to a harmonic function, we assume that the trajectory of motion in the X direction of the movable pulley is described by a sine function.
We assume that y 2 is a constant value and express the horizontal displacement of the pulley that represents the reciprocating motion x 2 (t) of the suspension of the double pendulum by a nonlinear function that describes the relationship between: (a) the position of the endeffector and the α angle at which the servo drive shaft is set; (b) constant values, i.e., the arm lengths l 1 and l 2 ; (c) the position of the motor axis at a point determined by the dimensions of the frame structure (x 0 , y 0 ); (d) and the coordinate of the axis of rotation of the endeffector's pulley at y 2 , measured in the Y direction.

Results and discussion
This section unveils the outcomes derived from both our numerical and experimental simulations of the multibody nonlinear system.In order to ascertain the system's complexity and non-smooth characteristics, various techniques have been employed, encompassing analysis of time histories, phase planes, bifurcation diagrams, and Lyapunov exponents.The assessment extends further to the application of vision-based measurements in the final stage.This innovative approach enables the capture of motion trajectories of the double pendulum, thus enriching our understanding of the system's dynamic complexity.
The mathematical model given by Eqs. ( 6) and ( 12) has been solved numerically using the Runge-Kutta procedure with adaptive step-size.The results describe the concept of a double variable-length pendulum with counterweight mass, kinematic excitation and electro- magnetic forcing as well, demonstrating the system behavior with a double pendulum and suspension system between the two pendulums excited by an electromagnet and armature.The excitation frequency has a crucial influence on the system dynamics; therefore, its impact will be discussed in next step.

Time trajectories
Figure 8 shows the numerically computed time responses of a variable-length pendulum subjected to external kinematic forcing at the suspension point.
In Fig. 8a, three interesting frictional responses of the damped oscillator are shown for different static coefficients of friction on the sliding surface between the cords strung over the pulleys.The green curve clearly illustrates the influence of this parameter, exhibiting a coexistence of stick-slip friction that results in temporary absence of length variation in the pendulum.The pendulum starts its evolution at a length of 0.49 m , then stabilizes its oscillations around the intended target length, with lower first harmonic frequency during transient and higher second harmonic frequency in steady state.An intriguing observation is the decreasing amplitude range with increasing static friction force, leading to progressively decreasing peak responses in blue, red, and green colors (as the friction coefficient increases).After a considerable time interval, the variable-length pendulum reaches the target length, and small amplitude oscillations are observed around this length, associated with the forcing frequency at the suspension point.
On the other hand, the trajectories in Fig. 8b look very similar, but with the fundamental difference that the influence of the viscous friction coefficient in the governing equation ( 15) is no longer significant.This is evidenced by the nearly overlapping time characteristics.The analysis of bifurcations in the eight-dimensional state space of the dynamical system is based on variation of the angular frequency of kinematic excitation of point O 2 , denoted by ω.Bifurcation diagram of the orbits of the investigated dynamical system shown in Fig. 9 reveals numerous periodic windows.Such windows result from multiple period doublings of the vibrations of the first pendulum (in the direction of l 1 ), which frequently evolve into multi-periodic and occasionally into chaotic orbits.
It can be observed that the transition from the 1period orbit ω 19 through successive period doublings, and sometimes reductions in the periodicity number, is quite extensive.Only in the final range of bifurcation parameters (approximately from ω 210 onwards), a substantial region of chaotic vibrations emerges, as evidenced by the high value of the first (largest) Lyapunov exponent.The Kaplan-Yorke dimension d K P in the range ω 210...250 falls within 1.72−3.33.The high complexity of this diagram indicates that the dynamic response of the studied system is unpredictable and strongly dependent on the angular frequency of external excitation.
Conducting a deeper dynamic analysis of this bifurcation diagram, let us discuss a few interesting scenarios of dynamic changes shown in Figs. 10 and 11.The left part of the discussed bifurcation diagram of the investigated multi-degree-of-freedom dynamic system determines a regular and periodic pattern of vibration dynamics.This is evident, for example, through the Poincaré maps and trajectories displayed in Fig. 10 for 2000 periods of the external excitation X 0 (t).
A chaotic orbit emerged as the bifurcation parameter increased, following a series of wraps around the quasi-periodic whale tail loop shown in Fig. 10a and subsequent splitting at growing values of the kinematic excitation frequency is visible in Fig. 10b.This presents a highly interesting instance of a transition to chaos observed within this dynamic system (it is noteworthy that the Coulomb friction coefficient μ c is very small-see model parameters under Fig. 9).However, we are unable to illustrate this winding and subsequent complex quasi-periodic loop splitting here, as it would require including numerous figures depicting this transition.The sudden emergence of multi-period motion within the domain of chaotic vibrations is observed within the chaotic window at ω 228 = 8.899 radians per second, as depicted in Fig. 11a.By either increasing or decreasing the excitation angular frequency by an increment ω on the diagram, we return to irregular behavior and even to chaos, as shown in Fig. 11b.
To sum up the numerical experiments, conditioning and precise solution at high accuracy of numerical integration is required.Only feasible solutions determined by mechanical constraints should be considered for the double variable-length spring pendulum.The length l 1 of the first pendulum must be within the permitted length range, which is determined by the length L of the cable connecting the pendulum mass m 1 and the counter mass M, minus the distance s+X 0 (t).Here, s is the constant separation between both supports of the system, and X 0 (t) represents the sinusoidal kinematic forcing of the amplitude f 0 .
If the solution along l 1 (t) = 0, i.e., it reaches the left boundary or L − s − X 0 (t), the right boundary of the constraint, then numerical integration is stopped, as our physical model does not account for hitting an obstacle.Consequently, certain maps are omitted in all bifurcation diagrams, which can sometimes result in gaps for certain bifurcation parameters, creating an interrupted bifurcation diagram that reflects real dynamical effects.
The qualitative assessment of the obtained trajectories was performed using the spectrum of Lyapunov exponents, the maximal value λ 1 of which is shown in Fig. 9.The conditions under which the spectrum of Lyapunov exponents for the analyzed dynamical system has been computed are as follows: i) The QRdecomposition on the parallelogram matrix is applied N times.This value is computed as 10 • h • N e , where h represents the step time of numerical integration and N e represents the length of the eight-dimensional time trajectory discretized at every h.ii) The growth rates along the trajectory are then averaged over N successive steps, yielding the Lyapunov exponent spectrum.The parallelogram is re-normalized at each step.iii) The resulting spectrum, which has a vector size dependent on the dimension of the state space, is sorted from the maximum to the minimum value.iv) The initial state u 0 for computation is determined by initializing the numerical model from the same initial conditions used to solve the system.v) An extra transient time t tr is applied to evolve the system before the algorithm is applied.vi) An individual evolution time of t = 0.1 is used.

Vision-based measurements on the experimental stand
To measure the real displacements of the pendulum bodies and the slider on the laboratory stand, a visionbased method was employed.A video was recorded at a frame rate of 59.94 fps, with markers of different colors attached to the elements to be tracked to facilitate data processing, see Fig. 14.Subsequently, a Python script utilizing OpenCV was developed to extract the data.Three sinusoidal excitation profiles with different amplitudes were used as presented in Tab. 1, and videos were recorded for each profile.
Recording and processing with the picture frames resulted in a total of 12 time histories of real responses, as shown in Figs. 12 and 13.It is noticeable that for cases where the electromagnetic coils were switched off, the system begins to stabilize around t = 15 seconds.
In Fig. 12e, the pendulum even stopped falling completely, despite the 18 mm amplitude of the excitation signal applied to the coils.However, this excitation resulted in larger values of the angular displacements ϕ 1 and ϕ 2 , and caused the distances l 1 to follow a smoother trajectory with lower amplitudes.In the case of the bodies m 1 and m 2 coupled with springs, shown in Fig. 13, the second type of excitation resulted in system stabilization after t = 15 seconds .Case III confirms the numerical observations in Fig. 8 of slow stabilization of the lower frequency time response at a steady-state value, accompanied by oscillations of the length l 1 at the excitation frequency of the right cable-suspended roller, which serves as the pivot point for the entire pendulum system.
It is worth mentioning another interesting behavior of the first pendulum, observed in Case V. Despite the kinematic excitation of the first pendulum, it exhibits almost no rotation relative to ϕ 1 .Instead, due to the synchronization of vibrations with the second pendulum, the first flexible link moves left and right, behaving like a rigid connector and almost not bending.This confirms the unprecedented diversity of possible combinations of dynamic responses of the studied system.
During experimental investigations, several challenges were encountered: 1) Achieving precise mechanical excitation control for the movable roller's position proved difficult.After several unsuccessful attempts, a ready-made servo mechanism with a dedicated controller and regulator was implemented.This solution enabled the imposition of nearly any displacement characteristic for the suspension point while maintaining high positioning precision, regardless of friction resistance in the closely fitted sliding guide.2) Contactless measurement of each pendulum mass position required a high-precision video camera, followed by the development of image recognition algorithms, specifically those accurately detecting the centers of masses.Notably, this was challenging due to the bodies rotating around the link, exhibiting variable geometry concerning the object, and undergoing color changes.Filtering and transformations on all registered image frames allowed measurement independence from these factors.3) The last challenge involved estimating the center of the extracted shapes and calculating the observed shape's center in the adopted coordinate system.Real positions of observed bodies were detected after two-dimensional interpolation and sharpening with adaptation using a dedicated program.

Conclusions
For the first time, a research stand with a broad cognitive horizon based on the double variable-length pendulum with counterweight mass at various sources of forcing has been presented.The investigated dynamic system is adapted to observe a wide range of dynamic behaviors related to variable-length pendulum models and their modifications.Several ways of kinematic and force excitation were considered to control the dynamic behavior as well as to observe irregular dynamics, which are the basis of difficulties encountered when controlling such systems, such as elastically suspended weights on cables driven by motors mounted on arms or moving cranes.Using the detailed physical, mechanical, and mathematical foundations provided, this system can be developed in many directions, such as controlling the angle of the first pendulum deflection when disturbed by the second pendulum by attaching a lowmass spring.We are dealing with a strongly nonlinear dynamic system that includes kinematic excitation of the cable suspension point, excitation caused by electromagnetic interaction, sliding friction, and precise measurement of the moving mass markers using a vision system supported by highly specialized software capable of measuring the position of multiple observed markers with an accuracy of about one micrometer.
Real and numerical trajectories primarily demonstrate and confirm the significantly low frequency of length variation in comparison with the higher frequencies of pendulum body oscillations.In both cases, there is an observed gradual descent of the elastically linked masses of the investigated pendulum until the system attains a stable dynamic equilibrium.Subsequently, with the specified counterweight mass, the vibrations evolve into a stable periodic orbit.The presentation of numerical solutions serves to verify the correctness of deriving the intricate dynamic model and to showcase the potential dynamic complexity inherent in the studied system.This work lays the foundation for future explorations, including deep learning applications to unveil hidden physical phenomena.

Fig. 7 A
Fig. 7 A physical model of the crank mechanism (a) and its time trajectory (b) kinematically forcing the suspension point O 2 (x 2 , y 2 )

Fig. 10 A
Fig. 10 A whale-shaped blue Poincaré map (a) on the orange background of phase-space projection on the plane φ1 (ϕ 1 ) and its transformation into chaotic undetermined map (b).The typical quasi-periodic orbit (a) of time-dependent angle of rotation of the first pendulum is observed in the window ω 1...210 mostly occupied by periodic dynamics (t 0 and t f is the initial and final time of observation

Fig. 11 A
Fig. 11 A scenario of sudden appearance of a quasi-periodic solution (a) in the chaotic window for ω 210...250 situated on the right-hand side of the bifurcation orbit diagram in Fig. 9.A representative blue Poincaré map (b) on the background of orange phase-space projection l1 (l 1 ) of time-varying length of the first pendulum is observed in the window mostly occupied by chaotic dynamics